src/project_surface_points_on_poly.cc

@@ -88,31 +88,24 @@ void ProjectSurfacePointsOnPoly::moveSurfacePoint(vtkUnstructuredGrid *tetMesh,
                              vtkIdType pointId,
                              std::vector<double> &motionVectors,
                              std::vector<unsigned> &numberOfAffectors,
-                             vtkKdTree *polyMeshKdTree,
+                             vtkKdTree *kdTree,
-                             vtkKdTree *tetMeshKdTree,
                              vtkStaticCellLinks *links) {
 	double point[3];
 	tetMesh->GetPoint(pointId, point);
 	double direction[3];
 	double distance;
 
-	closestPolyMeshPoint(
+	closestPolyMeshPoint(polyMesh, point, kdTree, links, direction, &distance);
-		polyMesh, point, polyMeshKdTree, links, direction, &distance);
 	vtkNew<vtkIdList> affectedPoints;
 	findPointsWithinRadius(distance * RadiusScale, tetMesh, point, affectedPoints);
-	// tetMeshKdTree->FindPointsWithinRadius(distance * RadiusScale, point, affectedPoints);
+	// 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);
-		double expected = (distance * RadiusScale) * (distance * RadiusScale);
+		// if(dist2 > (distance * radiusScale) * (distance * radiusScale)) continue;
-		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);
@@ -171,10 +164,8 @@ vtkTypeBool ProjectSurfacePointsOnPoly::RequestData(
 	vtkNew<vtkStaticCellLinks> links;
 	links->BuildLinks(polyMesh);
 
-	vtkNew<vtkKdTree> polyMeshKdTree;
+	vtkNew<vtkKdTree> kdTree;
-	polyMeshKdTree->BuildLocatorFromPoints(polyMesh->GetPoints());
+	kdTree->BuildLocatorFromPoints(polyMesh->GetPoints());
-	vtkNew<vtkKdTree> tetMeshKdTree;
-	tetMeshKdTree->BuildLocatorFromPoints(vtkPointSet::SafeDownCast(tetMesh));
 
 	vtkIdTypeArray *surfacePoints =
 		vtkIdTypeArray::SafeDownCast(tetMesh->GetFieldData()->GetArray("surface_points"));
@@ -185,8 +176,7 @@ vtkTypeBool ProjectSurfacePointsOnPoly::RequestData(
 		                 surfacePoints->GetValue(i),
 		                 motionVectors,
 		                 numberOfAffectors,
-		                 polyMeshKdTree, tetMeshKdTree,
+		                 kdTree, links);
-		                 links);
 	}
 	applyMotionVectors(tetMesh, surfacePoints,
 	                   motionVectors, numberOfAffectors);

src/project_surface_points_on_poly.h

@@ -32,8 +32,7 @@ protected:
                              vtkIdType pointId,
                              std::vector<double> &motionVectors,
                              std::vector<unsigned> &numberOfAffectors,
-                             vtkKdTree *polyMeshKdTree,
+                             vtkKdTree *kdTree,
-                             vtkKdTree *tetMeshKdTree,
                              vtkStaticCellLinks *links);
 };
 

src/relaxation_filter.cc

@@ -25,18 +25,16 @@ 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 < surfacePoints->GetNumberOfValues(); i++) {
+	for (vtkIdType i = 0; i < output->GetNumberOfPoints(); i++) {
-		vtkIdType id = surfacePoints->GetValue(i);
+		if (!isSurface->GetValue(i)) continue;
-
+		
-		output->GetPointCells(id, neighborCells);
+		output->GetPointCells(i, neighborCells);
 
 		if (neighborCells->GetNumberOfIds() != 0) {
 
@@ -54,7 +52,10 @@ 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) {
@@ -63,13 +64,13 @@ vtkTypeBool RelaxationFilter::RequestData(
 					vtkMath::Add(point, avg, avg);
 				}
 				vtkMath::MultiplyScalar(avg, 1. / neighbors.size());
-				newPoints->SetPoint(id, avg);
+				newPoints->SetPoint(i, avg);
 			}
 		}
 		neighborCells->Reset();
 		neighbors.clear();
 	}
-
+	
 	output->SetPoints(newPoints);
 	return true;
 }