From 2ac9f1ac585287f18e4cfe61c63ce981ed8f59b9 Mon Sep 17 00:00:00 2001 From: Kenneth Moreland Date: Fri, 28 Mar 2025 12:22:10 -0400 Subject: [PATCH] Enable extracting external faces from curvilinear data The external faces filter was not working with curvilinear data. It was throwing an exception about not being able to cast the coordinates array. This has been fixed to work with this condition. --- docs/changelog/curvilinear-external-faces.md | 7 + .../entity_extraction/ExternalFaces.cxx | 12 +- vtkm/filter/entity_extraction/ExternalFaces.h | 2 +- .../testing/UnitTestExternalFacesFilter.cxx | 85 +++- .../entity_extraction/worklet/ExternalFaces.h | 470 ++++++------------ 5 files changed, 239 insertions(+), 337 deletions(-) create mode 100644 docs/changelog/curvilinear-external-faces.md diff --git a/docs/changelog/curvilinear-external-faces.md b/docs/changelog/curvilinear-external-faces.md new file mode 100644 index 0000000000..7824d362f0 --- /dev/null +++ b/docs/changelog/curvilinear-external-faces.md @@ -0,0 +1,7 @@ +## Enable extracting external faces from curvilinear data + +The external faces filter was not working with curvilinear data. The +implementation for structured cells was relying on axis-aligned point +coordinates, which is not the case for curvilinear grids. The implementation now +only relies on the indices in the 3D grid, so it works on structured data +regardless of the point coordinates. This should also speed up the operation. diff --git a/vtkm/filter/entity_extraction/ExternalFaces.cxx b/vtkm/filter/entity_extraction/ExternalFaces.cxx index 9c2d4a5279..2541d2533d 100644 --- a/vtkm/filter/entity_extraction/ExternalFaces.cxx +++ b/vtkm/filter/entity_extraction/ExternalFaces.cxx @@ -38,7 +38,7 @@ void ExternalFaces::SetPassPolyData(bool value) //----------------------------------------------------------------------------- vtkm::cont::DataSet ExternalFaces::GenerateOutput(const vtkm::cont::DataSet& input, - vtkm::cont::CellSetExplicit<>& outCellSet) + vtkm::cont::UnknownCellSet& outCellSet) { //3. Check the fields of the dataset to see what kinds of fields are present, so // we can free the cell mapping array if it won't be needed. @@ -71,18 +71,16 @@ vtkm::cont::DataSet ExternalFaces::DoExecute(const vtkm::cont::DataSet& input) //2. using the policy convert the dynamic cell set, and run the // external faces worklet - vtkm::cont::CellSetExplicit<> outCellSet; + vtkm::cont::UnknownCellSet outCellSet; if (cells.CanConvert>()) { - this->Worklet->Run(cells.AsCellSet>(), - input.GetCoordinateSystem(this->GetActiveCoordinateSystemIndex()), - outCellSet); + outCellSet = this->Worklet->Run(cells.AsCellSet>()); } else { - this->Worklet->Run(cells.ResetCellSetList(), - outCellSet); + outCellSet = + this->Worklet->Run(cells.ResetCellSetList()); } // New Filter Design: we generate new output and map the fields first. diff --git a/vtkm/filter/entity_extraction/ExternalFaces.h b/vtkm/filter/entity_extraction/ExternalFaces.h index 02290d6c4e..5431525d98 100644 --- a/vtkm/filter/entity_extraction/ExternalFaces.h +++ b/vtkm/filter/entity_extraction/ExternalFaces.h @@ -64,7 +64,7 @@ private: VTKM_CONT vtkm::cont::DataSet DoExecute(const vtkm::cont::DataSet& input) override; vtkm::cont::DataSet GenerateOutput(const vtkm::cont::DataSet& input, - vtkm::cont::CellSetExplicit<>& outCellSet); + vtkm::cont::UnknownCellSet& outCellSet); VTKM_CONT bool MapFieldOntoOutput(vtkm::cont::DataSet& result, const vtkm::cont::Field& field); diff --git a/vtkm/filter/entity_extraction/testing/UnitTestExternalFacesFilter.cxx b/vtkm/filter/entity_extraction/testing/UnitTestExternalFacesFilter.cxx index fd29eb78cb..87ae350a4c 100644 --- a/vtkm/filter/entity_extraction/testing/UnitTestExternalFacesFilter.cxx +++ b/vtkm/filter/entity_extraction/testing/UnitTestExternalFacesFilter.cxx @@ -8,6 +8,7 @@ // PURPOSE. See the above copyright notice for more information. //============================================================================ +#include #include #include @@ -50,6 +51,21 @@ vtkm::cont::DataSet MakeDataTestSet5() return MakeTestDataSet().Make3DExplicitDataSet6(); } +vtkm::cont::DataSet MakeUniformDataTestSet() +{ + return MakeTestDataSet().Make3DUniformDataSet1(); +} + +vtkm::cont::DataSet MakeCurvilinearDataTestSet() +{ + vtkm::cont::DataSet data = MakeUniformDataTestSet(); + vtkm::cont::ArrayHandle coords; + vtkm::cont::CoordinateSystem oldCoords = data.GetCoordinateSystem(); + vtkm::cont::ArrayCopy(oldCoords.GetData(), coords); + data.AddCoordinateSystem(oldCoords.GetName(), coords); + return data; +} + void TestExternalFacesExplicitGrid(const vtkm::cont::DataSet& ds, bool compactPoints, vtkm::Id numExpectedExtFaces, @@ -63,9 +79,7 @@ void TestExternalFacesExplicitGrid(const vtkm::cont::DataSet& ds, vtkm::cont::DataSet resultds = externalFaces.Execute(ds); // verify cellset - vtkm::cont::CellSetExplicit<> new_cellSet = - resultds.GetCellSet().AsCellSet>(); - const vtkm::Id numOutputExtFaces = new_cellSet.GetNumberOfCells(); + const vtkm::Id numOutputExtFaces = resultds.GetNumberOfCells(); VTKM_TEST_ASSERT(numOutputExtFaces == numExpectedExtFaces, "Number of External Faces mismatch"); // verify fields @@ -135,6 +149,69 @@ void TestWithMixed2Dand3DMesh() TestExternalFacesExplicitGrid(ds, true, 6, 5, false); } +void TestExternalFacesStructuredGrid(const vtkm::cont::DataSet& ds, bool compactPoints) +{ + // Get the dimensions of the grid. + vtkm::cont::CellSetStructured<3> cellSet; + ds.GetCellSet().AsCellSet(cellSet); + vtkm::Id3 pointDims = cellSet.GetPointDimensions(); + vtkm::Id3 cellDims = cellSet.GetCellDimensions(); + + //Run the External Faces filter + vtkm::filter::entity_extraction::ExternalFaces externalFaces; + externalFaces.SetCompactPoints(compactPoints); + vtkm::cont::DataSet resultds = externalFaces.Execute(ds); + + // verify cellset + vtkm::Id numExpectedExtFaces = ((2 * cellDims[0] * cellDims[1]) + // x-y faces + (2 * cellDims[0] * cellDims[2]) + // x-z faces + (2 * cellDims[1] * cellDims[2])); // y-z faces + const vtkm::Id numOutputExtFaces = resultds.GetNumberOfCells(); + VTKM_TEST_ASSERT(numOutputExtFaces == numExpectedExtFaces, "Number of External Faces mismatch"); + + // verify fields + VTKM_TEST_ASSERT(resultds.HasField("pointvar"), "Point field not mapped successfully"); + VTKM_TEST_ASSERT(resultds.HasField("cellvar"), "Cell field not mapped successfully"); + + // verify CompactPoints + if (compactPoints) + { + vtkm::Id numExpectedPoints = ((2 * pointDims[0] * pointDims[1]) // x-y faces + + (2 * pointDims[0] * pointDims[2]) // x-z faces + + (2 * pointDims[1] * pointDims[2]) // y-z faces + - (4 * pointDims[0]) // overcounted x edges + - (4 * pointDims[1]) // overcounted y edges + - (4 * pointDims[2]) // overcounted z edges + + 8); // undercounted corners + vtkm::Id numOutputPoints = resultds.GetNumberOfPoints(); + VTKM_TEST_ASSERT(numOutputPoints == numExpectedPoints); + } + else + { + VTKM_TEST_ASSERT(resultds.GetNumberOfPoints() == ds.GetNumberOfPoints()); + } +} + +void TestWithUniformGrid() +{ + std::cout << "Testing with uniform grid\n"; + vtkm::cont::DataSet ds = MakeUniformDataTestSet(); + std::cout << "Compact Points Off\n"; + TestExternalFacesStructuredGrid(ds, false); + std::cout << "Compact Points On\n"; + TestExternalFacesStructuredGrid(ds, true); +} + +void TestWithCurvilinearGrid() +{ + std::cout << "Testing with curvilinear grid\n"; + vtkm::cont::DataSet ds = MakeCurvilinearDataTestSet(); + std::cout << "Compact Points Off\n"; + TestExternalFacesStructuredGrid(ds, false); + std::cout << "Compact Points On\n"; + TestExternalFacesStructuredGrid(ds, true); +} + void TestExternalFacesFilter() { TestWithHeterogeneousMesh(); @@ -142,6 +219,8 @@ void TestExternalFacesFilter() TestWithUniformMesh(); TestWithRectilinearMesh(); TestWithMixed2Dand3DMesh(); + TestWithUniformGrid(); + TestWithCurvilinearGrid(); } } // anonymous namespace diff --git a/vtkm/filter/entity_extraction/worklet/ExternalFaces.h b/vtkm/filter/entity_extraction/worklet/ExternalFaces.h index 8ccb9a9a85..7251cd4f27 100644 --- a/vtkm/filter/entity_extraction/worklet/ExternalFaces.h +++ b/vtkm/filter/entity_extraction/worklet/ExternalFaces.h @@ -43,269 +43,139 @@ namespace worklet struct ExternalFaces { - //Worklet that returns the number of external faces for each structured cell - class NumExternalFacesPerStructuredCell : public vtkm::worklet::WorkletVisitCellsWithPoints + struct ExtractStructuredFace : public vtkm::worklet::WorkletMapField { - public: - using ControlSignature = void(CellSetIn inCellSet, - FieldOut numFacesInCell, - FieldInPoint pointCoordinates); - using ExecutionSignature = _2(CellShape, _3); - using InputDomain = _1; - - VTKM_CONT - NumExternalFacesPerStructuredCell(const vtkm::Vec3f_64& min_point, - const vtkm::Vec3f_64& max_point) - : MinPoint(min_point) - , MaxPoint(max_point) + using ControlSignature = void(FieldIn indices, FieldOut cellConnections, FieldOut cellMap); + + vtkm::Id3 CellDimensions; + vtkm::Id3 PointDimensions; + vtkm::Id XYCellSize; + vtkm::Id XZCellSize; + vtkm::Id YZCellSize; + vtkm::Id XYPointSize; + vtkm::Id XZPointSize; + vtkm::Id YZPointSize; + + VTKM_CONT ExtractStructuredFace(vtkm::Id3 cellDimensions) + : CellDimensions(cellDimensions) + , PointDimensions(cellDimensions[0] + 1, cellDimensions[1] + 1, cellDimensions[2] + 1) + , XYCellSize(CellDimensions[0] * CellDimensions[1]) + , XZCellSize(CellDimensions[0] * CellDimensions[2]) + , YZCellSize(CellDimensions[1] * CellDimensions[2]) + , XYPointSize(PointDimensions[0] * PointDimensions[1]) + , XZPointSize(PointDimensions[0] * PointDimensions[2]) + , YZPointSize(PointDimensions[1] * PointDimensions[2]) { } - VTKM_EXEC - static inline vtkm::IdComponent CountExternalFacesOnDimension(vtkm::Float64 grid_min, - vtkm::Float64 grid_max, - vtkm::Float64 cell_min, - vtkm::Float64 cell_max) + VTKM_EXEC void operator()(vtkm::Id index, vtkm::Id4& connections, vtkm::Id& cellMap) const { - vtkm::IdComponent count = 0; - - bool cell_min_at_grid_boundary = cell_min <= grid_min; - bool cell_max_at_grid_boundary = cell_max >= grid_max; - - if (cell_min_at_grid_boundary && !cell_max_at_grid_boundary) + if (index < this->XYCellSize) { - count++; + this->GetXYLowCell(index, connections, cellMap); + return; } - else if (!cell_min_at_grid_boundary && cell_max_at_grid_boundary) + index -= this->XYCellSize; + if (index < this->XYCellSize) { - count++; + this->GetXYHighCell(index, connections, cellMap); + return; } - else if (cell_min_at_grid_boundary && cell_max_at_grid_boundary) + index -= this->XYCellSize; + if (index < this->XZCellSize) { - count += 2; + this->GetXZLowCell(index, connections, cellMap); + return; } - - return count; + index -= this->XZCellSize; + if (index < this->XZCellSize) + { + this->GetXZHighCell(index, connections, cellMap); + return; + } + index -= this->XZCellSize; + if (index < this->YZCellSize) + { + this->GetYZLowCell(index, connections, cellMap); + return; + } + index -= this->YZCellSize; + VTKM_ASSERT(index < this->YZCellSize); + this->GetYZHighCell(index, connections, cellMap); } - template - VTKM_EXEC vtkm::IdComponent operator()(CellShapeTag shape, - const PointCoordVecType& pointCoordinates) const + VTKM_EXEC void GetXYLowCell(vtkm::Id index, vtkm::Id4& connections, vtkm::Id& cellMap) const { - (void)shape; // C4100 false positive workaround - VTKM_ASSERT(shape.Id == CELL_SHAPE_HEXAHEDRON); - - vtkm::IdComponent count = 0; - - count += CountExternalFacesOnDimension( - MinPoint[0], MaxPoint[0], pointCoordinates[0][0], pointCoordinates[1][0]); - - count += CountExternalFacesOnDimension( - MinPoint[1], MaxPoint[1], pointCoordinates[0][1], pointCoordinates[3][1]); - - count += CountExternalFacesOnDimension( - MinPoint[2], MaxPoint[2], pointCoordinates[0][2], pointCoordinates[4][2]); - - return count; + vtkm::Id2 cellIndex{ index % this->CellDimensions[0], index / this->CellDimensions[0] }; + vtkm::Id pointIndex = cellIndex[0] + (cellIndex[1] * this->PointDimensions[0]); + connections[0] = pointIndex; + connections[1] = pointIndex + this->PointDimensions[0]; + connections[2] = pointIndex + this->PointDimensions[0] + 1; + connections[3] = pointIndex + 1; + cellMap = index; } - private: - vtkm::Vec3f_64 MinPoint; - vtkm::Vec3f_64 MaxPoint; - }; - - - //Worklet that finds face connectivity for each structured cell - class BuildConnectivityStructured : public vtkm::worklet::WorkletVisitCellsWithPoints - { - public: - using ControlSignature = void(CellSetIn inCellSet, - WholeCellSetIn<> inputCell, - FieldOut faceShapes, - FieldOut facePointCount, - FieldOut faceConnectivity, - FieldInPoint pointCoordinates); - using ExecutionSignature = void(CellShape, VisitIndex, InputIndex, _2, _3, _4, _5, _6); - using InputDomain = _1; - - using ScatterType = vtkm::worklet::ScatterCounting; - - VTKM_CONT - BuildConnectivityStructured(const vtkm::Vec3f_64& min_point, const vtkm::Vec3f_64& max_point) - : MinPoint(min_point) - , MaxPoint(max_point) + VTKM_EXEC void GetXYHighCell(vtkm::Id index, vtkm::Id4& connections, vtkm::Id& cellMap) const { + vtkm::Id2 cellIndex{ index % this->CellDimensions[0], index / this->CellDimensions[0] }; + vtkm::Id offset = this->XYPointSize * (this->PointDimensions[2] - 1); + vtkm::Id pointIndex = offset + cellIndex[0] + (cellIndex[1] * this->PointDimensions[0]); + connections[0] = pointIndex; + connections[1] = pointIndex + 1; + connections[2] = pointIndex + this->PointDimensions[0] + 1; + connections[3] = pointIndex + this->PointDimensions[0]; + cellMap = this->XYCellSize * (this->CellDimensions[2] - 1) + index; } - enum FaceType - { - FACE_GRID_MIN, - FACE_GRID_MAX, - FACE_GRID_MIN_AND_MAX, - FACE_NONE - }; - - VTKM_EXEC - static inline bool FoundFaceOnDimension(vtkm::Float64 grid_min, - vtkm::Float64 grid_max, - vtkm::Float64 cell_min, - vtkm::Float64 cell_max, - vtkm::IdComponent& faceIndex, - vtkm::IdComponent& count, - vtkm::IdComponent dimensionFaceOffset, - vtkm::IdComponent visitIndex) + VTKM_EXEC void GetXZLowCell(vtkm::Id index, vtkm::Id4& connections, vtkm::Id& cellMap) const { - bool cell_min_at_grid_boundary = cell_min <= grid_min; - bool cell_max_at_grid_boundary = cell_max >= grid_max; - - FaceType Faces = FaceType::FACE_NONE; - - if (cell_min_at_grid_boundary && !cell_max_at_grid_boundary) - { - Faces = FaceType::FACE_GRID_MIN; - } - else if (!cell_min_at_grid_boundary && cell_max_at_grid_boundary) - { - Faces = FaceType::FACE_GRID_MAX; - } - else if (cell_min_at_grid_boundary && cell_max_at_grid_boundary) - { - Faces = FaceType::FACE_GRID_MIN_AND_MAX; - } - - if (Faces == FaceType::FACE_NONE) - return false; - - if (Faces == FaceType::FACE_GRID_MIN) - { - if (visitIndex == count) - { - faceIndex = dimensionFaceOffset; - return true; - } - else - { - count++; - } - } - else if (Faces == FaceType::FACE_GRID_MAX) - { - if (visitIndex == count) - { - faceIndex = dimensionFaceOffset + 1; - return true; - } - else - { - count++; - } - } - else if (Faces == FaceType::FACE_GRID_MIN_AND_MAX) - { - if (visitIndex == count) - { - faceIndex = dimensionFaceOffset; - return true; - } - count++; - if (visitIndex == count) - { - faceIndex = dimensionFaceOffset + 1; - return true; - } - count++; - } - - return false; + vtkm::Id2 cellIndex{ index % this->CellDimensions[0], index / this->CellDimensions[0] }; + vtkm::Id pointIndex = cellIndex[0] + (cellIndex[1] * this->XYPointSize); + connections[0] = pointIndex; + connections[1] = pointIndex + this->XYPointSize; + connections[2] = pointIndex + this->XYPointSize + 1; + connections[3] = pointIndex + 1; + cellMap = cellIndex[0] + (cellIndex[1] * this->XYCellSize); } - template - VTKM_EXEC inline vtkm::IdComponent FindFaceIndexForVisit( - vtkm::IdComponent visitIndex, - const PointCoordVecType& pointCoordinates) const + VTKM_EXEC void GetXZHighCell(vtkm::Id index, vtkm::Id4& connections, vtkm::Id& cellMap) const { - vtkm::IdComponent count = 0; - vtkm::IdComponent faceIndex = 0; - // Search X dimension - if (!FoundFaceOnDimension(MinPoint[0], - MaxPoint[0], - pointCoordinates[0][0], - pointCoordinates[1][0], - faceIndex, - count, - 0, - visitIndex)) - { - // Search Y dimension - if (!FoundFaceOnDimension(MinPoint[1], - MaxPoint[1], - pointCoordinates[0][1], - pointCoordinates[3][1], - faceIndex, - count, - 2, - visitIndex)) - { - // Search Z dimension - FoundFaceOnDimension(MinPoint[2], - MaxPoint[2], - pointCoordinates[0][2], - pointCoordinates[4][2], - faceIndex, - count, - 4, - visitIndex); - } - } - return faceIndex; + vtkm::Id2 cellIndex{ index % this->CellDimensions[0], index / this->CellDimensions[0] }; + vtkm::Id offset = this->XYPointSize - this->PointDimensions[0]; + vtkm::Id pointIndex = offset + cellIndex[0] + (cellIndex[1] * this->XYPointSize); + connections[0] = pointIndex; + connections[1] = pointIndex + 1; + connections[2] = pointIndex + this->XYPointSize + 1; + connections[3] = pointIndex + this->XYPointSize; + cellMap = this->XYCellSize - this->CellDimensions[0] + cellIndex[0] + + (cellIndex[1] * this->XYCellSize); } - template - VTKM_EXEC void operator()(CellShapeTag shape, - vtkm::IdComponent visitIndex, - vtkm::Id inputIndex, - const CellSetType& cellSet, - vtkm::UInt8& shapeOut, - vtkm::IdComponent& numFacePointsOut, - ConnectivityType& faceConnectivity, - const PointCoordVecType& pointCoordinates) const + VTKM_EXEC void GetYZLowCell(vtkm::Id index, vtkm::Id4& connections, vtkm::Id& cellMap) const { - VTKM_ASSERT(shape.Id == CELL_SHAPE_HEXAHEDRON); - - vtkm::IdComponent faceIndex = FindFaceIndexForVisit(visitIndex, pointCoordinates); - - vtkm::IdComponent numFacePoints; - vtkm::exec::CellFaceNumberOfPoints(faceIndex, shape, numFacePoints); - VTKM_ASSERT(numFacePoints == faceConnectivity.GetNumberOfComponents()); - - typename CellSetType::IndicesType inCellIndices = cellSet.GetIndices(inputIndex); - - shapeOut = vtkm::CELL_SHAPE_QUAD; - numFacePointsOut = 4; - - for (vtkm::IdComponent facePointIndex = 0; facePointIndex < numFacePoints; facePointIndex++) - { - vtkm::IdComponent localFaceIndex; - vtkm::ErrorCode status = - vtkm::exec::CellFaceLocalIndex(facePointIndex, faceIndex, shape, localFaceIndex); - if (status == vtkm::ErrorCode::Success) - { - faceConnectivity[facePointIndex] = inCellIndices[localFaceIndex]; - } - else - { - // An error condition, but do we want to crash the operation? - faceConnectivity[facePointIndex] = 0; - } - } + vtkm::Id2 cellIndex{ index % this->CellDimensions[1], index / this->CellDimensions[1] }; + vtkm::Id pointIndex = + (cellIndex[0] * this->PointDimensions[0]) + (cellIndex[1] * this->XYPointSize); + connections[0] = pointIndex; + connections[1] = pointIndex + this->XYPointSize; + connections[2] = pointIndex + this->XYPointSize + this->PointDimensions[0]; + connections[3] = pointIndex + this->PointDimensions[0]; + cellMap = (cellIndex[0] * this->CellDimensions[0]) + (cellIndex[1] * this->XYCellSize); } - private: - vtkm::Vec3f_64 MinPoint; - vtkm::Vec3f_64 MaxPoint; + VTKM_EXEC void GetYZHighCell(vtkm::Id index, vtkm::Id4& connections, vtkm::Id& cellMap) const + { + vtkm::Id2 cellIndex{ index % this->CellDimensions[1], index / this->CellDimensions[1] }; + vtkm::Id offset = this->PointDimensions[0] - 1; + vtkm::Id pointIndex = + offset + (cellIndex[0] * this->PointDimensions[0]) + (cellIndex[1] * this->XYPointSize); + connections[0] = pointIndex; + connections[1] = pointIndex + this->PointDimensions[0]; + connections[2] = pointIndex + this->XYPointSize + this->PointDimensions[0]; + connections[3] = pointIndex + this->XYPointSize; + cellMap = (this->CellDimensions[0] - 1) + (cellIndex[0] * this->CellDimensions[0]) + + (cellIndex[1] * this->XYCellSize); + } }; // Worklet that returns the number of faces for each cell/shape @@ -754,6 +624,17 @@ public: T Bias; }; + VTKM_CONT vtkm::cont::CellSetExplicit<> MakeCellSetExplicit( + vtkm::Id numPoints, + const vtkm::cont::ArrayHandle& shapes, + const vtkm::cont::ArrayHandle& connectivity, + const vtkm::cont::ArrayHandle& offsets) + { + vtkm::cont::CellSetExplicit<> outCellSet; + outCellSet.Fill(numPoints, shapes, connectivity, offsets); + return outCellSet; + } + public: VTKM_CONT ExternalFaces() @@ -775,102 +656,39 @@ public: /// /// Faster Run() method for uniform and rectilinear grid types. /// Uses grid extents to find cells on the boundaries of the grid. - template - VTKM_CONT void Run( - const vtkm::cont::CellSetStructured<3>& inCellSet, - const vtkm::cont::CoordinateSystem& coord, - vtkm::cont::CellSetExplicit& outCellSet) + VTKM_CONT vtkm::cont::CellSetSingleType<> Run(const vtkm::cont::CellSetStructured<3>& inCellSet) { // create an invoker vtkm::cont::Invoker invoke; - vtkm::Vec3f_64 MinPoint; - vtkm::Vec3f_64 MaxPoint; - - vtkm::Id3 PointDimensions = inCellSet.GetPointDimensions(); - - using DefaultHandle = vtkm::cont::ArrayHandle; - using CartesianArrayHandle = - vtkm::cont::ArrayHandleCartesianProduct; - - auto coordData = coord.GetData(); - if (coordData.CanConvert()) - { - const auto vertices = coordData.AsArrayHandle(); - const auto vertsSize = vertices.GetNumberOfValues(); - const auto tmp = vtkm::cont::ArrayGetValues({ 0, vertsSize - 1 }, vertices); - MinPoint = tmp[0]; - MaxPoint = tmp[1]; - } - else - { - auto vertices = coordData.AsArrayHandle(); - auto Coordinates = vertices.ReadPortal(); - - MinPoint = Coordinates.GetOrigin(); - vtkm::Vec3f_64 spacing = Coordinates.GetSpacing(); - - vtkm::Vec3f_64 unitLength; - unitLength[0] = static_cast(PointDimensions[0] - 1); - unitLength[1] = static_cast(PointDimensions[1] - 1); - unitLength[2] = static_cast(PointDimensions[2] - 1); - MaxPoint = MinPoint + spacing * unitLength; - } - - // Count the number of external faces per cell - vtkm::cont::ArrayHandle numExternalFaces; - invoke(NumExternalFacesPerStructuredCell(MinPoint, MaxPoint), - inCellSet, - numExternalFaces, - coordData); - - vtkm::Id numberOfExternalFaces = - vtkm::cont::Algorithm::Reduce(numExternalFaces, 0, vtkm::Sum()); + vtkm::Id3 cellDimensions = inCellSet.GetCellDimensions(); - vtkm::worklet::ScatterCounting scatterCellToExternalFace(numExternalFaces); + ExtractStructuredFace worklet{ cellDimensions }; - // Maps output cells to input cells. Store this for cell field mapping. - this->CellIdMap = scatterCellToExternalFace.GetOutputToInputMap(); + vtkm::Id numOutCells = + (2 * worklet.XYCellSize) + (2 * worklet.XZCellSize) + (2 * worklet.YZCellSize); - numExternalFaces.ReleaseResources(); - - vtkm::Id connectivitySize = 4 * numberOfExternalFaces; - vtkm::cont::ArrayHandle faceConnectivity; - vtkm::cont::ArrayHandle faceShapes; - vtkm::cont::ArrayHandle facePointCount; - // Must pre allocate because worklet invocation will not have enough - // information to. - faceConnectivity.Allocate(connectivitySize); - - // Build connectivity for external faces - invoke(BuildConnectivityStructured(MinPoint, MaxPoint), - scatterCellToExternalFace, - inCellSet, - inCellSet, - faceShapes, - facePointCount, - vtkm::cont::make_ArrayHandleGroupVec<4>(faceConnectivity), - coordData); + vtkm::cont::ArrayHandle connections; - auto offsets = vtkm::cont::ConvertNumComponentsToOffsets(facePointCount); + invoke(worklet, + vtkm::cont::ArrayHandleIndex(numOutCells), + vtkm::cont::make_ArrayHandleGroupVec<4>(connections), + this->CellIdMap); - outCellSet.Fill(inCellSet.GetNumberOfPoints(), faceShapes, faceConnectivity, offsets); + vtkm::cont::CellSetSingleType<> outCellSet; + outCellSet.Fill(inCellSet.GetNumberOfPoints(), vtkm::CELL_SHAPE_QUAD, 4, connections); + return outCellSet; } /////////////////////////////////////////////////// /// \brief ExternalFaces: Extract Faces on outside of geometry - template - VTKM_CONT void Run( - const InCellSetType& inCellSet, - vtkm::cont::CellSetExplicit& outCellSet) + template + VTKM_CONT vtkm::cont::CellSetExplicit<> Run(const InCellSetType& inCellSet) { using PointCountArrayType = vtkm::cont::ArrayHandle; - using ShapeArrayType = vtkm::cont::ArrayHandle; - using OffsetsArrayType = vtkm::cont::ArrayHandle; - using ConnectivityArrayType = vtkm::cont::ArrayHandle; + using ShapeArrayType = vtkm::cont::ArrayHandle; + using OffsetsArrayType = vtkm::cont::ArrayHandle; + using ConnectivityArrayType = vtkm::cont::ArrayHandle; // create an invoker vtkm::cont::Invoker invoke; @@ -928,17 +746,17 @@ public: if (!polyDataConnectivitySize) { // Data has no faces. Output is empty. + vtkm::cont::CellSetExplicit<> outCellSet; outCellSet.PrepareToAddCells(0, 0); outCellSet.CompleteAddingCells(inCellSet.GetNumberOfPoints()); - return; + return outCellSet; } else { // Pass only input poly data to output - outCellSet.Fill( - inCellSet.GetNumberOfPoints(), polyDataShapes, polyDataConnectivity, polyDataOffsets); this->CellIdMap = polyDataCellIdMap; - return; + return MakeCellSetExplicit( + inCellSet.GetNumberOfPoints(), polyDataShapes, polyDataConnectivity, polyDataOffsets); } } @@ -1040,11 +858,11 @@ public: if (!polyDataConnectivitySize) { - outCellSet.Fill(inCellSet.GetNumberOfPoints(), - externalFacesShapes, - externalFacesConnectivity, - pointsPerExternalFaceOffsets); this->CellIdMap = faceToCellIdMap; + return MakeCellSetExplicit(inCellSet.GetNumberOfPoints(), + externalFacesShapes, + externalFacesConnectivity, + pointsPerExternalFaceOffsets); } else { @@ -1085,9 +903,9 @@ public: vtkm::cont::ArrayHandle joinedCellIdMap; vtkm::cont::ArrayCopy(cellIdMapArray, joinedCellIdMap); - outCellSet.Fill( - inCellSet.GetNumberOfPoints(), joinedShapesArray, joinedConnectivity, joinedOffsets); this->CellIdMap = joinedCellIdMap; + return MakeCellSetExplicit( + inCellSet.GetNumberOfPoints(), joinedShapesArray, joinedConnectivity, joinedOffsets); } } -- GitLab