start implementing proportional motion

src/dihedral_angles_filter.cc

@@ -24,11 +24,18 @@ 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()) {
@@ -52,8 +59,13 @@ 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,12 +65,9 @@ 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,
-		dihedralAnglesFilter->GetOutputPort());
+		tetMeshReader->GetOutputPort());
 	removeExternalCellsFilter->SetInputConnection(1,
 		polyMeshReader->GetOutputPort());
 
@@ -82,8 +79,11 @@ 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(project->GetOutputPort());
+	writer->SetInputConnection(dihedralAnglesFilter->GetOutputPort());
 	writer->SetDataModeToAscii();
 	writer->SetFileName("out.vtu");
 	writer->Write();

src/project_surface_points_on_poly.cc

@@ -12,6 +12,7 @@
 #include <vtkPolyData.h>
 #include <vtkKdTree.h>
 #include <vtkStaticCellLinks.h>
+#include <vtkDoubleArray.h>
 
 #include <limits>
 
@@ -36,6 +37,74 @@ 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,
@@ -59,33 +128,27 @@ 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++) {
-		double p[3];
+		moveSurfacePoint(tetMesh, polyMesh,
-		tetMesh->GetPoint(surfacePoints->GetTuple1(i), p);
+		                 surfacePoints->GetValue(i),
-		double dist2;
+		                 motionVectors,
-		vtkIdType closest = kdTree->FindClosestPoint(p, dist2);
+		                 numberOfAffectors,
-		vtkIdType *cellIds = links->GetCells(closest);
+		                 kdTree, links);
-		vtkIdType nCells = links->GetNcells(closest);
+	}
-		double min_dist = std::numeric_limits<double>::infinity();
+	for (vtkIdType i = 0; i < tetMesh->GetNumberOfPoints(); i++) {
-		double min_direction[3];
+		bool skip = false;
-
+		for (vtkIdType j = 0; j < surfacePoints->GetNumberOfTuples(); j++) {
-		/* For each tri the closest point belongs to. */
+			if (i == surfacePoints->GetValue(j)) skip = true;
-		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;
-		vtkMath::MultiplyScalar(min_direction, min_dist);
+		double point[3];
-		vtkMath::Subtract(p, min_direction, p);
+		tetMesh->GetPoint(i, point);
-		tetMesh->GetPoints()->SetPoint(surfacePoints->GetTuple1(i), p);
+		double *motion = motionVectors.data() + i;
+		vtkMath::MultiplyScalar(motion, 1. / numberOfAffectors[i]);
+		vtkMath::Subtract(point, motion, point);
+		tetMesh->GetPoints()->SetPoint(i, point);
 		tetMesh->Modified();
 	}
 

src/project_surface_points_on_poly.h

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