fix surface point projection kdtree logic

but still use the brute-force approach because it still doesn't work
This commit is contained in:
papush! 2022-03-17 15:18:19 +01:00
parent dacb2ffa04
commit 8ee09d686e
2 changed files with 19 additions and 8 deletions

View File

@ -88,24 +88,31 @@ void ProjectSurfacePointsOnPoly::moveSurfacePoint(vtkUnstructuredGrid *tetMesh,
vtkIdType pointId, vtkIdType pointId,
std::vector<double> &motionVectors, std::vector<double> &motionVectors,
std::vector<unsigned> &numberOfAffectors, std::vector<unsigned> &numberOfAffectors,
vtkKdTree *kdTree, vtkKdTree *polyMeshKdTree,
vtkKdTree *tetMeshKdTree,
vtkStaticCellLinks *links) { vtkStaticCellLinks *links) {
double point[3]; double point[3];
tetMesh->GetPoint(pointId, point); tetMesh->GetPoint(pointId, point);
double direction[3]; double direction[3];
double distance; double distance;
closestPolyMeshPoint(polyMesh, point, kdTree, links, direction, &distance); closestPolyMeshPoint(
polyMesh, point, polyMeshKdTree, links, direction, &distance);
vtkNew<vtkIdList> affectedPoints; vtkNew<vtkIdList> affectedPoints;
findPointsWithinRadius(distance * RadiusScale, tetMesh, point, 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++) { if(distance > 0) for (vtkIdType j = 0; j < affectedPoints->GetNumberOfIds(); j++) {
vtkIdType affectedPointId = affectedPoints->GetId(j); vtkIdType affectedPointId = affectedPoints->GetId(j);
double affectedPoint[3]; double affectedPoint[3];
tetMesh->GetPoint(affectedPointId, affectedPoint); tetMesh->GetPoint(affectedPointId, affectedPoint);
double dist2 = vtkMath::Distance2BetweenPoints(affectedPoint, point); 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 motion[3] = {direction[0], direction[1], direction[2]};
double factor = 1. - std::sqrt(dist2) / (distance * RadiusScale); double factor = 1. - std::sqrt(dist2) / (distance * RadiusScale);
@ -164,8 +171,10 @@ vtkTypeBool ProjectSurfacePointsOnPoly::RequestData(
vtkNew<vtkStaticCellLinks> links; vtkNew<vtkStaticCellLinks> links;
links->BuildLinks(polyMesh); links->BuildLinks(polyMesh);
vtkNew<vtkKdTree> kdTree; vtkNew<vtkKdTree> polyMeshKdTree;
kdTree->BuildLocatorFromPoints(polyMesh->GetPoints()); polyMeshKdTree->BuildLocatorFromPoints(polyMesh->GetPoints());
vtkNew<vtkKdTree> tetMeshKdTree;
tetMeshKdTree->BuildLocatorFromPoints(vtkPointSet::SafeDownCast(tetMesh));
vtkIdTypeArray *surfacePoints = vtkIdTypeArray *surfacePoints =
vtkIdTypeArray::SafeDownCast(tetMesh->GetFieldData()->GetArray("surface_points")); vtkIdTypeArray::SafeDownCast(tetMesh->GetFieldData()->GetArray("surface_points"));
@ -176,7 +185,8 @@ vtkTypeBool ProjectSurfacePointsOnPoly::RequestData(
surfacePoints->GetValue(i), surfacePoints->GetValue(i),
motionVectors, motionVectors,
numberOfAffectors, numberOfAffectors,
kdTree, links); polyMeshKdTree, tetMeshKdTree,
links);
} }
applyMotionVectors(tetMesh, surfacePoints, applyMotionVectors(tetMesh, surfacePoints,
motionVectors, numberOfAffectors); motionVectors, numberOfAffectors);

View File

@ -32,7 +32,8 @@ protected:
vtkIdType pointId, vtkIdType pointId,
std::vector<double> &motionVectors, std::vector<double> &motionVectors,
std::vector<unsigned> &numberOfAffectors, std::vector<unsigned> &numberOfAffectors,
vtkKdTree *kdTree, vtkKdTree *polyMeshKdTree,
vtkKdTree *tetMeshKdTree,
vtkStaticCellLinks *links); vtkStaticCellLinks *links);
}; };