fix relaxation + radius arg

src/main.cc

@@ -16,6 +16,7 @@
 #include <vtkXMLUnstructuredGridReader.h>
 #include <vtkXMLUnstructuredGridWriter.h>
 #include <vtksys/SystemTools.hxx>
+#include <string>
 
 #ifdef USE_VIEWER
 #include <vtkCamera.h>
@@ -58,11 +59,16 @@ vtkSmartPointer<vtkAlgorithm> readerFromFileName(const char *fileName) {
 
 
 int main(int argc, char **argv) {
-	if (argc != 3) {
-		std::cerr << "Usage: " << argv[0] << " tetmesh polydata" << std::endl;
+	if (argc < 3 || argc > 5) {
+		std::cerr << "Usage: " << argv[0] << " tetmesh polydata [radiusScale]" << std::endl;
 		return EXIT_FAILURE;
 	}
 
+	double radiusScale = 2.;
+	if(argc == 4) {
+		radiusScale = std::stod(argv[3]);
+	}
+
 	auto tetMeshReader = readerFromFileName(argv[1]);
 	auto polyMeshReader = readerFromFileName(argv[2]);
 
@@ -77,6 +83,7 @@ int main(int argc, char **argv) {
 		removeExternalCellsFilter->GetOutputPort());
 
 	vtkNew<ProjectSurfacePointsOnPoly> project;
+	project->SetRadiusScale(radiusScale);
 	project->SetInputConnection(0, surfacePointsFilter->GetOutputPort());
 	project->SetInputConnection(1, polyMeshReader->GetOutputPort());
 

src/project_surface_points_on_poly.cc

@@ -9,9 +9,6 @@
 #include <vtkCellIterator.h>
 #include <vtkInformation.h>
 #include <vtkInformationVector.h>
-#include <vtkPolyData.h>
-#include <vtkKdTree.h>
-#include <vtkStaticCellLinks.h>
 #include <vtkDoubleArray.h>
 
 #include <limits>
@@ -22,6 +19,8 @@ vtkStandardNewMacro(ProjectSurfacePointsOnPoly);
 
 ProjectSurfacePointsOnPoly::ProjectSurfacePointsOnPoly() {
 	SetNumberOfInputPorts(2);
+
+	RadiusScale = 2.;
 }
 
 
@@ -84,7 +83,7 @@ void findPointsWithinRadius(double radius, vtkPointSet *pointSet,
 }
 
 
-static void moveSurfacePoint(vtkUnstructuredGrid *tetMesh,
+void ProjectSurfacePointsOnPoly::moveSurfacePoint(vtkUnstructuredGrid *tetMesh,
                              vtkPolyData *polyMesh,
                              vtkIdType pointId,
                              std::vector<double> &motionVectors,
@@ -96,11 +95,9 @@ static void moveSurfacePoint(vtkUnstructuredGrid *tetMesh,
 	double direction[3];
 	double distance;
 
-	double radiusScale = 5.;
-
 	closestPolyMeshPoint(polyMesh, point, kdTree, links, direction, &distance);
 	vtkNew<vtkIdList> affectedPoints;
-	findPointsWithinRadius(distance * radiusScale, tetMesh, point, affectedPoints);
+	findPointsWithinRadius(distance * RadiusScale, tetMesh, point, affectedPoints);
 	// kdTree->FindPointsWithinRadius(distance * radiusScale, point, affectedPoints);
 
 	if(distance > 0) for (vtkIdType j = 0; j < affectedPoints->GetNumberOfIds(); j++) {
@@ -111,7 +108,7 @@ static void moveSurfacePoint(vtkUnstructuredGrid *tetMesh,
 		// if(dist2 > (distance * radiusScale) * (distance * radiusScale)) continue;
 
 		double motion[3] = {direction[0], direction[1], direction[2]};
-		double factor = 1. - std::sqrt(dist2) / (distance * radiusScale);
+		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];

src/project_surface_points_on_poly.h

@@ -3,7 +3,10 @@
 
 #include <vtkUnstructuredGridAlgorithm.h>
 #include <vtkIdList.h>
-
+#include <vector>
+#include <vtkPolyData.h>
+#include <vtkKdTree.h>
+#include <vtkStaticCellLinks.h>
 
 class ProjectSurfacePointsOnPoly : public vtkUnstructuredGridAlgorithm {
 public:
@@ -17,9 +20,20 @@ public:
 	vtkSetMacro(affectedNeighborsCount, int);
 	vtkGetMacro(affectedNeighborsCount, int);
 
+	vtkSetMacro(RadiusScale, double);
+
 protected:
 	int affectedNeighborsCount = 0;
+	double RadiusScale;
 	~ProjectSurfacePointsOnPoly() override = default;
+
+	void moveSurfacePoint(vtkUnstructuredGrid *tetMesh,
+                             vtkPolyData *polyMesh,
+                             vtkIdType pointId,
+                             std::vector<double> &motionVectors,
+                             std::vector<unsigned> &numberOfAffectors,
+                             vtkKdTree *kdTree,
+                             vtkStaticCellLinks *links);
 };
 
 

src/relaxation_filter.cc

@@ -30,13 +30,21 @@ vtkTypeBool RelaxationFilter::RequestData(
 	std::set<vtkIdType> neighbors;
 	output->BuildLinks();
 	double avg[3];
+	
 	for (vtkIdType i = 0; i < output->GetNumberOfPoints(); i++) {
+		if (!isSurface->GetValue(i)) continue;
+		
 		output->GetPointCells(i, neighborCells);
+
 		if (neighborCells->GetNumberOfIds() != 0) {
+
 			for (vtkIdType j = 0; j < neighborCells->GetNumberOfIds(); j++) {
+
 				vtkIdType cellId = neighborCells->GetId(j);
 				output->GetCellPoints(cellId, cellPoints);
+
 				for (vtkIdType k = 0; k < cellPoints->GetNumberOfIds(); k++) {
+
 					vtkIdType neighborId = cellPoints->GetId(k);
 					if (isSurface->GetValue(neighborId)) {
 						neighbors.insert(neighborId);
@@ -44,7 +52,10 @@ vtkTypeBool RelaxationFilter::RequestData(
 				}
 				cellPoints->Reset();
 			}
-			neighbors.erase(neighbors.find(i));
+
+			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) {
@@ -59,6 +70,7 @@ vtkTypeBool RelaxationFilter::RequestData(
 		neighborCells->Reset();
 		neighbors.clear();
 	}
+	
 	output->SetPoints(newPoints);
 	return true;
 }