*clap*

of relaxation of the tetrahedral mesh of the program of tetrahedral
mesh generation

CMakeLists.txt

@@ -12,6 +12,7 @@ set(VTK_COMPONENTS
   VTK::IOGeometry
   VTK::IOXML
   VTK::FiltersModeling
+  VTK::FiltersGeometry
   VTK::vtksys)
 set(ENABLE_VIEWER OFF CACHE BOOL "Enable the 3D viewer, depends on Qt.")
 if(ENABLE_VIEWER)
@@ -62,7 +63,11 @@ target_sources(pfe PRIVATE
   src/remove_external_cells_filter.cc
   src/remove_external_cells_filter.h
   src/relaxation_filter.cc
-  src/relaxation_filter.h)
+  src/relaxation_filter.h
+  src/max_distance_filter.cc
+  src/max_distance_filter.h
+  src/closest_polymesh_point.cc
+  src/closest_polymesh_point.h)
 
 target_link_libraries(pfe PRIVATE ${VTK_COMPONENTS})
 

data/tetrahedron.obj

@@ -0,0 +1,12 @@
+# Blender v2.79 (sub 0) OBJ File: ''
+# www.blender.org
+o Solid
+v 0.000000 0.769688 0.000000
+v 0.000000 0.000000 -0.769688
+v 0.769688 0.000000 0.000000
+v 0.000000 0.000000 0.000000
+s off
+f 1 2 3
+f 1 3 4
+f 1 4 2
+f 2 4 3

src/closest_polymesh_point.cc

@@ -0,0 +1,60 @@
+#include "closest_polymesh_point.h"
+#include "point_tris_dist.h"
+#include <limits>
+
+
+vtkIdType findClosestPoint(const double *point, vtkPointSet *pointSet) {
+	vtkIdType minPoint = 0;
+	double minDist2 = std::numeric_limits<double>::infinity();
+	for (vtkIdType i = 0; i < pointSet->GetNumberOfPoints(); i++) {
+		double point2[3];
+		pointSet->GetPoint(i, point2);
+		double dist2 = vtkMath::Distance2BetweenPoints(point, point2);
+		if (dist2 < minDist2) {
+			minDist2 = dist2;
+			minPoint = i;
+		}
+	}
+	return minPoint;
+}
+
+
+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);
+		closestPolyMeshPoint = findClosestPoint(point, polyMesh);
+	}
+	vtkIdType *cellIds = links->GetCells(closestPolyMeshPoint);
+	vtkIdType nCells = links->GetNcells(closestPolyMeshPoint);
+	if (nCells == 0) {
+		*distanceToClosestPoint = -1;
+		return;
+	}
+	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;
+}
+
+

src/closest_polymesh_point.h

@@ -0,0 +1,17 @@
+#ifndef CLOSEST_POLYMESH_POINT_H
+#define CLOSEST_POLYMESH_POINT_H
+
+#include <vtkPolyData.h>
+#include <vtkKdTree.h>
+#include <vtkStaticCellLinks.h>
+
+
+void closestPolyMeshPoint(vtkPolyData *polyMesh,
+                          const double *point,
+                          vtkKdTree *kdTree,
+                          vtkStaticCellLinks *links,
+                          double *toClosestPoint,
+                          double *distanceToClosestPoint);
+
+
+#endif

src/dihedral_angles_filter.cc

@@ -1,5 +1,7 @@
 #include "dihedral_angles_filter.h"
 
+#include <algorithm>
+#include <limits>
 #include <vtkUnstructuredGrid.h>
 #include <vtkPointData.h>
 #include <vtkCellData.h>
@@ -36,6 +38,9 @@ vtkTypeBool DihedralAnglesFilter::RequestData(
 	minDihedralAngleArray->SetNumberOfTuples(input->GetNumberOfCells());
 	double *minDihedralAngleBase = minDihedralAngleArray->GetPointer(0);
 
+	AverageMinDegrees = 0;
+	MinMinDegrees = std::numeric_limits<double>::infinity();
+
 	size_t i = 0;
 	vtkCellIterator *it = input->NewCellIterator();
 	for (it->InitTraversal(); !it->IsDoneWithTraversal(); it->GoToNextCell()) {
@@ -62,9 +67,14 @@ vtkTypeBool DihedralAnglesFilter::RequestData(
 		minDihedralAngleBase[i / 6] =
 			*std::min_element(dihedralAnglesBase + i,
 			                  dihedralAnglesBase + i + 6);
+		AverageMinDegrees += minDihedralAngleBase[i / 6];
+		MinMinDegrees = std::min(MinMinDegrees, minDihedralAngleBase[i / 6]);
 
 		i += 6;
 	}
+	AverageMinDegrees /= input->GetNumberOfCells();
+	AverageMinDegrees = AverageMinDegrees * 180. / 3.141592653589793;
+	MinMinDegrees = MinMinDegrees * 180. / 3.141592653589793;
 	cellData->AddArray((vtkAbstractArray *) dihedralAnglesArray);
 	cellData->AddArray((vtkAbstractArray *) minDihedralAngleArray);
 	return true;

src/dihedral_angles_filter.h

@@ -12,9 +12,13 @@ public:
 	vtkTypeBool RequestData(vtkInformation *request,
 							vtkInformationVector **inputVector,
 							vtkInformationVector *outputVector) override;
+	vtkGetMacro(AverageMinDegrees, double);
+	vtkGetMacro(MinMinDegrees, double);
 
 protected:
 	~DihedralAnglesFilter() override = default;
+	double AverageMinDegrees;
+	double MinMinDegrees;
 };
 
 

src/main.cc

@@ -3,13 +3,16 @@
 #include "surface_points_filter.h"
 #include "project_surface_points_on_poly.h"
 #include "relaxation_filter.h"
+#include "max_distance_filter.h"
 
+#include <vtkGeometryFilter.h>
 #include <vtkCellArrayIterator.h>
 #include <vtkCellData.h>
 #include <vtkDataSetReader.h>
 #include <vtkOBJReader.h>
 #include <vtkPolyData.h>
 #include <vtkPolyDataReader.h>
+#include <vtkSetGet.h>
 #include <vtkUnstructuredGrid.h>
 #include <vtkUnstructuredGridReader.h>
 #include <vtkXMLPolyDataReader.h>
@@ -59,14 +62,16 @@ vtkSmartPointer<vtkAlgorithm> readerFromFileName(const char *fileName) {
 
 
 int main(int argc, char **argv) {
-	if (argc < 3 || argc > 5) {
+	if (!(argc == 3 || argc == 5)) {
-		std::cerr << "Usage: " << argv[0] << " tetmesh polydata [radiusScale]" << std::endl;
+		std::cerr << "Usage: " << argv[0] << " tetmesh polydata [radiusScale relaxationIterCount]" << std::endl;
 		return EXIT_FAILURE;
 	}
 
 	double radiusScale = 2.;
-	if(argc == 4) {
+	int iterCount = 1;
+	if(argc > 3) {
 		radiusScale = std::stod(argv[3]);
+		iterCount = std::stoi(argv[4]);
 	}
 
 	auto tetMeshReader = readerFromFileName(argv[1]);
@@ -88,18 +93,31 @@ int main(int argc, char **argv) {
 	project->SetInputConnection(1, polyMeshReader->GetOutputPort());
 
 	vtkNew<RelaxationFilter> relaxationFilter;
+	relaxationFilter->SetIterCount(iterCount);
 	relaxationFilter->SetInputConnection(project->GetOutputPort());
 
 	vtkNew<DihedralAnglesFilter> dihedralAnglesFilter;
 	dihedralAnglesFilter->SetInputConnection(
 		relaxationFilter->GetOutputPort());
 
+	vtkNew<vtkGeometryFilter> geometryFilter;
+	geometryFilter->SetInputConnection(dihedralAnglesFilter->GetOutputPort());
+
+	vtkNew<MaxDistanceFilter> maxDistanceFilter;
+	maxDistanceFilter->SetInputConnection(0, geometryFilter->GetOutputPort());
+	maxDistanceFilter->SetInputConnection(1, polyMeshReader->GetOutputPort());
+
 	vtkNew<vtkXMLUnstructuredGridWriter> writer;
 	writer->SetInputConnection(dihedralAnglesFilter->GetOutputPort());
 	writer->SetDataModeToAscii();
 	writer->SetFileName("out.vtu");
 	writer->Write();
 
+	maxDistanceFilter->Update();
+	std::cerr << "Max distance: " << maxDistanceFilter->GetMaxDist() << "\n"
+	          << "Average min angle: " << dihedralAnglesFilter->GetAverageMinDegrees() << "\n"
+	          << "Min min angle: " << dihedralAnglesFilter->GetMinMinDegrees() << "\n";
+
 #ifdef USE_VIEWER
 	/* Volume rendering properties */
 	vtkNew<vtkOpenGLProjectedTetrahedraMapper> volumeMapper;

src/max_distance_filter.cc

@@ -0,0 +1,63 @@
+#include "max_distance_filter.h"
+#include "closest_polymesh_point.h"
+
+#include <limits>
+#include <vtkPolyData.h>
+#include <vtkInformation.h>
+#include <vtkKdTree.h>
+
+
+vtkStandardNewMacro(MaxDistanceFilter);
+
+
+MaxDistanceFilter::MaxDistanceFilter() {
+	SetNumberOfInputPorts(2);
+	SetNumberOfOutputPorts(0);
+}
+
+
+int MaxDistanceFilter::FillInputPortInformation(int port,
+                                                vtkInformation *info) {
+	if (port == 0) {
+		vtkPolyDataAlgorithm::FillInputPortInformation(port, info);
+	} else if (port == 1) {
+		info->Set(vtkDataObject::DATA_TYPE_NAME(), "vtkPolyData");
+	}
+	return 1;
+}
+
+
+vtkTypeBool MaxDistanceFilter::RequestData(vtkInformation *request,
+                                           vtkInformationVector **inputVector,
+                                           vtkInformationVector *outputVector) {
+	(void) request;
+	(void) outputVector;
+	vtkPolyData* input1 = vtkPolyData::GetData(inputVector[0]);
+	vtkPolyData* input2 = vtkPolyData::GetData(inputVector[1]);
+	vtkNew<vtkKdTree> tree1;
+	tree1->BuildLocatorFromPoints(input1->GetPoints());
+	vtkNew<vtkKdTree> tree2;
+	tree2->BuildLocatorFromPoints(input2->GetPoints());
+	vtkNew<vtkStaticCellLinks> links1;
+	links1->BuildLinks(input1);
+	vtkNew<vtkStaticCellLinks> links2;
+	links2->BuildLinks(input2);
+	MaxDist = 0;
+	for (vtkIdType i = 0; i < input1->GetNumberOfPoints(); i++) {
+		double point[3];
+		input1->GetPoint(i, point);
+		double vec[3];
+		double dist;
+		closestPolyMeshPoint(input2, point, tree2, links2, vec, &dist);
+		MaxDist = std::max(dist, MaxDist);
+	}
+	for (vtkIdType i = 0; i < input2->GetNumberOfPoints(); i++) {
+		double point[3];
+		input2->GetPoint(i, point);
+		double vec[3];
+		double dist;
+		closestPolyMeshPoint(input1, point, tree1, links1, vec, &dist);
+		MaxDist = std::max(dist, MaxDist);
+	}
+	return true;
+}

src/max_distance_filter.h

@@ -0,0 +1,23 @@
+#ifndef MAX_DISTANCE_FILTER_H
+#define MAX_DISTANCE_FILTER_H
+
+#include <vtkPolyDataAlgorithm.h>
+
+
+class MaxDistanceFilter : public vtkPolyDataAlgorithm {
+public:
+	static MaxDistanceFilter *New();
+	vtkTypeMacro(MaxDistanceFilter, vtkPolyDataAlgorithm);
+	MaxDistanceFilter();
+	int FillInputPortInformation(int, vtkInformation *info) override;
+	vtkTypeBool RequestData(vtkInformation *request,
+							vtkInformationVector **inputVector,
+							vtkInformationVector *outputVector) override;
+	vtkGetMacro(MaxDist, double);
+
+protected:
+	double MaxDist;
+};
+
+
+#endif

src/project_surface_points_on_poly.cc

@@ -1,6 +1,5 @@
-#include "point_tris_dist.h"
-
 #include "project_surface_points_on_poly.h"
+#include "closest_polymesh_point.h"
 
 #include <vtkUnstructuredGrid.h>
 #include <vtkPointData.h>
@@ -36,40 +35,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;
-}
-
-
 void findPointsWithinRadius(double radius, vtkPointSet *pointSet,
                             const double *point, vtkIdList *points) {
 	for (vtkIdType i = 0; i < pointSet->GetNumberOfPoints(); i++) {

src/relaxation_filter.cc

@@ -10,6 +10,10 @@
 
 vtkStandardNewMacro(RelaxationFilter);
 
+RelaxationFilter::RelaxationFilter()
+	:IterCount(1) {
+}
+
 vtkTypeBool RelaxationFilter::RequestData(
 	vtkInformation *request,
 	vtkInformationVector **inputVector,
@@ -33,43 +37,45 @@ vtkTypeBool RelaxationFilter::RequestData(
 	output->BuildLinks();
 	double avg[3];
 
-	for (vtkIdType i = 0; i < surfacePoints->GetNumberOfValues(); i++) {
+	for (int iter = 0; iter < IterCount; iter++) {
-		vtkIdType id = surfacePoints->GetValue(i);
+		for (vtkIdType i = 0; i < surfacePoints->GetNumberOfValues(); i++) {
+			vtkIdType id = surfacePoints->GetValue(i);
 
-		output->GetPointCells(id, neighborCells);
+			output->GetPointCells(id, neighborCells);
 
-		if (neighborCells->GetNumberOfIds() != 0) {
+			if (neighborCells->GetNumberOfIds() != 0) {
 
-			for (vtkIdType j = 0; j < neighborCells->GetNumberOfIds(); j++) {
+				for (vtkIdType j = 0; j < neighborCells->GetNumberOfIds(); j++) {
 
-				vtkIdType cellId = neighborCells->GetId(j);
+					vtkIdType cellId = neighborCells->GetId(j);
-				output->GetCellPoints(cellId, cellPoints);
+					output->GetCellPoints(cellId, cellPoints);
 
-				for (vtkIdType k = 0; k < cellPoints->GetNumberOfIds(); k++) {
+					for (vtkIdType k = 0; k < cellPoints->GetNumberOfIds(); k++) {
 
-					vtkIdType neighborId = cellPoints->GetId(k);
+						vtkIdType neighborId = cellPoints->GetId(k);
-					if (isSurface->GetValue(neighborId)) {
+						if (isSurface->GetValue(neighborId)) {
-						neighbors.insert(neighborId);
+							neighbors.insert(neighborId);
+						}
 					}
+					cellPoints->Reset();
 				}
-				cellPoints->Reset();
+				neighbors.erase(neighbors.find(id));
-			}
+				if (neighbors.size() != 0) {
-			neighbors.erase(neighbors.find(id));
+					avg[0] = 0; avg[1] = 0; avg[2] = 0;
-			if (neighbors.size() != 0) {
+					for (const vtkIdType &j : neighbors) {
-				avg[0] = 0; avg[1] = 0; avg[2] = 0;
+						double point[3];
-				for (const vtkIdType &j : neighbors) {
+						output->GetPoints()->GetPoint(j, point);
-					double point[3];
+						vtkMath::Add(point, avg, avg);
-					output->GetPoints()->GetPoint(j, point);
+					}
-					vtkMath::Add(point, avg, avg);
+					vtkMath::MultiplyScalar(avg, 1. / neighbors.size());
+					newPoints->SetPoint(id, avg);
 				}
-				vtkMath::MultiplyScalar(avg, 1. / neighbors.size());
-				newPoints->SetPoint(id, avg);
 			}
+			neighborCells->Reset();
+			neighbors.clear();
 		}
-		neighborCells->Reset();
-		neighbors.clear();
-	}
 
-	output->SetPoints(newPoints);
+		output->SetPoints(newPoints);
+	}
 	return true;
 }

src/relaxation_filter.h

@@ -9,12 +9,16 @@ class RelaxationFilter : public vtkUnstructuredGridAlgorithm {
 public:
 	static RelaxationFilter *New();
 	vtkTypeMacro(RelaxationFilter, vtkUnstructuredGridAlgorithm);
+	RelaxationFilter();
 	vtkTypeBool RequestData(vtkInformation *request,
 							vtkInformationVector **inputVector,
 							vtkInformationVector *outputVector) override;
+	vtkSetMacro(IterCount, int);
+	vtkGetMacro(IterCount, int);
 
 protected:
 	~RelaxationFilter() override = default;
+	int IterCount;
 };