From 56b426948c3780bdc4f225a80510638167704a7f Mon Sep 17 00:00:00 2001 From: papush! Date: Thu, 10 Mar 2022 16:51:05 +0100 Subject: [PATCH] start implementing proportional motion --- src/project_surface_points_on_poly.cc | 113 ++++++++++++++++++++------ src/project_surface_points_on_poly.h | 3 + 2 files changed, 91 insertions(+), 25 deletions(-) diff --git a/src/project_surface_points_on_poly.cc b/src/project_surface_points_on_poly.cc index 5375bda..971d292 100644 --- a/src/project_surface_points_on_poly.cc +++ b/src/project_surface_points_on_poly.cc @@ -12,6 +12,7 @@ #include #include #include +#include #include @@ -36,6 +37,74 @@ int ProjectSurfacePointsOnPoly::FillInputPortInformation(int port, vtkInformatio } +static void closestPolyMeshPoint(vtkPolyData *polyMesh, + const double *point, + vtkKdTree *kdTree, + vtkStaticCellLinks *links, + double *toClosestPoint, + double *distanceToClosestPoint) { + vtkIdType closestPolyMeshPoint; + { + double dummy; + closestPolyMeshPoint = kdTree->FindClosestPoint((double *) point, dummy); + } + vtkIdType *cellIds = links->GetCells(closestPolyMeshPoint); + vtkIdType nCells = links->GetNcells(closestPolyMeshPoint); + double minDist = std::numeric_limits::infinity(); + + /* For each tri the closest point belongs to. */ + for (vtkIdType i = 0; i < nCells; i++) { + vtkIdType cellId = cellIds[i]; + double direction[3]; + double dist = pointTriangleDistance( + (double *) point, polyMesh->GetCell(cellId), direction); + /* Find the closest one. */ + if (dist < minDist) { + minDist = dist; + toClosestPoint[0] = direction[0]; + toClosestPoint[1] = direction[1]; + toClosestPoint[2] = direction[2]; + } + } + vtkMath::MultiplyScalar(toClosestPoint, minDist); + *distanceToClosestPoint = minDist; +} + + +static void moveSurfacePoint(vtkUnstructuredGrid *tetMesh, + vtkPolyData *polyMesh, + vtkIdType pointId, + std::vector &motionVectors, + std::vector &numberOfAffectors, + vtkKdTree *kdTree, + vtkStaticCellLinks *links) { + double point[3]; + tetMesh->GetPoint(pointId, point); + double direction[3]; + double distance; + closestPolyMeshPoint(polyMesh, point, kdTree, links, direction, &distance); + vtkNew affectedPoints; + kdTree->FindPointsWithinRadius(distance * 2, point, affectedPoints); + for (vtkIdType j = 0; j < affectedPoints->GetNumberOfIds(); j++) { + vtkIdType affectedPointId = affectedPoints->GetId(j); + double affectedPoint[3]; + tetMesh->GetPoint(affectedPointId, affectedPoint); + double dist2 = vtkMath::Distance2BetweenPoints(affectedPoint, point); + double motion[3] = {direction[0], direction[1], direction[2]}; + double factor = 1 - std::sqrt(dist2) / (distance * 2); + vtkMath::MultiplyScalar(motion, factor); + motionVectors[affectedPointId * 3 + 0] += motion[0]; + motionVectors[affectedPointId * 3 + 1] += motion[1]; + motionVectors[affectedPointId * 3 + 2] += motion[2]; + numberOfAffectors[affectedPointId]++; + } + vtkMath::Subtract(point, direction, point); + tetMesh->GetPoints()->SetPoint(pointId, point); + tetMesh->Modified(); + affectedPoints->Reset(); +} + + vtkTypeBool ProjectSurfacePointsOnPoly::RequestData( vtkInformation *request, vtkInformationVector **inputVector, @@ -59,33 +128,27 @@ vtkTypeBool ProjectSurfacePointsOnPoly::RequestData( vtkIdTypeArray *surfacePoints = vtkIdTypeArray::SafeDownCast(tetMesh->GetFieldData()->GetArray("surface_points")); + std::vector motionVectors(3 * tetMesh->GetNumberOfPoints(), 0); + std::vector numberOfAffectors(tetMesh->GetNumberOfPoints(), 0); for (vtkIdType i = 0; i < surfacePoints->GetNumberOfTuples(); i++) { - double p[3]; - tetMesh->GetPoint(surfacePoints->GetTuple1(i), p); - double dist2; - vtkIdType closest = kdTree->FindClosestPoint(p, dist2); - vtkIdType *cellIds = links->GetCells(closest); - vtkIdType nCells = links->GetNcells(closest); - double min_dist = std::numeric_limits::infinity(); - double min_direction[3]; - - /* For each tri the closest point belongs to. */ - for (vtkIdType j = 0; j < nCells; j++) { - vtkIdType cellId = cellIds[j]; - double direction[3]; - double dist = pointTriangleDistance(p, polyMesh->GetCell(cellId), direction); - /* Find the closest one. */ - if (dist < min_dist) { - min_dist = dist; - min_direction[0] = direction[0]; - min_direction[1] = direction[1]; - min_direction[2] = direction[2]; - } + moveSurfacePoint(tetMesh, polyMesh, + surfacePoints->GetValue(i), + motionVectors, + numberOfAffectors, + kdTree, links); + } + for (vtkIdType i = 0; i < tetMesh->GetNumberOfPoints(); i++) { + bool skip = false; + for (vtkIdType j = 0; j < surfacePoints->GetNumberOfTuples(); j++) { + if (i == surfacePoints->GetValue(j)) skip = true; } - - vtkMath::MultiplyScalar(min_direction, min_dist); - vtkMath::Subtract(p, min_direction, p); - tetMesh->GetPoints()->SetPoint(surfacePoints->GetTuple1(i), p); + if (skip) continue; + double point[3]; + tetMesh->GetPoint(i, point); + double *motion = motionVectors.data() + i; + vtkMath::MultiplyScalar(motion, 1. / numberOfAffectors[i]); + vtkMath::Subtract(point, motion, point); + tetMesh->GetPoints()->SetPoint(i, point); tetMesh->Modified(); } diff --git a/src/project_surface_points_on_poly.h b/src/project_surface_points_on_poly.h index 926a562..9433f96 100644 --- a/src/project_surface_points_on_poly.h +++ b/src/project_surface_points_on_poly.h @@ -14,8 +14,11 @@ public: vtkTypeBool RequestData(vtkInformation *request, vtkInformationVector **inputVector, vtkInformationVector *outputVector) override; + vtkSetMacro(affectedNeighborsCount, int); + vtkGetMacro(affectedNeighborsCount, int); protected: + int affectedNeighborsCount = 0; ~ProjectSurfacePointsOnPoly() override = default; };