(Bug)fix to suppress off-block regular nodes in Augmented Tree

vtkm/filter/scalar_topology/ContourTreeUniformDistributed.cxx

@@ -1138,11 +1138,17 @@ VTKM_CONT void ContourTreeUniformDistributed::DoPostExecute(
   // ******** 3. Augment the hierarchical tree if requested ********
   if (this->AugmentHierarchicalTree)
   {
-    master.foreach (
-      [](DistributedContourTreeBlockData* blockData, const vtkmdiy::Master::ProxyWithLink&) {
-        blockData->HierarchicalAugmenter.Initialize(
-          blockData->GlobalBlockId, &blockData->HierarchicalTree, &blockData->AugmentedTree);
-      });
+    master.foreach ([globalPointDimensions](DistributedContourTreeBlockData* blockData,
+                                            const vtkmdiy::Master::ProxyWithLink&) {
+      blockData->HierarchicalAugmenter.Initialize(
+        blockData->GlobalBlockId,
+        &blockData->HierarchicalTree,
+        &blockData->AugmentedTree,
+        blockData->BlockOrigin, // Origin of the data block
+        blockData->BlockSize,   // Extends of the data block
+        globalPointDimensions   // global point dimensions
+      );
+    });
 
     timingsStream << "    " << std::setw(38) << std::left << "Initalize Hierarchical Trees"
                   << ": " << timer.GetElapsedTime() << " seconds" << std::endl;

vtkm/filter/scalar_topology/worklet/contourtree_distributed/HierarchicalAugmenter.h

@@ -127,6 +127,11 @@ template <typename FieldType>
 class HierarchicalAugmenter
 { // class HierarchicalAugmenter
 public:
+  /// base mesh variable needs to determine whether a vertex is inside or outside of the block
+  vtkm::Id3 MeshBlockOrigin;
+  vtkm::Id3 MeshBlockSize;
+  vtkm::Id3 MeshGlobalSize;
+
   /// the tree that it hypersweeps over
   vtkm::worklet::contourtree_distributed::HierarchicalContourTree<FieldType>* BaseTree;
   /// the tree that it is building
@@ -198,7 +203,10 @@ public:
   void Initialize(
     vtkm::Id blockId,
     vtkm::worklet::contourtree_distributed::HierarchicalContourTree<FieldType>* inBaseTree,
-    vtkm::worklet::contourtree_distributed::HierarchicalContourTree<FieldType>* inAugmentedTree);
+    vtkm::worklet::contourtree_distributed::HierarchicalContourTree<FieldType>* inAugmentedTree,
+    vtkm::Id3 meshBlockOrigin,
+    vtkm::Id3 meshBockSize,
+    vtkm::Id3 meshGlobalSize);
 
   /// routine to prepare the set of attachment points to transfer
   void PrepareOutAttachmentPoints(vtkm::Id round);
@@ -249,12 +257,18 @@ template <typename FieldType>
 void HierarchicalAugmenter<FieldType>::Initialize(
   vtkm::Id blockId,
   vtkm::worklet::contourtree_distributed::HierarchicalContourTree<FieldType>* baseTree,
-  vtkm::worklet::contourtree_distributed::HierarchicalContourTree<FieldType>* augmentedTree)
+  vtkm::worklet::contourtree_distributed::HierarchicalContourTree<FieldType>* augmentedTree,
+  vtkm::Id3 meshBlockOrigin,
+  vtkm::Id3 meshBockSize,
+  vtkm::Id3 meshGlobalSize)
 { // Initialize()
   // copy the parameters for use
   this->BlockId = blockId;
   this->BaseTree = baseTree;
   this->AugmentedTree = augmentedTree;
+  this->MeshBlockOrigin = meshBlockOrigin;
+  this->MeshBlockSize = meshBockSize;
+  this->MeshGlobalSize = meshGlobalSize;
 
   // now construct a list of all attachment points on the block
   // to do this, we construct an index array with all supernode ID's that satisfy:
@@ -728,12 +742,13 @@ void HierarchicalAugmenter<FieldType>::CopyBaseRegularStructure()
     vtkm::worklet::contourtree_augmented::IdArrayType tempRegularNodesNeeded;
     // create the worklet
     vtkm::worklet::contourtree_distributed::hierarchical_augmenter::
-      FindSuperparentForNecessaryNodesWorklet findSuperparentForNecessaryNodesWorklet;
+      FindSuperparentForNecessaryNodesWorklet findSuperparentForNecessaryNodesWorklet(
+        this->MeshBlockOrigin, this->MeshBlockSize, this->MeshGlobalSize);
     // Get a FindRegularByGlobal and FindSuperArcForUnknownNode execution object for our worklet
     auto findRegularByGlobal = this->AugmentedTree->GetFindRegularByGlobal();
     auto findSuperArcForUnknownNode = this->AugmentedTree->GetFindSuperArcForUnknownNode();
 
-    // excute the worklet
+    // execute the worklet
     this->Invoke(findSuperparentForNecessaryNodesWorklet, // the worklet to call
                  // inputs
                  this->BaseTree->RegularNodeGlobalIds, // input domain
@@ -831,6 +846,12 @@ void HierarchicalAugmenter<FieldType>::CopyBaseRegularStructure()
     );
   }
 
+
+  // Reset the number of regular nodes in round 0
+  vtkm::Id regularNodesInRound0 =
+    numTotalRegular - this->AugmentedTree->NumRegularNodesInRound.ReadPortal().Get(1);
+  this->AugmentedTree->NumRegularNodesInRound.WritePortal().Set(0, regularNodesInRound0);
+
   //  Finally, we resort the regular node sort order
   {
     vtkm::worklet::contourtree_distributed::PermuteComparator // hierarchical_contour_tree::

vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_augmenter/FindSuperparentForNecessaryNodesWorklet.h

@@ -81,9 +81,17 @@ public:
   using ExecutionSignature = void(InputIndex, _1, _2, _3, _4, _5, _6, _7, _8, _9);
   using InputDomain = _1;
 
+
   /// Default Constructor
   VTKM_EXEC_CONT
-  FindSuperparentForNecessaryNodesWorklet() {}
+  FindSuperparentForNecessaryNodesWorklet(vtkm::Id3 meshBlockOrigin,
+                                          vtkm::Id3 meshBlockSize,
+                                          vtkm::Id3 meshGlobalSize)
+    : MeshBlockOrigin(meshBlockOrigin)
+    , MeshBlockSize(meshBlockSize)
+    , MeshGlobalSize(meshGlobalSize)
+  {
+  }
 
   /// operator() of the workelt
   template <typename InFieldPortalType,
@@ -111,8 +119,16 @@ public:
     // first check to see if it is already present (newRegularId set on input)
     vtkm::Id newRegularId = findRegularByGlobal.FindRegularByGlobal(globalRegularId);
 
+    // Explicitly check whether the vertex belongs to the base block. If it doesn't, we ignore it
+    if (!this->IsInMesh(globalRegularId))
+    {
+      // Set to NO_SUCH_ELEMENT by default. By doing this in the worklet we an avoid having to
+      // initialize the output arrays first and we can use FieldIn instead of FieldInOut
+      regularSuperparentsValue = vtkm::worklet::contourtree_augmented::NO_SUCH_ELEMENT;
+      regularNodesNeededValue = vtkm::worklet::contourtree_augmented::NO_SUCH_ELEMENT;
+    }
     // if it fails this test, then it's already in tree
-    if (vtkm::worklet::contourtree_augmented::NoSuchElement(newRegularId))
+    else if (vtkm::worklet::contourtree_augmented::NoSuchElement(newRegularId))
     { // not yet in tree
       // since it's not in the tree, we want to find where it belongs
       // to do so, we need to find an "above" and "below" node for it. Since it exists in the old tree, it belongs to a superarc, and we
@@ -200,7 +216,80 @@ public:
     */
   } // operator()()
 
-}; // FindSuperparentForNecessaryNodesWorklet
+
+private:
+  // Mesh data
+  vtkm::Id3 MeshBlockOrigin;
+  vtkm::Id3 MeshBlockSize;
+  vtkm::Id3 MeshGlobalSize;
+
+  VTKM_EXEC
+  bool IsInMesh(vtkm::Id globalId) const
+  {                                  // IsInMesh()
+    if (this->MeshGlobalSize[2] > 1) // 3D
+    {
+      // convert from global ID to global coords
+      vtkm::Id globalSliceSize = this->MeshGlobalSize[0] * this->MeshGlobalSize[1];
+      vtkm::Id globalSlice = globalId / globalSliceSize;
+      vtkm::Id globalRow = globalId / this->MeshGlobalSize[0];
+      vtkm::Id globalCol = globalId % this->MeshGlobalSize[0];
+
+      // test validity
+      if (globalSlice < this->MeshBlockOrigin[2])
+      {
+        return false;
+      }
+      if (globalSlice >= this->MeshBlockOrigin[2] + this->MeshBlockSize[2])
+      {
+        return false;
+      }
+      if (globalRow < this->MeshBlockOrigin[1])
+      {
+        return false;
+      }
+      if (globalRow >= this->MeshBlockOrigin[1] + this->MeshBlockSize[1])
+      {
+        return false;
+      }
+      if (globalCol < this->MeshBlockOrigin[0])
+      {
+        return false;
+      }
+      if (globalCol >= this->MeshBlockOrigin[0] + this->MeshBlockSize[0])
+      {
+        return false;
+      }
+      // it's in the block - return true
+      return true;
+    }    // end if 3D
+    else // 2D mesh
+    {
+      // convert from global ID to global coords
+      vtkm::Id globalRow = globalId / this->MeshGlobalSize[0];
+      vtkm::Id globalCol = globalId % this->MeshGlobalSize[0];
+
+      // test validity
+      if (globalRow < this->MeshBlockOrigin[1])
+      {
+        return false;
+      }
+      if (globalRow >= this->MeshBlockOrigin[1] + this->MeshBlockSize[1])
+      {
+        return false;
+      }
+      if (globalCol < this->MeshBlockOrigin[0])
+      {
+        return false;
+      }
+      if (globalCol >= this->MeshBlockOrigin[0] + this->MeshBlockSize[0])
+      {
+        return false;
+      }
+      // it's in the block - return true
+      return true;
+    }
+  } // IsInMesh()
+};  // FindSuperparentForNecessaryNodesWorklet
 
 } // namespace hierarchical_augmenter
 } // namespace contourtree_distributed