diff --git a/docs/changelog/tetrahedralize-triangulate-check.md b/docs/changelog/tetrahedralize-triangulate-check.md new file mode 100644 index 000000000..de95e8b6c --- /dev/null +++ b/docs/changelog/tetrahedralize-triangulate-check.md @@ -0,0 +1,3 @@ +# Tetrahedralize and Triangulate filters now check if the input is already tetrahedral/triangular + +Previously, tetrahedralize/triangulate would blindly convert all the cells to tetrahedra/triangles, even when they were already. Now, the dataset is directly returned if the CellSet is a CellSetSingleType of tetras/triangles, and no further processing is done in the worklets for CellSetExplicit when all shapes are tetras or triangles. diff --git a/vtkm/filter/geometry_refinement/Tetrahedralize.cxx b/vtkm/filter/geometry_refinement/Tetrahedralize.cxx index bec934e50..979eae5c1 100644 --- a/vtkm/filter/geometry_refinement/Tetrahedralize.cxx +++ b/vtkm/filter/geometry_refinement/Tetrahedralize.cxx @@ -8,6 +8,7 @@ // PURPOSE. See the above copyright notice for more information. //============================================================================ +#include #include #include #include @@ -41,6 +42,18 @@ VTKM_CONT bool DoMapField(vtkm::cont::DataSet& result, return false; } } + +struct IsShapeTetra +{ + VTKM_EXEC_CONT + bool operator()(vtkm::UInt8 shape) const { return shape == vtkm::CELL_SHAPE_TETRA; } +}; + +struct BinaryAnd +{ + VTKM_EXEC_CONT + bool operator()(bool u, bool v) const { return u && v; } +}; } // anonymous namespace namespace vtkm @@ -53,14 +66,58 @@ VTKM_CONT vtkm::cont::DataSet Tetrahedralize::DoExecute(const vtkm::cont::DataSe { const vtkm::cont::UnknownCellSet& inCellSet = input.GetCellSet(); - vtkm::cont::CellSetSingleType<> outCellSet; - vtkm::worklet::Tetrahedralize worklet; - vtkm::cont::CastAndCall(inCellSet, - [&](const auto& concrete) { outCellSet = worklet.Run(concrete); }); + // In case we already have a CellSetSingleType of tetras, + // don't call the worklet and return the input DataSet directly + if (inCellSet.CanConvert>() && + inCellSet.AsCellSet>().GetCellShapeAsId() == + vtkm::CellShapeTagTetra::Id) + { + return input; + } - auto mapper = [&](auto& result, const auto& f) { DoMapField(result, f, worklet); }; - // create the output dataset (without a CoordinateSystem). - vtkm::cont::DataSet output = this->CreateResult(input, outCellSet, mapper); + vtkm::cont::CellSetSingleType<> outCellSet; + vtkm::cont::DataSet output; + + // Optimization in case we only have tetras in the CellSet + bool allTetras = false; + if (inCellSet.CanConvert>()) + { + vtkm::cont::CellSetExplicit<> inCellSetExplicit = + inCellSet.AsCellSet>(); + + auto shapeArray = inCellSetExplicit.GetShapesArray(vtkm::TopologyElementTagCell(), + vtkm::TopologyElementTagPoint()); + auto isCellTetraArray = vtkm::cont::make_ArrayHandleTransform(shapeArray, IsShapeTetra{}); + + allTetras = vtkm::cont::Algorithm::Reduce(isCellTetraArray, true, BinaryAnd{}); + + if (allTetras) + { + // Reuse the input's connectivity array + outCellSet.Fill(inCellSet.GetNumberOfPoints(), + vtkm::CellShapeTagTetra::Id, + 4, + inCellSetExplicit.GetConnectivityArray(vtkm::TopologyElementTagCell(), + vtkm::TopologyElementTagPoint())); + + // Copy all fields from the input + output = this->CreateResult(input, outCellSet, [&](auto& result, const auto& f) { + result.AddField(f); + return true; + }); + } + } + + if (!allTetras) + { + vtkm::worklet::Tetrahedralize worklet; + vtkm::cont::CastAndCall(inCellSet, + [&](const auto& concrete) { outCellSet = worklet.Run(concrete); }); + + auto mapper = [&](auto& result, const auto& f) { DoMapField(result, f, worklet); }; + // create the output dataset (without a CoordinateSystem). + output = this->CreateResult(input, outCellSet, mapper); + } // We did not change the geometry of the input dataset at all. Just attach coordinate system // of input dataset to output dataset. diff --git a/vtkm/filter/geometry_refinement/Triangulate.cxx b/vtkm/filter/geometry_refinement/Triangulate.cxx index e8eba2509..913c485a6 100644 --- a/vtkm/filter/geometry_refinement/Triangulate.cxx +++ b/vtkm/filter/geometry_refinement/Triangulate.cxx @@ -8,6 +8,7 @@ // PURPOSE. See the above copyright notice for more information. //============================================================================ +#include #include #include #include @@ -42,6 +43,18 @@ VTKM_CONT bool DoMapField(vtkm::cont::DataSet& result, return false; } } + +struct IsShapeTriangle +{ + VTKM_EXEC_CONT + bool operator()(vtkm::UInt8 shape) const { return shape == vtkm::CELL_SHAPE_TRIANGLE; } +}; + +struct BinaryAnd +{ + VTKM_EXEC_CONT + bool operator()(bool u, bool v) const { return u && v; } +}; } // anonymous namespace namespace vtkm @@ -54,15 +67,57 @@ VTKM_CONT vtkm::cont::DataSet Triangulate::DoExecute(const vtkm::cont::DataSet& { const vtkm::cont::UnknownCellSet& inCellSet = input.GetCellSet(); + // In case we already have a CellSetSingleType of tetras, + // don't call the worklet and return the input DataSet directly + if (inCellSet.CanConvert>() && + inCellSet.AsCellSet>().GetCellShapeAsId() == + vtkm::CellShapeTagTriangle::Id) + { + return input; + } + vtkm::cont::CellSetSingleType<> outCellSet; - vtkm::worklet::Triangulate worklet; + vtkm::cont::DataSet output; - vtkm::cont::CastAndCall(inCellSet, - [&](const auto& concrete) { outCellSet = worklet.Run(concrete); }); + // Optimization in case we only have triangles in the CellSet + bool allTriangles = false; + if (inCellSet.CanConvert>()) + { + vtkm::cont::CellSetExplicit<> inCellSetExplicit = + inCellSet.AsCellSet>(); - auto mapper = [&](auto& result, const auto& f) { DoMapField(result, f, worklet); }; - // create the output dataset (without a CoordinateSystem). - vtkm::cont::DataSet output = this->CreateResult(input, outCellSet, mapper); + auto shapeArray = inCellSetExplicit.GetShapesArray(vtkm::TopologyElementTagCell(), + vtkm::TopologyElementTagPoint()); + auto isCellTriangleArray = vtkm::cont::make_ArrayHandleTransform(shapeArray, IsShapeTriangle{}); + allTriangles = vtkm::cont::Algorithm::Reduce(isCellTriangleArray, true, BinaryAnd{}); + + if (allTriangles) + { + // Reuse the input's connectivity array + outCellSet.Fill(inCellSet.GetNumberOfPoints(), + vtkm::CellShapeTagTriangle::Id, + 3, + inCellSetExplicit.GetConnectivityArray(vtkm::TopologyElementTagCell(), + vtkm::TopologyElementTagPoint())); + + // Copy all fields from the input + output = this->CreateResult(input, outCellSet, [&](auto& result, const auto& f) { + result.AddField(f); + return true; + }); + } + } + + if (!allTriangles) + { + vtkm::worklet::Triangulate worklet; + vtkm::cont::CastAndCall(inCellSet, + [&](const auto& concrete) { outCellSet = worklet.Run(concrete); }); + + auto mapper = [&](auto& result, const auto& f) { DoMapField(result, f, worklet); }; + // create the output dataset (without a CoordinateSystem). + output = this->CreateResult(input, outCellSet, mapper); + } // We did not change the geometry of the input dataset at all. Just attach coordinate system // of input dataset to output dataset. diff --git a/vtkm/filter/geometry_refinement/testing/UnitTestTetrahedralizeFilter.cxx b/vtkm/filter/geometry_refinement/testing/UnitTestTetrahedralizeFilter.cxx index cae7c2413..0e24e868f 100644 --- a/vtkm/filter/geometry_refinement/testing/UnitTestTetrahedralizeFilter.cxx +++ b/vtkm/filter/geometry_refinement/testing/UnitTestTetrahedralizeFilter.cxx @@ -8,6 +8,7 @@ // PURPOSE. See the above copyright notice for more information. //============================================================================ +#include #include #include @@ -67,10 +68,56 @@ public: VTKM_TEST_ASSERT(outData.ReadPortal().Get(10) == 130.5f, "Wrong cell field data"); } + void TestCellSetSingleTypeTetra() const + { + vtkm::cont::DataSet dataset; + vtkm::cont::CellSetSingleType<> cellSet; + + auto connectivity = vtkm::cont::make_ArrayHandle({ 0, 1, 2, 3, 3, 2, 1, 4 }); + cellSet.Fill(5, vtkm::CELL_SHAPE_TETRA, 4, connectivity); + + dataset.SetCellSet(cellSet); + + vtkm::filter::geometry_refinement::Tetrahedralize tetrahedralize; + vtkm::cont::DataSet output = tetrahedralize.Execute(dataset); + + VTKM_TEST_ASSERT(dataset.GetCellSet().GetCellSetBase() == output.GetCellSet().GetCellSetBase(), + "Pointer to the CellSetSingleType has changed."); + } + + void TestCellSetExplicitTetra() const + { + std::vector coords{ + vtkm::Vec3f_32(0.0f, 0.0f, 0.0f), vtkm::Vec3f_32(2.0f, 0.0f, 0.0f), + vtkm::Vec3f_32(2.0f, 4.0f, 0.0f), vtkm::Vec3f_32(0.0f, 4.0f, 0.0f), + vtkm::Vec3f_32(1.0f, 0.0f, 3.0f), + }; + std::vector shapes{ vtkm::CELL_SHAPE_TETRA, vtkm::CELL_SHAPE_TETRA }; + std::vector indices{ 4, 4 }; + std::vector connectivity{ 0, 1, 2, 3, 1, 2, 3, 4 }; + + vtkm::cont::DataSetBuilderExplicit dsb; + vtkm::cont::DataSet dataset = dsb.Create(coords, shapes, indices, connectivity); + + vtkm::filter::geometry_refinement::Tetrahedralize tetrahedralize; + vtkm::cont::DataSet output = tetrahedralize.Execute(dataset); + vtkm::cont::UnknownCellSet outputCellSet = output.GetCellSet(); + + VTKM_TEST_ASSERT(outputCellSet.IsType>(), + "Output CellSet is not CellSetSingleType"); + VTKM_TEST_ASSERT(output.GetNumberOfCells() == 2, "Wrong number of cells"); + VTKM_TEST_ASSERT(outputCellSet.GetCellShape(0) == vtkm::CellShapeTagTetra::Id, + "Cell is not tetra"); + VTKM_TEST_ASSERT(outputCellSet.GetCellShape(1) == vtkm::CellShapeTagTetra::Id, + "Cell is not tetra"); + } + void operator()() const { this->TestStructured(); this->TestExplicit(); + this->TestCellSetSingleTypeTetra(); + this->TestCellSetExplicitTetra(); } }; } diff --git a/vtkm/filter/geometry_refinement/testing/UnitTestTriangulateFilter.cxx b/vtkm/filter/geometry_refinement/testing/UnitTestTriangulateFilter.cxx index c7f2c5cf2..7a167710e 100644 --- a/vtkm/filter/geometry_refinement/testing/UnitTestTriangulateFilter.cxx +++ b/vtkm/filter/geometry_refinement/testing/UnitTestTriangulateFilter.cxx @@ -8,6 +8,7 @@ // PURPOSE. See the above copyright notice for more information. //============================================================================ +#include #include #include @@ -61,10 +62,55 @@ public: VTKM_TEST_ASSERT(outData.ReadPortal().Get(6) == 3.f, "Wrong cell field data"); } + void TestCellSetSingleTypeTriangle() const + { + vtkm::cont::DataSet dataset; + vtkm::cont::CellSetSingleType<> cellSet; + + auto connectivity = vtkm::cont::make_ArrayHandle({ 0, 1, 2, 1, 2, 3 }); + cellSet.Fill(4, vtkm::CELL_SHAPE_TRIANGLE, 3, connectivity); + + dataset.SetCellSet(cellSet); + + vtkm::filter::geometry_refinement::Triangulate triangulate; + vtkm::cont::DataSet output = triangulate.Execute(dataset); + + VTKM_TEST_ASSERT(dataset.GetCellSet().GetCellSetBase() == output.GetCellSet().GetCellSetBase(), + "Pointer to the CellSetSingleType has changed."); + } + + void TestCellSetExplicitTriangle() const + { + std::vector coords{ vtkm::Vec3f_32(0.0f, 0.0f, 0.0f), + vtkm::Vec3f_32(2.0f, 0.0f, 0.0f), + vtkm::Vec3f_32(2.0f, 4.0f, 0.0f), + vtkm::Vec3f_32(0.0f, 4.0f, 0.0f) }; + std::vector shapes{ vtkm::CELL_SHAPE_TRIANGLE, vtkm::CELL_SHAPE_TRIANGLE }; + std::vector indices{ 3, 3 }; + std::vector connectivity{ 0, 1, 2, 1, 2, 3 }; + + vtkm::cont::DataSetBuilderExplicit dsb; + vtkm::cont::DataSet dataset = dsb.Create(coords, shapes, indices, connectivity); + + vtkm::filter::geometry_refinement::Triangulate triangulate; + vtkm::cont::DataSet output = triangulate.Execute(dataset); + vtkm::cont::UnknownCellSet outputCellSet = output.GetCellSet(); + + VTKM_TEST_ASSERT(outputCellSet.IsType>(), + "Output CellSet is not CellSetSingleType"); + VTKM_TEST_ASSERT(output.GetNumberOfCells() == 2, "Wrong number of cells"); + VTKM_TEST_ASSERT(outputCellSet.GetCellShape(0) == vtkm::CellShapeTagTriangle::Id, + "Cell is not triangular"); + VTKM_TEST_ASSERT(outputCellSet.GetCellShape(1) == vtkm::CellShapeTagTriangle::Id, + "Cell is not triangular"); + } + void operator()() const { this->TestStructured(); this->TestExplicit(); + this->TestCellSetSingleTypeTriangle(); + this->TestCellSetExplicitTriangle(); } }; }