diff --git a/src/project_surface_points_on_poly.cc b/src/project_surface_points_on_poly.cc index 8a9e38e..0e9891f 100644 --- a/src/project_surface_points_on_poly.cc +++ b/src/project_surface_points_on_poly.cc @@ -88,24 +88,31 @@ void ProjectSurfacePointsOnPoly::moveSurfacePoint(vtkUnstructuredGrid *tetMesh, vtkIdType pointId, std::vector &motionVectors, std::vector &numberOfAffectors, - vtkKdTree *kdTree, + vtkKdTree *polyMeshKdTree, + vtkKdTree *tetMeshKdTree, vtkStaticCellLinks *links) { double point[3]; tetMesh->GetPoint(pointId, point); double direction[3]; double distance; - closestPolyMeshPoint(polyMesh, point, kdTree, links, direction, &distance); + closestPolyMeshPoint( + polyMesh, point, polyMeshKdTree, links, direction, &distance); vtkNew affectedPoints; findPointsWithinRadius(distance * RadiusScale, tetMesh, point, affectedPoints); - // kdTree->FindPointsWithinRadius(distance * radiusScale, point, affectedPoints); + // tetMeshKdTree->FindPointsWithinRadius(distance * RadiusScale, point, affectedPoints); if(distance > 0) 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); - // if(dist2 > (distance * radiusScale) * (distance * radiusScale)) continue; + double expected = (distance * RadiusScale) * (distance * RadiusScale); + if (dist2 > 2 * expected) { + std::cerr << dist2 << " " << expected << "\n"; + // throw std::runtime_error(""); + continue; + } double motion[3] = {direction[0], direction[1], direction[2]}; double factor = 1. - std::sqrt(dist2) / (distance * RadiusScale); @@ -164,8 +171,10 @@ vtkTypeBool ProjectSurfacePointsOnPoly::RequestData( vtkNew links; links->BuildLinks(polyMesh); - vtkNew kdTree; - kdTree->BuildLocatorFromPoints(polyMesh->GetPoints()); + vtkNew polyMeshKdTree; + polyMeshKdTree->BuildLocatorFromPoints(polyMesh->GetPoints()); + vtkNew tetMeshKdTree; + tetMeshKdTree->BuildLocatorFromPoints(vtkPointSet::SafeDownCast(tetMesh)); vtkIdTypeArray *surfacePoints = vtkIdTypeArray::SafeDownCast(tetMesh->GetFieldData()->GetArray("surface_points")); @@ -176,7 +185,8 @@ vtkTypeBool ProjectSurfacePointsOnPoly::RequestData( surfacePoints->GetValue(i), motionVectors, numberOfAffectors, - kdTree, links); + polyMeshKdTree, tetMeshKdTree, + links); } applyMotionVectors(tetMesh, surfacePoints, motionVectors, numberOfAffectors); diff --git a/src/project_surface_points_on_poly.h b/src/project_surface_points_on_poly.h index f73a0c4..7c0439c 100644 --- a/src/project_surface_points_on_poly.h +++ b/src/project_surface_points_on_poly.h @@ -32,7 +32,8 @@ protected: vtkIdType pointId, std::vector &motionVectors, std::vector &numberOfAffectors, - vtkKdTree *kdTree, + vtkKdTree *polyMeshKdTree, + vtkKdTree *tetMeshKdTree, vtkStaticCellLinks *links); };