pfe/src/surface_points_filter.cc

94 lines
2.8 KiB
C++

#include "surface_points_filter.h"
#include <vtkIdList.h>
#include <vtkUnstructuredGrid.h>
#include <vtkPointData.h>
#include <vtkCellData.h>
#include <vtkIdTypeArray.h>
#include <vtkCellIterator.h>
#include <vtkTriangle.h>
#include <vtkStaticCellLinks.h>
#include <vtkIntArray.h>
/* 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);
vtkTypeBool SurfacePointsFilter::RequestData(
vtkInformation *request,
vtkInformationVector **inputVector,
vtkInformationVector *outputVector) {
(void) request;
vtkUnstructuredGrid* input =
vtkUnstructuredGrid::GetData(inputVector[0]);
vtkUnstructuredGrid* output =
vtkUnstructuredGrid::GetData(outputVector);
/* Copy point and cell data unmodified. */
output->CopyStructure(input);
output->GetPointData()->PassData(input->GetPointData());
output->GetCellData()->PassData(input->GetCellData());
/* Create an array to store the IDs of external points. */
vtkNew<vtkIdTypeArray> externalPoints;
externalPoints->SetName("external_points");
externalPoints->SetNumberOfComponents(1);
vtkNew<vtkIntArray> isExternal;
isExternal->SetName("bl");
isExternal->SetNumberOfComponents(1);
isExternal->SetNumberOfTuples(input->GetNumberOfPoints());
isExternal->FillValue(0);
/* Generate cell links, these tell us which cell a point belongs to. */
vtkNew<vtkStaticCellLinks> links;
links->BuildLinks(input);
// TODO: try to use vtkDataSet::GetCellNeighbors instead
auto *it = input->NewCellIterator();
for (it->InitTraversal(); !it->IsDoneWithTraversal(); it->GoToNextCell()) {
vtkIdList *pointIds = it->GetPointIds();
for (vtkIdType faceId = 0; faceId < 4; faceId++) {
vtkIdType firstPoint = pointIds->GetId(faceId);
vtkIdType n_cells = links->GetNcells(firstPoint);
vtkIdType *cells = links->GetCells(firstPoint);
std::vector<vtkIdType> intersection(cells, cells + n_cells);
for (vtkIdType i = 1; i < 3; i++) {
vtkIdType pointId = pointIds->GetId((faceId + i) % 4);
keep_only_dupes(intersection,
links->GetNcells(pointId), links->GetCells(pointId));
}
if (intersection.size() == 1) { // The face only belongs to one cell.
for (vtkIdType i = 0; i < 3; i++) {
externalPoints->InsertNextValue(pointIds->GetId((faceId + i) % 4));
isExternal->SetValue(pointIds->GetId((faceId + i) % 4), 1);
}
}
}
}
output->GetPointData()->SetScalars(isExternal);
output->GetFieldData()->AddArray(externalPoints);
return true;
}