From 6c2e3e8a897b18c00c7315e6855b9caa9a389b2b Mon Sep 17 00:00:00 2001
From: CookieKastanie <jerem0205@yahoo.fr>
Date: Thu, 10 Mar 2022 22:09:21 +0100
Subject: [PATCH] fix

---
 src/project_surface_points_on_poly.cc | 31 ++++++++++++++++++---------
 1 file changed, 21 insertions(+), 10 deletions(-)

diff --git a/src/project_surface_points_on_poly.cc b/src/project_surface_points_on_poly.cc
index 971d292..9e6ebca 100644
--- a/src/project_surface_points_on_poly.cc
+++ b/src/project_surface_points_on_poly.cc
@@ -82,17 +82,24 @@ static void moveSurfacePoint(vtkUnstructuredGrid *tetMesh,
 	tetMesh->GetPoint(pointId, point);
 	double direction[3];
 	double distance;
+
+	double radiusScale = 5.;
+
 	closestPolyMeshPoint(polyMesh, point, kdTree, links, direction, &distance);
 	vtkNew<vtkIdList> affectedPoints;
-	kdTree->FindPointsWithinRadius(distance * 2, point, affectedPoints);
-	for (vtkIdType j = 0; j < affectedPoints->GetNumberOfIds(); j++) {
+	kdTree->FindPointsWithinRadius(distance * radiusScale, point, affectedPoints);
+	
+	if(distance > 0) 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);
+		if(dist2 > (distance * radiusScale) * (distance * radiusScale)) continue;
+
 		double motion[3] = {direction[0], direction[1], direction[2]};
-		double factor = 1 - std::sqrt(dist2) / (distance * 2);
+		double factor = 1. - std::sqrt(dist2) / (distance * radiusScale);
 		vtkMath::MultiplyScalar(motion, factor);
+		//std::cout << std::sqrt(dist2) << " " << (distance * radiusScale) << " " << factor << "\n";
 		motionVectors[affectedPointId * 3 + 0] += motion[0];
 		motionVectors[affectedPointId * 3 + 1] += motion[1];
 		motionVectors[affectedPointId * 3 + 2] += motion[2];
@@ -143,13 +150,17 @@ vtkTypeBool ProjectSurfacePointsOnPoly::RequestData(
 			if (i == surfacePoints->GetValue(j)) skip = true;
 		}
 		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();
+		
+		if(numberOfAffectors[i] != 0) {
+			double point[3];
+			tetMesh->GetPoint(i, point);
+			double *motion = motionVectors.data() + (i * 3);
+
+			vtkMath::MultiplyScalar(motion, 1. / numberOfAffectors[i]);
+			vtkMath::Subtract(point, motion, point);
+			tetMesh->GetPoints()->SetPoint(i, point);
+			tetMesh->Modified();
+		}
 	}
 
 	return true;