Merge topic 'fix-mc-connectivity'

27002d020 Fix distributed contour tree for marching cubes without using full boundary

Acked-by: Kitware Robot <kwrobot@kitware.com>
Acked-by: Kenneth Moreland <morelandkd@ornl.gov>
Merge-request: !3218

examples/contour_tree_distributed/ContourTreeApp.cxx

@@ -228,14 +228,6 @@ int main(int argc, char* argv[])
   if (parser.hasOption("--mc"))
   {
     useMarchingCubes = true;
-    if (useBoundaryExtremaOnly)
-    {
-      VTKM_LOG_S(vtkm::cont::LogLevel::Warn,
-                 "Warning: Marching cubes connectivity currently only works when "
-                 "using full boundary. Enabling the --useFullBoundary option "
-                 "to ensure that the app works.");
-      useBoundaryExtremaOnly = false;
-    }
   }
   bool preSplitFiles = false;
   if (parser.hasOption("--preSplitFiles"))

vtkm/filter/scalar_topology/testing/TestingContourTreeUniformDistributedFilter.h

@@ -289,9 +289,7 @@ inline vtkm::cont::PartitionedDataSet RunContourTreeDUniformDistributed(
   }
 
   filter.SetUseMarchingCubes(useMarchingCubes);
-  // Freudenthal: Only use boundary extrema; MC: use all points on boundary
+  filter.SetUseBoundaryExtremaOnly(true);
-  // TODO/FIXME: Figure out why MC does not work when only using boundary extrema
-  filter.SetUseBoundaryExtremaOnly(!useMarchingCubes);
   filter.SetAugmentHierarchicalTree(augmentHierarchicalTree);
   filter.SetActiveField(fieldName);
   auto result = filter.Execute(pds);

vtkm/filter/scalar_topology/worklet/contourtree_augmented/meshtypes/DataSetMeshTriangulation3DFreudenthal.h

@@ -89,11 +89,11 @@ public:
 
   DataSetMeshTriangulation3DFreudenthal(vtkm::Id3 meshSize);
 
-  MeshBoundary3DExec GetMeshBoundaryExecutionObject() const;
+  MeshBoundary3DExec<false> GetMeshBoundaryExecutionObject() const;
 
   void GetBoundaryVertices(IdArrayType& boundaryVertexArray,    // output
                            IdArrayType& boundarySortIndexArray, // output
-                           MeshBoundary3DExec* meshBoundaryExecObj =
+                           MeshBoundary3DExec<false>* meshBoundaryExecObj =
                              NULL // optional input, included for consistency with ContourTreeMesh
   ) const;
 
@@ -148,16 +148,16 @@ inline MeshStructureFreudenthal3D DataSetMeshTriangulation3DFreudenthal::Prepare
                                     token);
 }
 
-inline MeshBoundary3DExec DataSetMeshTriangulation3DFreudenthal::GetMeshBoundaryExecutionObject()
+inline MeshBoundary3DExec<false>
-  const
+DataSetMeshTriangulation3DFreudenthal::GetMeshBoundaryExecutionObject() const
 {
-  return MeshBoundary3DExec(this->MeshSize, this->SortIndices);
+  return MeshBoundary3DExec<false>(this->MeshSize, this->SortIndices);
 }
 
 inline void DataSetMeshTriangulation3DFreudenthal::GetBoundaryVertices(
-  IdArrayType& boundaryVertexArray,       // output
+  IdArrayType& boundaryVertexArray,              // output
-  IdArrayType& boundarySortIndexArray,    // output
+  IdArrayType& boundarySortIndexArray,           // output
-  MeshBoundary3DExec* meshBoundaryExecObj // input
+  MeshBoundary3DExec<false>* meshBoundaryExecObj // input
 ) const
 {
   vtkm::Id numBoundary = 2 * this->MeshSize[1] * this->MeshSize[0] // xy faces

vtkm/filter/scalar_topology/worklet/contourtree_augmented/meshtypes/DataSetMeshTriangulation3DMarchingCubes.h

@@ -92,12 +92,12 @@ public:
 
   DataSetMeshTriangulation3DMarchingCubes(vtkm::Id3 meshSize);
 
-  MeshBoundary3DExec GetMeshBoundaryExecutionObject() const;
+  MeshBoundary3DExec<true> GetMeshBoundaryExecutionObject() const;
 
   void GetBoundaryVertices(
     IdArrayType& boundaryVertexArray,    // output
     IdArrayType& boundarySortIndexArray, // output
-    MeshBoundary3DExec* meshBoundaryExecObj =
+    MeshBoundary3DExec<true>* meshBoundaryExecObj =
       nullptr // optional input, included for consistency with ContourTreeMesh
   ) const;
 
@@ -180,16 +180,16 @@ inline MeshStructureMarchingCubes DataSetMeshTriangulation3DMarchingCubes::Prepa
                                     token);
 }
 
-inline MeshBoundary3DExec DataSetMeshTriangulation3DMarchingCubes::GetMeshBoundaryExecutionObject()
+inline MeshBoundary3DExec<true>
-  const
+DataSetMeshTriangulation3DMarchingCubes::GetMeshBoundaryExecutionObject() const
 {
-  return MeshBoundary3DExec(this->MeshSize, this->SortOrder);
+  return MeshBoundary3DExec<true>(this->MeshSize, this->SortIndices);
 }
 
 inline void DataSetMeshTriangulation3DMarchingCubes::GetBoundaryVertices(
-  IdArrayType& boundaryVertexArray,       // output
+  IdArrayType& boundaryVertexArray,             // output
-  IdArrayType& boundarySortIndexArray,    // output
+  IdArrayType& boundarySortIndexArray,          // output
-  MeshBoundary3DExec* meshBoundaryExecObj // input
+  MeshBoundary3DExec<true>* meshBoundaryExecObj // input
 ) const
 {
   vtkm::Id numBoundary = 2 * this->MeshSize[1] * this->MeshSize[0] // xy faces

vtkm/filter/scalar_topology/worklet/contourtree_augmented/meshtypes/mesh_boundary/MeshBoundary3D.h

@@ -73,6 +73,10 @@ namespace worklet
 namespace contourtree_augmented
 {
 
+// TODO/FIXME: Consider making marching cubes connectivity its own class. However
+// at the moment making this a boolean template parameter makes it easuer to avoid
+// code duplication. However, if we add more mesh types we should refactor.
+template <bool MarchingCubesConnectivity>
 class MeshBoundary3D
 {
 public:
@@ -106,22 +110,69 @@ public:
   }
 
   VTKM_EXEC_CONT
-  vtkm::Id CountLinkComponentsIn2DSlice(const vtkm::Id meshIndex, const vtkm::Id2 strides) const
+  vtkm::Id CountLinkComponentsIn2DSlice4Neighborhood(const vtkm::Id meshIndex,
+                                                     const vtkm::Id2 strides) const
   {
-    // IMPORTANT: We assume that function is called only for *interior* vertices (i.e., neither row nor col
+    // IMPORTANT: We assume that function is called only for *interior* vertices (i.e.,
-    // within slice is 0 and we do not need to check for boundary cases).
+    // neither row nor col within slice is 0 and we do not need to check for boundary cases).
     vtkm::Id sortIndex = this->SortIndicesPortal.Get(meshIndex);
     bool prevWasInUpperLink = false;
-    vtkm::Id numComponents = 0;
+    vtkm::Id numLinkTypeChanges = 0;
 
-    const int N_INCIDENT_EDGES_2D = 6;
+    for (vtkm::Id edgeNo = 0; edgeNo < 4 + 1; edgeNo++)
-    for (vtkm::Id edgeNo = 0; edgeNo < N_INCIDENT_EDGES_2D; edgeNo++)
+    {
-    { // per edge
       VTKM_ASSERT(meshIndex + strides[1] + strides[0] <
                   this->SortIndicesPortal.GetNumberOfValues());
       VTKM_ASSERT(meshIndex - strides[1] - strides[0] >= 0);
       vtkm::Id nbrSortIndex;
-      switch (edgeNo)
+      switch (edgeNo % 4)
+      {
+        case 0:
+          nbrSortIndex = this->SortIndicesPortal.Get(meshIndex + strides[0]);
+          break; // [1]    , [0] + 1
+        case 1:
+          nbrSortIndex = this->SortIndicesPortal.Get(meshIndex + strides[1]);
+          break; // [1] + 1, [0]
+        case 2:
+          nbrSortIndex = this->SortIndicesPortal.Get(meshIndex - strides[0]);
+          break; // [1]    , [0] - 1
+        case 3:
+          nbrSortIndex = this->SortIndicesPortal.Get(meshIndex - strides[1]);
+          break; // [1] - 1, [0]
+        default:
+          // Should not occur, since switch statement is over edgeNo % 4
+          VTKM_ASSERT(false);
+          // Initialize nbrSortIndex to something anyway to prevent compiler warning
+          // Set to the sort index of the vertex itself since there is "no" edge so
+          // that it contains a "sane" value if it should ever be reached.
+          nbrSortIndex = this->SortIndicesPortal.Get(meshIndex);
+          break;
+      } // switch edgeNo
+
+      bool currIsInUpperLink = (nbrSortIndex > sortIndex);
+      numLinkTypeChanges += (edgeNo != 0 && currIsInUpperLink != prevWasInUpperLink) ? 1 : 0;
+      prevWasInUpperLink = currIsInUpperLink;
+    } // per edge
+    return numLinkTypeChanges == 0 ? 1 : numLinkTypeChanges;
+  }
+
+  VTKM_EXEC_CONT
+  vtkm::Id CountLinkComponentsIn2DSlice6Neighborhood(const vtkm::Id meshIndex,
+                                                     const vtkm::Id2 strides) const
+  {
+    // IMPORTANT: We assume that function is called only for *interior* vertices (i.e.,
+    // neither row nor col within slice is 0 and we do not need to check for boundary cases).
+    vtkm::Id sortIndex = this->SortIndicesPortal.Get(meshIndex);
+    bool prevWasInUpperLink = false;
+    vtkm::Id numLinkTypeChanges = 0;
+
+    for (vtkm::Id edgeNo = 0; edgeNo < 6 + 1; edgeNo++)
+    {
+      VTKM_ASSERT(meshIndex + strides[1] + strides[0] <
+                  this->SortIndicesPortal.GetNumberOfValues());
+      VTKM_ASSERT(meshIndex - strides[1] - strides[0] >= 0);
+      vtkm::Id nbrSortIndex;
+      switch (edgeNo % 6)
       {
         case 0:
           nbrSortIndex = this->SortIndicesPortal.Get(meshIndex + strides[0]);
@@ -142,20 +193,80 @@ public:
           nbrSortIndex = this->SortIndicesPortal.Get(meshIndex - strides[1]);
           break; // [1] - 1, [0]
         default:
-          // Due to CUDA we cannot throw an exception here, which would make the most
+          // Should not occur, since switch statement is over edgeNo % 6
-          // sense
+          VTKM_ASSERT(false);
-          VTKM_ASSERT(false); // Should not occur, edgeNo < N_INCIDENT_EDGES_2D = 6
           // Initialize nbrSortIndex to something anyway to prevent compiler warning
           // Set to the sort index of the vertex itself since there is "no" edge so
           // that it contains a "sane" value if it should ever be reached.
           nbrSortIndex = this->SortIndicesPortal.Get(meshIndex);
           break;
-      }
+      } // seitch edgeNo
 
       bool currIsInUpperLink = (nbrSortIndex > sortIndex);
-      numComponents += (edgeNo != 0 && currIsInUpperLink != prevWasInUpperLink) ? 1 : 0;
+      numLinkTypeChanges += (edgeNo != 0 && currIsInUpperLink != prevWasInUpperLink) ? 1 : 0;
+      prevWasInUpperLink = currIsInUpperLink;
     } // per edge
-    return numComponents;
+    return numLinkTypeChanges == 0 ? 1 : numLinkTypeChanges;
+  }
+
+  VTKM_EXEC_CONT
+  vtkm::Id CountLinkComponentsIn2DSlice8Neighborhood(const vtkm::Id meshIndex,
+                                                     const vtkm::Id2 strides) const
+  {
+    // IMPORTANT: We assume that function is called only for *interior* vertices (i.e.,
+    // neither row nor col within slice is 0 and we do not need to check for boundary cases).
+    vtkm::Id sortIndex = this->SortIndicesPortal.Get(meshIndex);
+    bool prevWasInUpperLink = false;
+    vtkm::Id numLinkTypeChanges = 0;
+
+    for (vtkm::Id edgeNo = 0; edgeNo < 8 + 1; edgeNo++)
+    {
+      VTKM_ASSERT(meshIndex + strides[1] + strides[0] <
+                  this->SortIndicesPortal.GetNumberOfValues());
+      VTKM_ASSERT(meshIndex - strides[1] - strides[0] >= 0);
+      vtkm::Id nbrSortIndex;
+
+      switch (edgeNo % 8)
+      {
+        case 0:
+          nbrSortIndex = this->SortIndicesPortal.Get(meshIndex + strides[0]);
+          break; // [1]    , [0] + 1
+        case 1:
+          nbrSortIndex = this->SortIndicesPortal.Get(meshIndex + strides[1] + strides[0]);
+          break; // [1] + 1, [0] + 1
+        case 2:
+          nbrSortIndex = this->SortIndicesPortal.Get(meshIndex + strides[1]);
+          break; // [1] + 1, [0]
+        case 3:
+          nbrSortIndex = this->SortIndicesPortal.Get(meshIndex + strides[1] - strides[0]);
+          break; // [1] + 1, [0] -1
+        case 4:
+          nbrSortIndex = this->SortIndicesPortal.Get(meshIndex - strides[0]);
+          break; // [1]    , [0] - 1
+        case 5:
+          nbrSortIndex = this->SortIndicesPortal.Get(meshIndex - strides[1] - strides[0]);
+          break; // [1] - 1, [0] - 1
+        case 6:
+          nbrSortIndex = this->SortIndicesPortal.Get(meshIndex - strides[1]);
+          break; // [1] - 1, [0]
+        case 7:
+          nbrSortIndex = this->SortIndicesPortal.Get(meshIndex - strides[1] + strides[0]);
+          break; // [1] - 1, [0] + 1
+        default:
+          // Should not occur, since switch statement is over edgeNo % 8
+          VTKM_ASSERT(false);
+          // Initialize nbrSortIndex to something anyway to prevent compiler warning
+          // Set to the sort index of the vertex itself since there is "no" edge so
+          // that it contains a "sane" value if it should ever be reached.
+          nbrSortIndex = this->SortIndicesPortal.Get(meshIndex);
+          break;
+      } // switch edgeNo
+
+      bool currIsInUpperLink = (nbrSortIndex > sortIndex);
+      numLinkTypeChanges += (edgeNo != 0 && currIsInUpperLink != prevWasInUpperLink) ? 1 : 0;
+      prevWasInUpperLink = currIsInUpperLink;
+    } // per edge
+    return numLinkTypeChanges == 0 ? 1 : numLinkTypeChanges;
   }
 
   VTKM_EXEC_CONT
@@ -192,11 +303,20 @@ public:
             (pos[1] == this->MeshStructure.MeshSize[1] - 1 &&
              pos[2] == this->MeshStructure.MeshSize[2] - 1))
         {
-          VTKM_ASSERT(meshIndex >= 1);
+          if (MarchingCubesConnectivity)
-          vtkm::Id sp = this->SortIndicesPortal.Get(meshIndex - 1);
+          {
-          VTKM_ASSERT(meshIndex + 1 < this->SortIndicesPortal.GetNumberOfValues());
+            // TODO/FIXME: Be more selective about what points to include, but will likely
-          vtkm::Id sn = this->SortIndicesPortal.Get(meshIndex + 1);
+            // require an additional layer of "ghost cells"
-          return (sortIndex < sp && sortIndex < sn) || (sortIndex > sp && sortIndex > sn);
+            return true;
+          }
+          else
+          {
+            VTKM_ASSERT(meshIndex >= 1);
+            vtkm::Id sp = this->SortIndicesPortal.Get(meshIndex - 1);
+            VTKM_ASSERT(meshIndex + 1 < this->SortIndicesPortal.GetNumberOfValues());
+            vtkm::Id sn = this->SortIndicesPortal.Get(meshIndex + 1);
+            return (sortIndex < sp && sortIndex < sn) || (sortIndex > sp && sortIndex > sn);
+          }
         }
         // Edges in [1] directtion
         else if ((pos[0] == 0 && pos[2] == 0) ||
@@ -205,13 +325,22 @@ public:
                  (pos[0] == this->MeshStructure.MeshSize[0] - 1 &&
                   pos[2] == this->MeshStructure.MeshSize[2] - 1))
         {
-          VTKM_ASSERT(pos[1] > 0 && pos[1] < this->MeshStructure.MeshSize[1] - 1);
+          if (MarchingCubesConnectivity)
-          VTKM_ASSERT(meshIndex >= this->MeshStructure.MeshSize[0]);
+          {
-          vtkm::Id sp = this->SortIndicesPortal.Get(meshIndex - this->MeshStructure.MeshSize[0]);
+            // TODO/FIXME: Be more selective about what points to include, but will likely
-          VTKM_ASSERT(meshIndex + this->MeshStructure.MeshSize[0] <
+            // require an additional layer of "ghost cells"
-                      this->SortIndicesPortal.GetNumberOfValues());
+            return true;
-          vtkm::Id sn = this->SortIndicesPortal.Get(meshIndex + this->MeshStructure.MeshSize[0]);
+          }
-          return (sortIndex < sp && sortIndex < sn) || (sortIndex > sp && sortIndex > sn);
+          else
+          {
+            VTKM_ASSERT(pos[1] > 0 && pos[1] < this->MeshStructure.MeshSize[1] - 1);
+            VTKM_ASSERT(meshIndex >= this->MeshStructure.MeshSize[0]);
+            vtkm::Id sp = this->SortIndicesPortal.Get(meshIndex - this->MeshStructure.MeshSize[0]);
+            VTKM_ASSERT(meshIndex + this->MeshStructure.MeshSize[0] <
+                        this->SortIndicesPortal.GetNumberOfValues());
+            vtkm::Id sn = this->SortIndicesPortal.Get(meshIndex + this->MeshStructure.MeshSize[0]);
+            return (sortIndex < sp && sortIndex < sn) || (sortIndex > sp && sortIndex > sn);
+          }
         }
         // Edges in [2] direction
         else if ((pos[1] == 0 && pos[0] == 0) ||
@@ -220,11 +349,20 @@ public:
                  (pos[1] == this->MeshStructure.MeshSize[1] - 1 &&
                   pos[0] == this->MeshStructure.MeshSize[0] - 1))
         {
-          VTKM_ASSERT(meshIndex >= nPerSlice);
+          if (MarchingCubesConnectivity)
-          vtkm::Id sp = this->SortIndicesPortal.Get(meshIndex - nPerSlice);
+          {
-          VTKM_ASSERT(meshIndex + nPerSlice < this->SortIndicesPortal.GetNumberOfValues());
+            // TODO/FIXME: Be more selective about what points to include, but will likely
-          vtkm::Id sn = this->SortIndicesPortal.Get(meshIndex + nPerSlice);
+            // require an additional layer of "ghost cells"
-          return (sortIndex < sp && sortIndex < sn) || (sortIndex > sp && sortIndex > sn);
+            return true;
+          }
+          else
+          {
+            VTKM_ASSERT(meshIndex >= nPerSlice);
+            vtkm::Id sp = this->SortIndicesPortal.Get(meshIndex - nPerSlice);
+            VTKM_ASSERT(meshIndex + nPerSlice < this->SortIndicesPortal.GetNumberOfValues());
+            vtkm::Id sn = this->SortIndicesPortal.Get(meshIndex + nPerSlice);
+            return (sortIndex < sp && sortIndex < sn) || (sortIndex > sp && sortIndex > sn);
+          }
         }
         else
         {
@@ -233,22 +371,52 @@ public:
           { // On [2]-perpendicular face
             VTKM_ASSERT(pos[0] != 0 && pos[0] != this->MeshStructure.MeshSize[0]);
             VTKM_ASSERT(pos[1] != 0 && pos[1] != this->MeshStructure.MeshSize[1]);
-            return CountLinkComponentsIn2DSlice(meshIndex,
+            if (MarchingCubesConnectivity)
-                                                vtkm::Id2(1, this->MeshStructure.MeshSize[0])) != 2;
+            {
+              return CountLinkComponentsIn2DSlice4Neighborhood(
+                       meshIndex, vtkm::Id2(1, this->MeshStructure.MeshSize[0])) != 2 ||
+                CountLinkComponentsIn2DSlice8Neighborhood(
+                  meshIndex, vtkm::Id2(1, this->MeshStructure.MeshSize[0])) != 2;
+            }
+            else
+            {
+              return CountLinkComponentsIn2DSlice6Neighborhood(
+                       meshIndex, vtkm::Id2(1, this->MeshStructure.MeshSize[0])) != 2;
+            }
           }
           else if (pos[1] == 0 || pos[1] == this->MeshStructure.MeshSize[1] - 1)
           { // On [1]-perpendicular face
             VTKM_ASSERT(pos[0] != 0 && pos[0] != this->MeshStructure.MeshSize[0]);
             VTKM_ASSERT(pos[2] != 0 && pos[2] != this->MeshStructure.MeshSize[2]);
-            return CountLinkComponentsIn2DSlice(meshIndex, vtkm::Id2(1, nPerSlice)) != 2;
+            if (MarchingCubesConnectivity)
+            {
+              return CountLinkComponentsIn2DSlice4Neighborhood(meshIndex,
+                                                               vtkm::Id2(1, nPerSlice)) != 2 ||
+                CountLinkComponentsIn2DSlice8Neighborhood(meshIndex, vtkm::Id2(1, nPerSlice)) != 2;
+            }
+            else
+            {
+              return CountLinkComponentsIn2DSlice6Neighborhood(meshIndex,
+                                                               vtkm::Id2(1, nPerSlice)) != 2;
+            }
           }
           else
           { // On [0]-perpendicular face
             VTKM_ASSERT(pos[0] == 0 || pos[0] == this->MeshStructure.MeshSize[0] - 1);
             VTKM_ASSERT(pos[1] != 0 && pos[1] != this->MeshStructure.MeshSize[1]);
             VTKM_ASSERT(pos[2] != 0 && pos[2] != this->MeshStructure.MeshSize[2]);
-            return CountLinkComponentsIn2DSlice(
+            if (MarchingCubesConnectivity)
-                     meshIndex, vtkm::Id2(nPerSlice, this->MeshStructure.MeshSize[0])) != 2;
+            {
+              return CountLinkComponentsIn2DSlice4Neighborhood(
+                       meshIndex, vtkm::Id2(nPerSlice, this->MeshStructure.MeshSize[0])) != 2 ||
+                CountLinkComponentsIn2DSlice8Neighborhood(
+                  meshIndex, vtkm::Id2(nPerSlice, this->MeshStructure.MeshSize[0])) != 2;
+            }
+            else
+            {
+              return CountLinkComponentsIn2DSlice6Neighborhood(
+                       meshIndex, vtkm::Id2(nPerSlice, this->MeshStructure.MeshSize[0])) != 2;
+            }
           }
         }
       }
@@ -269,6 +437,10 @@ protected:
 };
 
 
+// TODO/FIXME: Consider making marching cubes connectivity its own class. However
+// at the moment making this a boolean template parameter makes it easuer to avoid
+// code duplication. However, if we add more mesh types we should refactor.
+template <bool MarchingCubesConnectivity>
 class MeshBoundary3DExec : public vtkm::cont::ExecutionObjectBase
 {
 public:
@@ -279,10 +451,12 @@ public:
   {
   }
 
-  VTKM_CONT MeshBoundary3D PrepareForExecution(vtkm::cont::DeviceAdapterId device,
+  VTKM_CONT MeshBoundary3D<MarchingCubesConnectivity> PrepareForExecution(
-                                               vtkm::cont::Token& token) const
+    vtkm::cont::DeviceAdapterId device,
+    vtkm::cont::Token& token) const
   {
-    return MeshBoundary3D(this->MeshSize, this->SortIndices, device, token);
+    return MeshBoundary3D<MarchingCubesConnectivity>(
+      this->MeshSize, this->SortIndices, device, token);
   }
 
 protected: