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

@@ -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;
+	vtkNew<vtkKdTree> polyMeshKdTree;
-	kdTree->BuildLocatorFromPoints(polyMesh->GetPoints());
+	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);

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

src/relaxation_filter.cc

@@ -25,16 +25,18 @@ vtkTypeBool RelaxationFilter::RequestData(
 	newPoints->DeepCopy(output->GetPoints());
 	vtkNew<vtkIdList> neighborCells;
 	vtkNew<vtkIdList> cellPoints;
+	vtkIdTypeArray *surfacePoints = vtkIdTypeArray::SafeDownCast(
+		output->GetFieldData()->GetArray("surface_points"));
 	vtkIntArray *isSurface = vtkIntArray::SafeDownCast(
 		output->GetPointData()->GetArray("isSurface"));
 	std::set<vtkIdType> neighbors;
 	output->BuildLinks();
 	double avg[3];
 
-	for (vtkIdType i = 0; i < output->GetNumberOfPoints(); i++) {
+	for (vtkIdType i = 0; i < surfacePoints->GetNumberOfValues(); i++) {
-		if (!isSurface->GetValue(i)) continue;
+		vtkIdType id = surfacePoints->GetValue(i);
 
-		output->GetPointCells(i, neighborCells);
+		output->GetPointCells(id, neighborCells);
 
 		if (neighborCells->GetNumberOfIds() != 0) {
 
@@ -52,10 +54,7 @@ vtkTypeBool RelaxationFilter::RequestData(
 				}
 				cellPoints->Reset();
 			}
-
+			neighbors.erase(neighbors.find(id));
-			if (neighbors.find(i) != neighbors.end()) {
-				neighbors.erase(neighbors.find(i));
-			}
 			if (neighbors.size() != 0) {
 				avg[0] = 0; avg[1] = 0; avg[2] = 0;
 				for (const vtkIdType &j : neighbors) {
@@ -64,7 +63,7 @@ vtkTypeBool RelaxationFilter::RequestData(
 					vtkMath::Add(point, avg, avg);
 				}
 				vtkMath::MultiplyScalar(avg, 1. / neighbors.size());
-				newPoints->SetPoint(i, avg);
+				newPoints->SetPoint(id, avg);
 			}
 		}
 		neighborCells->Reset();