From 8ee09d686e6cbc3098428cd37f516675f497d2a5 Mon Sep 17 00:00:00 2001
From: papush! <contact@ccolin.fr>
Date: Thu, 17 Mar 2022 15:18:19 +0100
Subject: [PATCH] fix surface point projection kdtree logic

but still use the brute-force approach because it still doesn't work
---
 src/project_surface_points_on_poly.cc | 24 +++++++++++++++++-------
 src/project_surface_points_on_poly.h  |  3 ++-
 2 files changed, 19 insertions(+), 8 deletions(-)

diff --git a/src/project_surface_points_on_poly.cc b/src/project_surface_points_on_poly.cc
index 8a9e38e..0e9891f 100644
--- a/src/project_surface_points_on_poly.cc
+++ b/src/project_surface_points_on_poly.cc
@@ -88,24 +88,31 @@ void ProjectSurfacePointsOnPoly::moveSurfacePoint(vtkUnstructuredGrid *tetMesh,
                              vtkIdType pointId,
                              std::vector<double> &motionVectors,
                              std::vector<unsigned> &numberOfAffectors,
-                             vtkKdTree *kdTree,
+                             vtkKdTree *polyMeshKdTree,
+                             vtkKdTree *tetMeshKdTree,
                              vtkStaticCellLinks *links) {
 	double point[3];
 	tetMesh->GetPoint(pointId, point);
 	double direction[3];
 	double distance;
 
-	closestPolyMeshPoint(polyMesh, point, kdTree, links, direction, &distance);
+	closestPolyMeshPoint(
+		polyMesh, point, polyMeshKdTree, links, direction, &distance);
 	vtkNew<vtkIdList> affectedPoints;
 	findPointsWithinRadius(distance * RadiusScale, tetMesh, point, affectedPoints);
-	// kdTree->FindPointsWithinRadius(distance * radiusScale, point, affectedPoints);
+	// tetMeshKdTree->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 expected = (distance * RadiusScale) * (distance * RadiusScale);
+		if (dist2 > 2 * expected) {
+			std::cerr << dist2 << " " << expected << "\n";
+			// throw std::runtime_error("");
+			continue;
+		}
 
 		double motion[3] = {direction[0], direction[1], direction[2]};
 		double factor = 1. - std::sqrt(dist2) / (distance * RadiusScale);
@@ -164,8 +171,10 @@ vtkTypeBool ProjectSurfacePointsOnPoly::RequestData(
 	vtkNew<vtkStaticCellLinks> links;
 	links->BuildLinks(polyMesh);
 
-	vtkNew<vtkKdTree> kdTree;
-	kdTree->BuildLocatorFromPoints(polyMesh->GetPoints());
+	vtkNew<vtkKdTree> polyMeshKdTree;
+	polyMeshKdTree->BuildLocatorFromPoints(polyMesh->GetPoints());
+	vtkNew<vtkKdTree> tetMeshKdTree;
+	tetMeshKdTree->BuildLocatorFromPoints(vtkPointSet::SafeDownCast(tetMesh));
 
 	vtkIdTypeArray *surfacePoints =
 		vtkIdTypeArray::SafeDownCast(tetMesh->GetFieldData()->GetArray("surface_points"));
@@ -176,7 +185,8 @@ vtkTypeBool ProjectSurfacePointsOnPoly::RequestData(
 		                 surfacePoints->GetValue(i),
 		                 motionVectors,
 		                 numberOfAffectors,
-		                 kdTree, links);
+		                 polyMeshKdTree, tetMeshKdTree,
+		                 links);
 	}
 	applyMotionVectors(tetMesh, surfacePoints,
 	                   motionVectors, numberOfAffectors);
diff --git a/src/project_surface_points_on_poly.h b/src/project_surface_points_on_poly.h
index f73a0c4..7c0439c 100644
--- a/src/project_surface_points_on_poly.h
+++ b/src/project_surface_points_on_poly.h
@@ -32,7 +32,8 @@ protected:
                              vtkIdType pointId,
                              std::vector<double> &motionVectors,
                              std::vector<unsigned> &numberOfAffectors,
-                             vtkKdTree *kdTree,
+                             vtkKdTree *polyMeshKdTree,
+                             vtkKdTree *tetMeshKdTree,
                              vtkStaticCellLinks *links);
 };