ccolin/pfe
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

use vtkDataSet::GetCellNeighbors to find surface points

Browse Source
This commit is contained in:
papush! 2022-03-08 14:59:15 +01:00
parent dcc99754b6
commit d0381b7552
1 changed files with 30 additions and 47 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

77
src/surface_points_filter.cc
View File

@ -12,26 +12,6 @@
/* Remove from ids_a the elements that are not also in ids_b. */
static void keep_only_dupes(std::vector<vtkIdType> &ids_a,
vtkIdType n_ids_b, const vtkIdType *ids_b) {
for (auto it = ids_a.begin(); it != ids_a.end();) {
bool duplicate = false;
for (vtkIdType i = 0; i < n_ids_b; i++) {
if (*it == ids_b[i]) {
duplicate = true;
break;
}
}
if (!duplicate) {
it = ids_a.erase(it);
} else {
++it;
}
}
}
vtkStandardNewMacro(SurfacePointsFilter); vtkStandardNewMacro(SurfacePointsFilter);
@ -50,44 +30,47 @@ vtkTypeBool SurfacePointsFilter::RequestData(
output->GetPointData()->PassData(input->GetPointData()); output->GetPointData()->PassData(input->GetPointData());
output->GetCellData()->PassData(input->GetCellData()); output->GetCellData()->PassData(input->GetCellData());
/* Create an array to store the IDs of external points. */ /* Create an array to store the IDs of surface points. */
vtkNew<vtkIdTypeArray> externalPoints; vtkNew<vtkIdTypeArray> surfacePoints;
externalPoints->SetName("surface_points"); surfacePoints->SetName("surface_points");
externalPoints->SetNumberOfComponents(1); surfacePoints->SetNumberOfComponents(1);
vtkNew<vtkIntArray> isExternal; vtkNew<vtkIntArray> isSurface;
isExternal->SetName("bl"); isSurface->SetName("isSurface");
isExternal->SetNumberOfComponents(1); isSurface->SetNumberOfComponents(1);
isExternal->SetNumberOfTuples(input->GetNumberOfPoints()); isSurface->SetNumberOfTuples(input->GetNumberOfPoints());
isExternal->FillValue(0); isSurface->FillValue(0);
/* Generate cell links, these tell us which cell a point belongs to. */ /* Generate cell links, these tell us which cell a point belongs to. */
vtkNew<vtkStaticCellLinks> links; vtkNew<vtkStaticCellLinks> links;
links->BuildLinks(input); links->BuildLinks(input);
// TODO: try to use vtkDataSet::GetCellNeighbors instead vtkNew<vtkIdList> neighborCells;
vtkNew<vtkIdList> facePoints;
auto *it = input->NewCellIterator(); auto *it = input->NewCellIterator();
for (it->InitTraversal(); !it->IsDoneWithTraversal(); it->GoToNextCell()) { for (it->InitTraversal(); !it->IsDoneWithTraversal(); it->GoToNextCell()) {
vtkIdList *pointIds = it->GetPointIds(); vtkIdList *cellPointIds = it->GetPointIds();
for (vtkIdType faceId = 0; faceId < 4; faceId++) { for (vtkIdType faceId = 0; faceId < 4; faceId++) {
vtkIdType firstPoint = pointIds->GetId(faceId); vtkIdType idA = cellPointIds->GetId((faceId + 0) % 4);
vtkIdType n_cells = links->GetNcells(firstPoint); vtkIdType idB = cellPointIds->GetId((faceId + 1) % 4);
vtkIdType *cells = links->GetCells(firstPoint); vtkIdType idC = cellPointIds->GetId((faceId + 2) % 4);
std::vector<vtkIdType> intersection(cells, cells + n_cells); facePoints->InsertNextId(idA);
for (vtkIdType i = 1; i < 3; i++) { facePoints->InsertNextId(idB);
vtkIdType pointId = pointIds->GetId((faceId + i) % 4); facePoints->InsertNextId(idC);
keep_only_dupes(intersection, input->GetCellNeighbors(it->GetCellId(), facePoints,
links->GetNcells(pointId), links->GetCells(pointId)); neighborCells);
} facePoints->Reset();
if (intersection.size() == 1) { // The face only belongs to one cell. if (neighborCells->GetNumberOfIds() == 1) {
for (vtkIdType i = 0; i < 3; i++) { surfacePoints->InsertNextValue(idA);
externalPoints->InsertNextValue(pointIds->GetId((faceId + i) % 4)); isSurface->SetValue(idA, 1);
isExternal->SetValue(pointIds->GetId((faceId + i) % 4), 1); surfacePoints->InsertNextValue(idB);
} isSurface->SetValue(idB, 1);
surfacePoints->InsertNextValue(idC);
isSurface->SetValue(idC, 1);
} }
} }
} }
output->GetPointData()->SetScalars(isExternal); output->GetPointData()->SetScalars(isSurface);
output->GetFieldData()->AddArray(externalPoints); output->GetFieldData()->AddArray(surfacePoints);
return true; return true;
} }