diff --git a/src/closest_polymesh_point.cc b/src/closest_polymesh_point.cc index 7053093..2f552ce 100644 --- a/src/closest_polymesh_point.cc +++ b/src/closest_polymesh_point.cc @@ -15,6 +15,7 @@ vtkIdType findClosestPoint(const double *point, vtkPointSet *pointSet) { minPoint = i; } } + std::cerr << "closest: " << minPoint << "\n"; return minPoint; } diff --git a/src/project_surface_points_on_poly.cc b/src/project_surface_points_on_poly.cc index 5ae8850..6988c7d 100644 --- a/src/project_surface_points_on_poly.cc +++ b/src/project_surface_points_on_poly.cc @@ -48,23 +48,38 @@ void findPointsWithinRadius(double radius, vtkPointSet *pointSet, } -void ProjectSurfacePointsOnPoly::moveSurfacePoint(vtkUnstructuredGrid *tetMesh, - vtkPolyData *polyMesh, - vtkIdType pointId, - std::vector &motionVectors, - std::vector &numberOfAffectors, - vtkKdTree *polyMeshKdTree, - vtkKdTree *tetMeshKdTree, - vtkStaticCellLinks *links) { +vtkSmartPointer removeSurfacePoints(vtkIdList *list, + vtkIntArray *isSurface) { + vtkNew affectedPoints; + affectedPoints->Allocate(list->GetNumberOfIds()); + for (vtkIdType i = 0; i < list->GetNumberOfIds(); i++) { + vtkIdType id = list->GetId(i); + if (!isSurface->GetValue(id)) { + affectedPoints->InsertNextId(id); + } + } + return affectedPoints; +} + + +void ProjectSurfacePointsOnPoly::moveSurfacePoint( + vtkUnstructuredGrid *tetMesh, + vtkPolyData *polyMesh, + vtkIdType pointId, + std::vector &motionVectors, + std::vector &numberOfAffectors, + vtkKdTree *polyMeshKdTree, + vtkKdTree *tetMeshKdTree, + vtkStaticCellLinks *links) { double point[3]; tetMesh->GetPoint(pointId, point); double direction[3]; double distance; - closestPolyMeshPoint( polyMesh, point, polyMeshKdTree, links, direction, &distance); - vtkNew affectedPoints; - findPointsWithinRadius(distance * RadiusScale, tetMesh, point, affectedPoints); + vtkNew pointsInRadius; + findPointsWithinRadius(distance * RadiusScale, tetMesh, point, pointsInRadius); + auto affectedPoints = removeSurfacePoints(pointsInRadius, isSurface); // tetMeshKdTree->FindPointsWithinRadius(distance * RadiusScale, point, affectedPoints); if(distance > 0) for (vtkIdType j = 0; j < affectedPoints->GetNumberOfIds(); j++) { @@ -96,20 +111,13 @@ void ProjectSurfacePointsOnPoly::moveSurfacePoint(vtkUnstructuredGrid *tetMesh, static void applyMotionVectors(vtkUnstructuredGrid *tetMesh, - vtkIdTypeArray *surfacePoints, std::vector &motionVectors, const std::vector &numberOfAffectors) { 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; - } - if (skip) continue; if(numberOfAffectors[i] == 0) continue; double point[3]; tetMesh->GetPoint(i, point); double *motion = motionVectors.data() + (i * 3); - vtkMath::MultiplyScalar(motion, 1. / numberOfAffectors[i]); vtkMath::Subtract(point, motion, point); tetMesh->GetPoints()->SetPoint(i, point); @@ -142,7 +150,9 @@ vtkTypeBool ProjectSurfacePointsOnPoly::RequestData( tetMeshKdTree->BuildLocatorFromPoints(vtkPointSet::SafeDownCast(tetMesh)); vtkIdTypeArray *surfacePoints = - vtkIdTypeArray::SafeDownCast(tetMesh->GetFieldData()->GetArray("surface_points")); + vtkIdTypeArray::SafeDownCast(tetMesh->GetFieldData()->GetArray("surfacePoints")); + isSurface = vtkIntArray::SafeDownCast( + output->GetPointData()->GetArray("isSurface")); std::vector motionVectors(3 * tetMesh->GetNumberOfPoints(), 0); std::vector numberOfAffectors(tetMesh->GetNumberOfPoints(), 0); for (vtkIdType i = 0; i < surfacePoints->GetNumberOfTuples(); i++) { @@ -153,8 +163,7 @@ vtkTypeBool ProjectSurfacePointsOnPoly::RequestData( polyMeshKdTree, tetMeshKdTree, links); } - applyMotionVectors(tetMesh, surfacePoints, - motionVectors, numberOfAffectors); + applyMotionVectors(tetMesh, motionVectors, numberOfAffectors); return true; } diff --git a/src/project_surface_points_on_poly.h b/src/project_surface_points_on_poly.h index 7c0439c..5e0eb31 100644 --- a/src/project_surface_points_on_poly.h +++ b/src/project_surface_points_on_poly.h @@ -25,6 +25,8 @@ public: protected: int affectedNeighborsCount = 0; double RadiusScale; + vtkIntArray *isSurface = nullptr; + ~ProjectSurfacePointsOnPoly() override = default; void moveSurfacePoint(vtkUnstructuredGrid *tetMesh, diff --git a/src/relaxation_filter.cc b/src/relaxation_filter.cc index 9e48f26..c658d42 100644 --- a/src/relaxation_filter.cc +++ b/src/relaxation_filter.cc @@ -30,7 +30,7 @@ vtkTypeBool RelaxationFilter::RequestData( vtkNew neighborCells; vtkNew cellPoints; vtkIdTypeArray *surfacePoints = vtkIdTypeArray::SafeDownCast( - output->GetFieldData()->GetArray("surface_points")); + output->GetFieldData()->GetArray("surfacePoints")); vtkIntArray *isSurface = vtkIntArray::SafeDownCast( output->GetPointData()->GetArray("isSurface")); std::set neighbors; diff --git a/src/surface_points_filter.cc b/src/surface_points_filter.cc index 2de5261..7d1358b 100644 --- a/src/surface_points_filter.cc +++ b/src/surface_points_filter.cc @@ -32,7 +32,7 @@ vtkTypeBool SurfacePointsFilter::RequestData( /* Create an array to store the IDs of surface points. */ vtkNew surfacePoints; - surfacePoints->SetName("surface_points"); + surfacePoints->SetName("surfacePoints"); surfacePoints->SetNumberOfComponents(1); vtkNew isSurface; isSurface->SetName("isSurface");