start implementing proportional motion

src/project_surface_points_on_poly.cc

@@ -12,6 +12,7 @@
 #include <vtkPolyData.h>
 #include <vtkKdTree.h>
 #include <vtkStaticCellLinks.h>
+#include <vtkDoubleArray.h>
 
 #include <limits>
 
@@ -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<double>::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<double> &motionVectors,
+                             std::vector<unsigned> &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<vtkIdList> 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<double> motionVectors(3 * tetMesh->GetNumberOfPoints(), 0);
+	std::vector<unsigned> numberOfAffectors(tetMesh->GetNumberOfPoints(), 0);
 	for (vtkIdType i = 0; i < surfacePoints->GetNumberOfTuples(); i++) {
-		double p[3];
+		moveSurfacePoint(tetMesh, polyMesh,
-		tetMesh->GetPoint(surfacePoints->GetTuple1(i), p);
+		                 surfacePoints->GetValue(i),
-		double dist2;
+		                 motionVectors,
-		vtkIdType closest = kdTree->FindClosestPoint(p, dist2);
+		                 numberOfAffectors,
-		vtkIdType *cellIds = links->GetCells(closest);
+		                 kdTree, links);
-		vtkIdType nCells = links->GetNcells(closest);
+	}
-		double min_dist = std::numeric_limits<double>::infinity();
+	for (vtkIdType i = 0; i < tetMesh->GetNumberOfPoints(); i++) {
-		double min_direction[3];
+		bool skip = false;
-
+		for (vtkIdType j = 0; j < surfacePoints->GetNumberOfTuples(); j++) {
-		/* For each tri the closest point belongs to. */
+			if (i == surfacePoints->GetValue(j)) skip = true;
-		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];
-			}
 		}
-
+		if (skip) continue;
-		vtkMath::MultiplyScalar(min_direction, min_dist);
+		double point[3];
-		vtkMath::Subtract(p, min_direction, p);
+		tetMesh->GetPoint(i, point);
-		tetMesh->GetPoints()->SetPoint(surfacePoints->GetTuple1(i), p);
+		double *motion = motionVectors.data() + i;
+		vtkMath::MultiplyScalar(motion, 1. / numberOfAffectors[i]);
+		vtkMath::Subtract(point, motion, point);
+		tetMesh->GetPoints()->SetPoint(i, point);
 		tetMesh->Modified();
 	}
 

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;
 };