Tetrahedralize and Triangulate filters check if input is already triangulated

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.

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.

vtkm/filter/geometry_refinement/Tetrahedralize.cxx

@@ -8,6 +8,7 @@
 //  PURPOSE.  See the above copyright notice for more information.
 //============================================================================
 
+#include <vtkm/cont/Algorithm.h>
 #include <vtkm/filter/MapFieldPermutation.h>
 #include <vtkm/filter/geometry_refinement/Tetrahedralize.h>
 #include <vtkm/filter/geometry_refinement/worklet/Tetrahedralize.h>
@@ -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<vtkm::cont::CellSetSingleType<>>() &&
+      inCellSet.AsCellSet<vtkm::cont::CellSetSingleType<>>().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<>>())
+  {
+    vtkm::cont::CellSetExplicit<> inCellSetExplicit =
+      inCellSet.AsCellSet<vtkm::cont::CellSetExplicit<>>();
+
+    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.

vtkm/filter/geometry_refinement/Triangulate.cxx

@@ -8,6 +8,7 @@
 //  PURPOSE.  See the above copyright notice for more information.
 //============================================================================
 
+#include <vtkm/cont/Algorithm.h>
 #include <vtkm/filter/MapFieldPermutation.h>
 #include <vtkm/filter/geometry_refinement/Triangulate.h>
 #include <vtkm/filter/geometry_refinement/worklet/Triangulate.h>
@@ -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<vtkm::cont::CellSetSingleType<>>() &&
+      inCellSet.AsCellSet<vtkm::cont::CellSetSingleType<>>().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<>>())
+  {
+    vtkm::cont::CellSetExplicit<> inCellSetExplicit =
+      inCellSet.AsCellSet<vtkm::cont::CellSetExplicit<>>();
 
-  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.

vtkm/filter/geometry_refinement/testing/UnitTestTetrahedralizeFilter.cxx

@@ -8,6 +8,7 @@
 //  PURPOSE.  See the above copyright notice for more information.
 //============================================================================
 
+#include <vtkm/cont/DataSetBuilderExplicit.h>
 #include <vtkm/cont/testing/MakeTestDataSet.h>
 #include <vtkm/cont/testing/Testing.h>
 
@@ -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<vtkm::Id>({ 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<vtkm::Vec3f_32> 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<vtkm::UInt8> shapes{ vtkm::CELL_SHAPE_TETRA, vtkm::CELL_SHAPE_TETRA };
+    std::vector<vtkm::IdComponent> indices{ 4, 4 };
+    std::vector<vtkm::Id> 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<vtkm::cont::CellSetSingleType<>>(),
+                     "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();
   }
 };
 }

vtkm/filter/geometry_refinement/testing/UnitTestTriangulateFilter.cxx

@@ -8,6 +8,7 @@
 //  PURPOSE.  See the above copyright notice for more information.
 //============================================================================
 
+#include <vtkm/cont/DataSetBuilderExplicit.h>
 #include <vtkm/cont/testing/MakeTestDataSet.h>
 #include <vtkm/cont/testing/Testing.h>
 
@@ -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<vtkm::Id>({ 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<vtkm::Vec3f_32> 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<vtkm::UInt8> shapes{ vtkm::CELL_SHAPE_TRIANGLE, vtkm::CELL_SHAPE_TRIANGLE };
+    std::vector<vtkm::IdComponent> indices{ 3, 3 };
+    std::vector<vtkm::Id> 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<vtkm::cont::CellSetSingleType<>>(),
+                     "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();
   }
 };
 }