From a60191b0683ff127a7d04ed7759f3c9262487b55 Mon Sep 17 00:00:00 2001 From: CookieKastanie Date: Tue, 15 Mar 2022 15:24:30 +0100 Subject: [PATCH] fix relaxation + radius arg --- src/main.cc | 11 +++++++++-- src/project_surface_points_on_poly.cc | 13 +++++-------- src/project_surface_points_on_poly.h | 16 +++++++++++++++- src/relaxation_filter.cc | 14 +++++++++++++- 4 files changed, 42 insertions(+), 12 deletions(-) diff --git a/src/main.cc b/src/main.cc index f6ef659..066b498 100644 --- a/src/main.cc +++ b/src/main.cc @@ -16,6 +16,7 @@ #include #include #include +#include #ifdef USE_VIEWER #include @@ -58,11 +59,16 @@ vtkSmartPointer 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 project; + project->SetRadiusScale(radiusScale); project->SetInputConnection(0, surfacePointsFilter->GetOutputPort()); project->SetInputConnection(1, polyMeshReader->GetOutputPort()); diff --git a/src/project_surface_points_on_poly.cc b/src/project_surface_points_on_poly.cc index 6502ce3..8a9e38e 100644 --- a/src/project_surface_points_on_poly.cc +++ b/src/project_surface_points_on_poly.cc @@ -9,9 +9,6 @@ #include #include #include -#include -#include -#include #include #include @@ -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 &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 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]; diff --git a/src/project_surface_points_on_poly.h b/src/project_surface_points_on_poly.h index 9433f96..f73a0c4 100644 --- a/src/project_surface_points_on_poly.h +++ b/src/project_surface_points_on_poly.h @@ -3,7 +3,10 @@ #include #include - +#include +#include +#include +#include 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 &motionVectors, + std::vector &numberOfAffectors, + vtkKdTree *kdTree, + vtkStaticCellLinks *links); }; diff --git a/src/relaxation_filter.cc b/src/relaxation_filter.cc index 13c2011..7927073 100644 --- a/src/relaxation_filter.cc +++ b/src/relaxation_filter.cc @@ -30,13 +30,21 @@ vtkTypeBool RelaxationFilter::RequestData( std::set 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; }