src/dihedral_angles_filter.cc

@@ -24,18 +24,11 @@ vtkTypeBool DihedralAnglesFilter::RequestData(
 	output->GetPointData()->PassData(input->GetPointData());
 	vtkCellData *cellData = output->GetCellData();
 	cellData->PassData(input->GetCellData());
-
 	vtkNew<vtkDoubleArray> dihedralAnglesArray;
 	dihedralAnglesArray->SetName("dihedral_angles");
 	dihedralAnglesArray->SetNumberOfComponents(6);
 	dihedralAnglesArray->SetNumberOfTuples(input->GetNumberOfCells());
 	double *dihedralAnglesBase = dihedralAnglesArray->GetPointer(0);
-
-	vtkNew<vtkDoubleArray> minDihedralAngleArray;
-	minDihedralAngleArray->SetName("min_dihedral_angle");
-	minDihedralAngleArray->SetNumberOfTuples(input->GetNumberOfCells());
-	double *minDihedralAngleBase = minDihedralAngleArray->GetPointer(0);
-
 	size_t i = 0;
 	vtkCellIterator *it = input->NewCellIterator();
 	for (it->InitTraversal(); !it->IsDoneWithTraversal(); it->GoToNextCell()) {
@@ -59,13 +52,8 @@ vtkTypeBool DihedralAnglesFilter::RequestData(
 		dihedralAnglesBase[i+4] = vtkMath::AngleBetweenVectors(nbdc, nadc);
 		dihedralAnglesBase[i+5] = vtkMath::AngleBetweenVectors(nbdc, nabd);
 
-		minDihedralAngleBase[i / 6] =
-			*std::min_element(dihedralAnglesBase + i,
-			                  dihedralAnglesBase + i + 6);
-
 		i += 6;
 	}
 	cellData->AddArray((vtkAbstractArray *) dihedralAnglesArray);
-	cellData->AddArray((vtkAbstractArray *) minDihedralAngleArray);
 	return true;
 }

src/main.cc

@@ -65,9 +65,12 @@ int main(int argc, char **argv) {
 	auto tetMeshReader = readerFromFileName(argv[1]);
 	auto polyMeshReader = readerFromFileName(argv[2]);
 
+	vtkNew<DihedralAnglesFilter> dihedralAnglesFilter;
+	dihedralAnglesFilter->SetInputConnection(tetMeshReader->GetOutputPort());
+
 	vtkNew<RemoveExternalCellsFilter> removeExternalCellsFilter;
 	removeExternalCellsFilter->SetInputConnection(0,
-		tetMeshReader->GetOutputPort());
+		dihedralAnglesFilter->GetOutputPort());
 	removeExternalCellsFilter->SetInputConnection(1,
 		polyMeshReader->GetOutputPort());
 
@@ -79,11 +82,8 @@ int main(int argc, char **argv) {
 	project->SetInputConnection(0, surfacePointsFilter->GetOutputPort());
 	project->SetInputConnection(1, polyMeshReader->GetOutputPort());
 
-	vtkNew<DihedralAnglesFilter> dihedralAnglesFilter;
-	dihedralAnglesFilter->SetInputConnection(project->GetOutputPort());
-
 	vtkNew<vtkXMLUnstructuredGridWriter> writer;
-	writer->SetInputConnection(dihedralAnglesFilter->GetOutputPort());
+	writer->SetInputConnection(project->GetOutputPort());
 	writer->SetDataModeToAscii();
 	writer->SetFileName("out.vtu");
 	writer->Write();

src/project_surface_points_on_poly.cc

@@ -12,7 +12,6 @@
 #include <vtkPolyData.h>
 #include <vtkKdTree.h>
 #include <vtkStaticCellLinks.h>
-#include <vtkDoubleArray.h>
 
 #include <limits>
 
@@ -37,74 +36,6 @@ int ProjectSurfacePointsOnPoly::FillInputPortInformation(int port, vtkInformatio
 }
 
 
-static void closestPolyMeshPoint(vtkPolyData *polyMesh,
-                                 const double *point,
-                                 vtkKdTree *kdTree,
-                                 vtkStaticCellLinks *links,
-                                 double *toClosestPoint,
-                                 double *distanceToClosestPoint) {
-	vtkIdType closestPolyMeshPoint;
-	{
-		double dummy;
-		closestPolyMeshPoint = kdTree->FindClosestPoint((double *) point, dummy);
-	}
-	vtkIdType *cellIds = links->GetCells(closestPolyMeshPoint);
-	vtkIdType nCells = links->GetNcells(closestPolyMeshPoint);
-	double minDist = std::numeric_limits<double>::infinity();
-
-	/* For each tri the closest point belongs to. */
-	for (vtkIdType i = 0; i < nCells; i++) {
-		vtkIdType cellId = cellIds[i];
-		double direction[3];
-		double dist = pointTriangleDistance(
-			(double *) point, polyMesh->GetCell(cellId), direction);
-		/* Find the closest one. */
-		if (dist < minDist) {
-			minDist = dist;
-			toClosestPoint[0] = direction[0];
-			toClosestPoint[1] = direction[1];
-			toClosestPoint[2] = direction[2];
-		}
-	}
-	vtkMath::MultiplyScalar(toClosestPoint, minDist);
-	*distanceToClosestPoint = minDist;
-}
-
-
-static void moveSurfacePoint(vtkUnstructuredGrid *tetMesh,
-                             vtkPolyData *polyMesh,
-                             vtkIdType pointId,
-                             std::vector<double> &motionVectors,
-                             std::vector<unsigned> &numberOfAffectors,
-                             vtkKdTree *kdTree,
-                             vtkStaticCellLinks *links) {
-	double point[3];
-	tetMesh->GetPoint(pointId, point);
-	double direction[3];
-	double distance;
-	closestPolyMeshPoint(polyMesh, point, kdTree, links, direction, &distance);
-	vtkNew<vtkIdList> affectedPoints;
-	kdTree->FindPointsWithinRadius(distance * 2, point, affectedPoints);
-	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 motion[3] = {direction[0], direction[1], direction[2]};
-		double factor = 1 - std::sqrt(dist2) / (distance * 2);
-		vtkMath::MultiplyScalar(motion, factor);
-		motionVectors[affectedPointId * 3 + 0] += motion[0];
-		motionVectors[affectedPointId * 3 + 1] += motion[1];
-		motionVectors[affectedPointId * 3 + 2] += motion[2];
-		numberOfAffectors[affectedPointId]++;
-	}
-	vtkMath::Subtract(point, direction, point);
-	tetMesh->GetPoints()->SetPoint(pointId, point);
-	tetMesh->Modified();
-	affectedPoints->Reset();
-}
-
-
 vtkTypeBool ProjectSurfacePointsOnPoly::RequestData(
 	vtkInformation *request,
 	vtkInformationVector **inputVector,
@@ -128,27 +59,33 @@ vtkTypeBool ProjectSurfacePointsOnPoly::RequestData(
 
 	vtkIdTypeArray *surfacePoints =
 		vtkIdTypeArray::SafeDownCast(tetMesh->GetFieldData()->GetArray("surface_points"));
-	std::vector<double> motionVectors(3 * tetMesh->GetNumberOfPoints(), 0);
-	std::vector<unsigned> numberOfAffectors(tetMesh->GetNumberOfPoints(), 0);
 	for (vtkIdType i = 0; i < surfacePoints->GetNumberOfTuples(); i++) {
-		moveSurfacePoint(tetMesh, polyMesh,
-		                 surfacePoints->GetValue(i),
-		                 motionVectors,
-		                 numberOfAffectors,
-		                 kdTree, links);
-	}
-	for (vtkIdType i = 0; i < tetMesh->GetNumberOfPoints(); i++) {
-		bool skip = false;
-		for (vtkIdType j = 0; j < surfacePoints->GetNumberOfTuples(); j++) {
-			if (i == surfacePoints->GetValue(j)) skip = true;
+		double p[3];
+		tetMesh->GetPoint(surfacePoints->GetTuple1(i), p);
+		double dist2;
+		vtkIdType closest = kdTree->FindClosestPoint(p, dist2);
+		vtkIdType *cellIds = links->GetCells(closest);
+		vtkIdType nCells = links->GetNcells(closest);
+		double min_dist = std::numeric_limits<double>::infinity();
+		double min_direction[3];
+
+		/* For each tri the closest point belongs to. */
+		for (vtkIdType j = 0; j < nCells; j++) {
+			vtkIdType cellId = cellIds[j];
+			double direction[3];
+			double dist = pointTriangleDistance(p, polyMesh->GetCell(cellId), direction);
+			/* Find the closest one. */
+			if (dist < min_dist) {
+				min_dist = dist;
+				min_direction[0] = direction[0];
+				min_direction[1] = direction[1];
+				min_direction[2] = direction[2];
+			}
 		}
-		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);
+
+		vtkMath::MultiplyScalar(min_direction, min_dist);
+		vtkMath::Subtract(p, min_direction, p);
+		tetMesh->GetPoints()->SetPoint(surfacePoints->GetTuple1(i), p);
 		tetMesh->Modified();
 	}
 

src/project_surface_points_on_poly.h

@@ -14,11 +14,8 @@ public:
 	vtkTypeBool RequestData(vtkInformation *request,
 							vtkInformationVector **inputVector,
 							vtkInformationVector *outputVector) override;
-	vtkSetMacro(affectedNeighborsCount, int);
-	vtkGetMacro(affectedNeighborsCount, int);
 
 protected:
-	int affectedNeighborsCount = 0;
 	~ProjectSurfacePointsOnPoly() override = default;
 };