Merge branch 'add/presimplification' into 'master'

Draft: Add pre-simplification to address scaling performance issues

See merge request vtk/vtk-m!3227

docs/changelog/env-options.md

@@ -0,0 +1,9 @@
+## Load options from environment variables
+
+Some common VTK-m options such as the device and log level could be
+specified on the command line but not through environment variables. It is
+not always possible to set VTK-m command line options, so environment
+variables are added.
+
+Also added documentation to the user's guide about what options are
+available and how to set them.

docs/users-guide/initialization.rst

@@ -19,6 +19,44 @@ But it can also optionally take the ``argc`` and ``argv`` arguments to the ``mai
 |VTKm| accepts arguments that, for example, configure the compute device to use or establish logging levels.
 Any arguments that are handled by |VTKm| are removed from the ``argc``/``argv`` list so that your program can then respond to the remaining arguments.
 
+Many options can also be set with environment variables.
+If both the environment variable and command line argument are provided, the command line argument is used.
+The following table lists the currently supported options.
+
+.. list-table:: |VTKm| command line arguments and environment variable options.
+   :widths: 23 22 15 40
+   :header-rows: 1
+
+   * - Command Line Argument
+     - Environment Variable
+     - Default Value
+     - Description
+   * - ``--vtkm-help``
+     -
+     -
+     - Causes the program to print information about |VTKm| command line arguments and then exits the program.
+   * - ``--vtkm-log-level``
+     - ``VTKM_LOG_LEVEL``
+     - ``WARNING``
+     - Specifies the logging level.
+       Valid values are ``INFO``, ``WARNING``, ``ERROR``, ``FATAL``, and ``OFF``.
+       This can also be set to a numeric value for the logging level.
+   * - ``--vtkm-device``
+     - ``VTKM_DEVICE``
+     -
+     - Force |VTKm| to use the specified device.
+       If not specified or ``Any`` given, then any available device may be used.
+   * - ``--vtkm-num-threads``
+     - ``VTKM_NUM_THREADS``
+     -
+     - Set the number of threads to use on a multi-core device.
+       If not specified, the device will use the cores available in the system.
+   * - ``--vtkm-device-instance``
+     - ``VTKM_DEVICE_INSTANCE``
+     -
+     - Selects the device to use when more than one device device of a given type is available.
+       The device is specified with a numbered index.
+
 :func:`vtkm::cont::Initialize` returns a :struct:`vtkm::cont::InitializeResult` structure.
 This structure contains information about the supported arguments and options selected during initialization.
 

examples/contour_tree_distributed/ContourTreeApp.cxx

@@ -220,6 +220,12 @@ int main(int argc, char* argv[])
     }
   }
 
+  vtkm::Id presimplifyThreshold = 0; // Do not presimplify the hierachical contour tree
+  if (parser.hasOption("--presimplifyThreshold"))
+  {
+    presimplifyThreshold = std::stoi(parser.getOption("--presimplifyThreshold"));
+  }
+
   bool useBoundaryExtremaOnly = true;
   if (parser.hasOption("--useFullBoundary"))
   {
@@ -368,6 +374,10 @@ int main(int argc, char* argv[])
       std::cout
         << "--eps=<float>   Floating point offset awary from the critical point. (default=0.00001)"
         << std::endl;
+      std::cout << "--presimplifyThreshold   Integer volume threshold for presimplifying the tree"
+                << std::endl;
+      std::cout << "                 Default value is 0, indicating no presimplification"
+                << std::endl;
       std::cout << "--preSplitFiles  Input data is already pre-split into blocks." << std::endl;
       std::cout << "--saveDot        Save DOT files of the distributed contour tree " << std::endl
                 << "                 computation (Default=False). " << std::endl;
@@ -411,6 +421,7 @@ int main(int argc, char* argv[])
                  << "    augmentHierarchicalTree=" << augmentHierarchicalTree << std::endl
                  << "    computeVolumetricBranchDecomposition="
                  << computeHierarchicalVolumetricBranchDecomposition << std::endl
+                 << "    presimplifyThreshold=" << presimplifyThreshold << std::endl
                  << "    saveOutputData=" << saveOutputData << std::endl
                  << "    forwardSummary=" << forwardSummary << std::endl
                  << "    nblocks=" << numBlocks << std::endl
@@ -648,6 +659,10 @@ int main(int argc, char* argv[])
   filter.SetUseBoundaryExtremaOnly(useBoundaryExtremaOnly);
   filter.SetUseMarchingCubes(useMarchingCubes);
   filter.SetAugmentHierarchicalTree(augmentHierarchicalTree);
+  if (presimplifyThreshold > 0)
+  {
+    filter.SetPresimplifyThreshold(presimplifyThreshold);
+  }
   filter.SetSaveDotFiles(saveDotFiles);
   filter.SetActiveField("values");
 
@@ -923,6 +938,7 @@ int main(int argc, char* argv[])
     close(save_err);
   }
 
+  std::cout << "DONE!!!" << std::endl;
   // Finalize and finish the execution
   MPI_Finalize();
   return EXIT_SUCCESS;

vtkm/cont/Initialize.cxx

@@ -17,6 +17,7 @@
 
 #include <vtkm/thirdparty/diy/environment.h>
 
+#include <cstdlib>
 #include <memory>
 #include <sstream>
 
@@ -123,7 +124,7 @@ InitializeResult Initialize(int& argc, char* argv[], InitializeOptions opts)
   }
   else
   {
-    vtkm::cont::InitLogging(argc, argv, loggingFlag);
+    vtkm::cont::InitLogging(argc, argv, loggingFlag, "VTKM_LOG_LEVEL");
   }
   if (!vtkmdiy::mpi::environment::initialized())
   {
@@ -225,37 +226,70 @@ InitializeResult Initialize(int& argc, char* argv[], InitializeOptions opts)
         vtkm::cont::DeviceAdapterTagAny{}, runtimeDeviceOptions, argc, argv);
     }
 
+    // Check for device on command line.
     if (options[opt::OptionIndex::DEVICE])
     {
       const char* arg = options[opt::OptionIndex::DEVICE].arg;
-      auto id = vtkm::cont::make_DeviceAdapterId(arg);
-      if (id != vtkm::cont::DeviceAdapterTagAny{})
+      config.Device = vtkm::cont::make_DeviceAdapterId(arg);
+    }
+    // If not on command line, check for device in environment variable.
+    if (config.Device == vtkm::cont::DeviceAdapterTagUndefined{})
+    {
+      const char* deviceEnv = std::getenv("VTKM_DEVICE");
+      if (deviceEnv != nullptr)
       {
-        vtkm::cont::GetRuntimeDeviceTracker().ForceDevice(id);
+        auto id = vtkm::cont::make_DeviceAdapterId(std::getenv("VTKM_DEVICE"));
+        if (VtkmDeviceArg::DeviceIsAvailable(id))
+        {
+          config.Device = id;
+        }
+        else
+        {
+          // Got invalid device. Log an error, but continue to do the default action for
+          // the device (i.e., ignore the environment variable setting).
+          VTKM_LOG_S(vtkm::cont::LogLevel::Error,
+                     "Invalid device `"
+                       << deviceEnv
+                       << "` specified in VTKM_DEVICE environment variable. Ignoring.");
+          VTKM_LOG_S(vtkm::cont::LogLevel::Error,
+                     "Valid devices are: " << VtkmDeviceArg::GetValidDeviceNames());
+        }
+      }
+    }
+    // If still not defined, check to see if "any" device should be added.
+    if ((config.Device == vtkm::cont::DeviceAdapterTagUndefined{}) &&
+        (opts & InitializeOptions::DefaultAnyDevice) != InitializeOptions::None)
+    {
+      config.Device = vtkm::cont::DeviceAdapterTagAny{};
+    }
+    // Set the state for the device selected.
+    if (config.Device == vtkm::cont::DeviceAdapterTagUndefined{})
+    {
+      if ((opts & InitializeOptions::RequireDevice) != InitializeOptions::None)
+      {
+        auto devices = VtkmDeviceArg::GetValidDeviceNames();
+        VTKM_LOG_S(vtkm::cont::LogLevel::Fatal, "Device not given on command line.");
+        std::cerr << "Target device must be specified via --vtkm-device.\n"
+                     "Valid devices: "
+                  << devices << std::endl;
+        if ((opts & InitializeOptions::AddHelp) != InitializeOptions::None)
+        {
+          std::cerr << config.Usage;
+        }
+        exit(1);
       }
       else
       {
-        vtkm::cont::GetRuntimeDeviceTracker().Reset();
+        // No device specified. Do nothing and let VTK-m decide what it is going to do.
       }
-      config.Device = id;
     }
-    else if ((opts & InitializeOptions::DefaultAnyDevice) != InitializeOptions::None)
+    else if (config.Device == vtkm::cont::DeviceAdapterTagAny{})
     {
       vtkm::cont::GetRuntimeDeviceTracker().Reset();
-      config.Device = vtkm::cont::DeviceAdapterTagAny{};
     }
-    else if ((opts & InitializeOptions::RequireDevice) != InitializeOptions::None)
+    else
     {
-      auto devices = VtkmDeviceArg::GetValidDeviceNames();
-      VTKM_LOG_S(vtkm::cont::LogLevel::Error, "Device not given on command line.");
-      std::cerr << "Target device must be specified via --vtkm-device.\n"
-                   "Valid devices: "
-                << devices << std::endl;
-      if ((opts & InitializeOptions::AddHelp) != InitializeOptions::None)
-      {
-        std::cerr << config.Usage;
-      }
-      exit(1);
+      vtkm::cont::GetRuntimeDeviceTracker().ForceDevice(config.Device);
     }
 
 

vtkm/cont/Logging.cxx

@@ -30,7 +30,7 @@
 
 #endif // VTKM_ENABLE_LOGGING
 
-#include <cassert>
+#include <cstdlib>
 #include <iomanip>
 #include <sstream>
 #include <stdexcept>
@@ -108,7 +108,10 @@ namespace cont
 {
 
 VTKM_CONT
-void InitLogging(int& argc, char* argv[], const std::string& loggingFlag)
+void InitLogging(int& argc,
+                 char* argv[],
+                 const std::string& loggingFlag,
+                 const std::string& loggingEnv)
 {
   SetLogLevelName(vtkm::cont::LogLevel::Off, "Off");
   SetLogLevelName(vtkm::cont::LogLevel::Fatal, "FATL");
@@ -130,8 +133,16 @@ void InitLogging(int& argc, char* argv[], const std::string& loggingFlag)
     loguru::set_verbosity_to_name_callback(&verbosityToNameCallback);
     loguru::set_name_to_verbosity_callback(&nameToVerbosityCallback);
 
-    // Set the default log level to warning
-    SetStderrLogLevel(vtkm::cont::LogLevel::Warn);
+    const char* envLevel = std::getenv(loggingEnv.c_str());
+    if (envLevel != nullptr)
+    {
+      SetStderrLogLevel(envLevel);
+    }
+    else
+    {
+      // Set the default log level to warning
+      SetStderrLogLevel(vtkm::cont::LogLevel::Warn);
+    }
     loguru::init(argc, argv, loggingFlag.c_str());
   }
 #else  // VTKM_ENABLE_LOGGING

vtkm/cont/Logging.h

@@ -371,7 +371,10 @@ enum class LogLevel
  */
 VTKM_CONT_EXPORT
 VTKM_CONT
-void InitLogging(int& argc, char* argv[], const std::string& loggingFlag = "--vtkm-log-level");
+void InitLogging(int& argc,
+                 char* argv[],
+                 const std::string& loggingFlag = "--vtkm-log-level",
+                 const std::string& loggingEnv = "VTKM_LOG_LEVEL");
 VTKM_CONT_EXPORT
 VTKM_CONT
 void InitLogging();

vtkm/filter/scalar_topology/ContourTreeUniformDistributed.cxx

@@ -169,6 +169,7 @@ ContourTreeUniformDistributed::ContourTreeUniformDistributed(vtkm::cont::LogLeve
   : UseBoundaryExtremaOnly(true)
   , UseMarchingCubes(false)
   , AugmentHierarchicalTree(false)
+  , PresimplifyThreshold(0)
   , SaveDotFiles(false)
   , TimingsLogLevel(timingsLogLevel)
   , TreeLogLevel(treeLogLevel)
@@ -585,7 +586,10 @@ inline VTKM_CONT void ContourTreeUniformDistributed::ComputeVolumeMetric(
   vtkmdiy::RegularSwapPartners& partners,
   const FieldType&, // dummy parameter to get the type
   std::stringstream& timingsStream,
-  std::vector<vtkm::cont::DataSet>& hierarchicalTreeOutputDataSet)
+  const vtkm::cont::PartitionedDataSet& input,
+  bool useAugmentedTree,
+  std::vector<vtkm::cont::ArrayHandle<vtkm::Id>>& intrinsicVolumes,
+  std::vector<vtkm::cont::ArrayHandle<vtkm::Id>>& dependentVolumes)
 {
   // TODO/FIXME: CONSIDER MOVING CONTENTS OF THIS METHOD TO SEPARATE FILTER
   vtkm::cont::Timer timer;
@@ -610,30 +614,35 @@ inline VTKM_CONT void ContourTreeUniformDistributed::ComputeVolumeMetric(
   inputContourTreeMaster.foreach (
     [&](DistributedContourTreeBlockData* currInBlock, const vtkmdiy::Master::ProxyWithLink&) {
       vtkm::Id blockNo = currInBlock->LocalBlockNo;
-      const vtkm::cont::DataSet& currDS = hierarchicalTreeOutputDataSet[blockNo];
+      //const vtkm::cont::DataSet& currDS = hierarchicalTreeOutputDataSet[blockNo];
+      auto currOriginalBlock = input.GetPartition(static_cast<vtkm::Id>(blockNo));
       // The block size and origin may be modified during the FanIn so we need to use the
       // size and origin from the original decomposition instead of looking it up in the currInBlock
       vtkm::Id3 pointDimensions, globalPointDimensions, globalPointIndexStart;
-      currDS.GetCellSet().CastAndCallForTypes<VTKM_DEFAULT_CELL_SET_LIST_STRUCTURED>(
+      currOriginalBlock.GetCellSet().CastAndCallForTypes<VTKM_DEFAULT_CELL_SET_LIST_STRUCTURED>(
         vtkm::worklet::contourtree_augmented::GetLocalAndGlobalPointDimensions(),
         pointDimensions,
         globalPointDimensions,
         globalPointIndexStart);
-
       // NOTE: Use dummy link to make DIY happy. The dummy link is never used, since all
       //       communication is via RegularDecomposer, which sets up its own links. No need
       //       to keep the pointer, as DIY will "own" it and delete it when no longer needed.
       // NOTE: Since we passed a "Destroy" function to DIY master, it will own the local data
       //       blocks and delete them when done.
-      hierarchical_hyper_sweep_master.add(
-        currInBlock->GlobalBlockId,
-        new HyperSweepBlock(blockNo,
-                            currInBlock->GlobalBlockId,
-                            globalPointIndexStart,
-                            pointDimensions,
-                            globalPointDimensions,
-                            *currInBlock->HierarchicalAugmenter.AugmentedTree),
-        new vtkmdiy::Link());
+
+      // If we are pre-simplifying the tree then we need to use the base tree and if we compute the
+      // final volume, then we need to use the augmented tree
+      auto hierarchicalTreeToProcess = useAugmentedTree
+        ? currInBlock->HierarchicalAugmenter.AugmentedTree
+        : currInBlock->HierarchicalAugmenter.BaseTree;
+      hierarchical_hyper_sweep_master.add(currInBlock->GlobalBlockId,
+                                          new HyperSweepBlock(blockNo,
+                                                              currInBlock->GlobalBlockId,
+                                                              globalPointIndexStart,
+                                                              pointDimensions,
+                                                              globalPointDimensions,
+                                                              *hierarchicalTreeToProcess),
+                                          new vtkmdiy::Link());
     });
 
   // Log time to copy the data to the HyperSweepBlock data objects
@@ -675,7 +684,6 @@ inline VTKM_CONT void ContourTreeUniformDistributed::ComputeVolumeMetric(
     vtkm::worklet::contourtree_augmented::mesh_dem::IdRelabeler idRelabeler{ b->Origin,
                                                                              b->Size,
                                                                              b->GlobalSize };
-
     if (b->GlobalSize[2] <= 1)
     {
       vtkm::worklet::contourtree_augmented::DataSetMeshTriangulation2DFreudenthal mesh(
@@ -690,7 +698,6 @@ inline VTKM_CONT void ContourTreeUniformDistributed::ComputeVolumeMetric(
       hyperSweeper.InitializeIntrinsicVertexCount(
         b->HierarchicalContourTree, mesh, idRelabeler, b->IntrinsicVolume);
     }
-
 #ifdef DEBUG_PRINT
     VTKM_LOG_S(this->TreeLogLevel, "Block " << b->GlobalBlockId);
     VTKM_LOG_S(this->TreeLogLevel,
@@ -719,7 +726,7 @@ inline VTKM_CONT void ContourTreeUniformDistributed::ComputeVolumeMetric(
     hyperSweeper.LocalHyperSweep();
 
 #ifdef DEBUG_PRINT
-    VTKM_LOG_S(this->TreeLogLevel, "Block " << b->GlobalBlockId);
+    VTKM_LOG_S(this->Tree LogLevel, "Block " << b->GlobalBlockId);
     VTKM_LOG_S(this->TreeLogLevel,
                b->HierarchicalContourTree.DebugPrint("After local hypersweep", __FILE__, __LINE__));
 #endif
@@ -755,16 +762,14 @@ inline VTKM_CONT void ContourTreeUniformDistributed::ComputeVolumeMetric(
                 << ": " << timer.GetElapsedTime() << " seconds" << std::endl;
   timer.Start();
 
-  // Print & add to output data set
-  //std::vector<vtkm::cont::DataSet> hierarchicalTreeAndVolumeOutputDataSet(localDataBlocks.size());
+  // Add the intrinsic and dependent volumes to the output vectors
+  intrinsicVolumes.resize(inputContourTreeMaster.size());
+  dependentVolumes.resize(inputContourTreeMaster.size());
   hierarchical_hyper_sweep_master.foreach (
     [&](HyperSweepBlock* b, const vtkmdiy::Master::ProxyWithLink&) {
-      vtkm::cont::Field intrinsicVolumeField(
-        "IntrinsicVolume", vtkm::cont::Field::Association::WholeDataSet, b->IntrinsicVolume);
-      hierarchicalTreeOutputDataSet[b->LocalBlockNo].AddField(intrinsicVolumeField);
-      vtkm::cont::Field dependentVolumeField(
-        "DependentVolume", vtkm::cont::Field::Association::WholeDataSet, b->DependentVolume);
-      hierarchicalTreeOutputDataSet[b->LocalBlockNo].AddField(dependentVolumeField);
+      intrinsicVolumes[b->LocalBlockNo] = b->IntrinsicVolume;
+      dependentVolumes[b->LocalBlockNo] = b->DependentVolume;
+
 #ifdef DEBUG_PRINT
       VTKM_LOG_S(this->TreeLogLevel, "Block " << b->GlobalBlockId);
       VTKM_LOG_S(
@@ -779,9 +784,6 @@ inline VTKM_CONT void ContourTreeUniformDistributed::ComputeVolumeMetric(
         "Dependent Volume", b->DependentVolume, -1, volumeStream);
       VTKM_LOG_S(this->TreeLogLevel, volumeStream.str());
 #endif
-      // Log the time for adding hypersweep data to the output dataset
-      timingsStream << "    " << std::setw(38) << std::left << "Create Output Data (Hypersweep)"
-                    << ": " << timer.GetElapsedTime() << " seconds" << std::endl;
     });
 }
 
@@ -1135,14 +1137,60 @@ VTKM_CONT void ContourTreeUniformDistributed::DoPostExecute(
                 << ": " << timer.GetElapsedTime() << " seconds" << std::endl;
   timer.Start();
 
+
+  // Compute the volume for pre-simplification if we want to pre-simplify
+  // The dependent volumes from the unaugemented hierarchical tree are used for the pre-simplification
+  // as part of HierarchicalAugmenter.Initialize.
+  std::vector<vtkm::cont::ArrayHandle<vtkm::Id>> unaugmentedDependentVolumes;
+  if (this->PresimplifyThreshold > 0)
+  {
+    // we don't need the unaugemented intrinsic Volumes for the pre-simplification, so we
+    // use a local variable that is being deleted automatically after the context
+    std::vector<vtkm::cont::ArrayHandle<vtkm::Id>> unaugmentedIntrinsicVolumes;
+    // Compute the volume for the base hierarchical tree before augmentation in order to allow for pre-simplification.
+    this->ComputeVolumeMetric(
+      master,
+      assigner,
+      partners,
+      FieldType{},
+      timingsStream,
+      input,
+      false, //   use the unaugmented hierarchical tree (i.e., the base tree) for the volume computation
+      unaugmentedIntrinsicVolumes,
+      unaugmentedDependentVolumes);
+    timingsStream << "    " << std::setw(38) << std::left << "Compute Volume for Presimplication"
+                  << ": " << timer.GetElapsedTime() << " seconds" << std::endl;
+    timer.Start();
+  }
+
   // ******** 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);
-      });
+    vtkm::Id localPresimplifyThreshold = this->PresimplifyThreshold;
+    master.foreach ([globalPointDimensions, localPresimplifyThreshold, unaugmentedDependentVolumes](
+                      DistributedContourTreeBlockData* blockData,
+                      const vtkmdiy::Master::ProxyWithLink&) {
+      // if we don't presimplify then use a NULL pointer for the dependent volume used for pre-simplification
+      vtkm::worklet::contourtree_augmented::IdArrayType* volumeArrayForPresimplifiction = NULL;
+      // if we presimplify then get a pointer for the dependent volume for the current block
+      if (localPresimplifyThreshold > 0)
+      {
+        volumeArrayForPresimplifiction =
+          const_cast<vtkm::worklet::contourtree_augmented::IdArrayType*>(
+            &unaugmentedDependentVolumes[blockData->LocalBlockNo]);
+      }
+      // Initialize the hierarchical augmenter
+      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
+        volumeArrayForPresimplifiction, // DependentVolume if we computed it or NULL if no presimplification is used
+        localPresimplifyThreshold // presimplify if threshold is > 0
+      );
+    });
 
     timingsStream << "    " << std::setw(38) << std::left << "Initalize Hierarchical Trees"
                   << ": " << timer.GetElapsedTime() << " seconds" << std::endl;
@@ -1259,8 +1307,34 @@ VTKM_CONT void ContourTreeUniformDistributed::DoPostExecute(
 
   if (this->AugmentHierarchicalTree)
   {
-    this->ComputeVolumeMetric(
-      master, assigner, partners, FieldType{}, timingsStream, hierarchicalTreeOutputDataSet);
+    std::vector<vtkm::cont::ArrayHandle<vtkm::Id>> augmentedIntrinsicVolumes;
+    std::vector<vtkm::cont::ArrayHandle<vtkm::Id>> augmentedDependentVolumes;
+    this->ComputeVolumeMetric(master,
+                              assigner,
+                              partners,
+                              FieldType{},
+                              timingsStream,
+                              input,
+                              true, //   use the augmented tree
+                              augmentedIntrinsicVolumes,
+                              augmentedDependentVolumes);
+    timer.Start();
+
+    master.foreach (
+      [&](DistributedContourTreeBlockData* blockData, const vtkmdiy::Master::ProxyWithLink&) {
+        // Add the intrinsic and dependent volumes to the output data set
+        vtkm::cont::Field intrinsicVolumeField("IntrinsicVolume",
+                                               vtkm::cont::Field::Association::WholeDataSet,
+                                               augmentedIntrinsicVolumes[blockData->LocalBlockNo]);
+        hierarchicalTreeOutputDataSet[blockData->LocalBlockNo].AddField(intrinsicVolumeField);
+        vtkm::cont::Field dependentVolumeField("DependentVolume",
+                                               vtkm::cont::Field::Association::WholeDataSet,
+                                               augmentedDependentVolumes[blockData->LocalBlockNo]);
+        hierarchicalTreeOutputDataSet[blockData->LocalBlockNo].AddField(dependentVolumeField);
+        // Log the time for adding hypersweep data to the output dataset
+        timingsStream << "    " << std::setw(38) << std::left << "Add Volume Output Data"
+                      << ": " << timer.GetElapsedTime() << " seconds" << std::endl;
+      });
   }
 
   VTKM_LOG_S(this->TimingsLogLevel,

vtkm/filter/scalar_topology/ContourTreeUniformDistributed.h

@@ -118,6 +118,11 @@ public:
     this->AugmentHierarchicalTree = augmentHierarchicalTree;
   }
 
+  VTKM_CONT void SetPresimplifyThreshold(vtkm::Id presimplifyThreshold)
+  {
+    this->PresimplifyThreshold = presimplifyThreshold;
+  }
+
   VTKM_CONT void SetBlockIndices(vtkm::Id3 blocksPerDim,
                                  const vtkm::cont::ArrayHandle<vtkm::Id3>& localBlockIndices)
   {
@@ -127,6 +132,8 @@ public:
 
   VTKM_CONT bool GetAugmentHierarchicalTree() { return this->AugmentHierarchicalTree; }
 
+  VTKM_CONT vtkm::Id GetPresimplifyThreshold() { return this->PresimplifyThreshold; }
+
   VTKM_CONT void SetSaveDotFiles(bool saveDotFiles) { this->SaveDotFiles = saveDotFiles; }
 
   VTKM_CONT bool GetSaveDotFiles() { return this->SaveDotFiles; }
@@ -166,7 +173,10 @@ private:
     vtkmdiy::RegularSwapPartners& partners,
     const FieldType&, // dummy parameter to get the type
     std::stringstream& timingsStream,
-    std::vector<vtkm::cont::DataSet>& hierarchicalTreeOutputDataSet);
+    const vtkm::cont::PartitionedDataSet& input,
+    bool useAugmentedTree,
+    std::vector<vtkm::cont::ArrayHandle<vtkm::Id>>& intrinsicVolumes,
+    std::vector<vtkm::cont::ArrayHandle<vtkm::Id>>& dependentVolumes);
 
   ///
   /// Internal helper function that implements the actual functionality of PostExecute
@@ -188,6 +198,9 @@ private:
   /// Augment hierarchical tree
   bool AugmentHierarchicalTree;
 
+  /// Threshold to use for volume pre-simplification
+  vtkm::Id PresimplifyThreshold;
+
   /// Save dot files for all tree computations
   bool SaveDotFiles;
 

vtkm/filter/scalar_topology/worklet/branch_decomposition/HierarchicalVolumetricBranchDecomposer.h

@@ -506,12 +506,19 @@ inline void HierarchicalVolumetricBranchDecomposer::CollapseBranches(
     vtkm::worklet::contourtree_distributed::FindSuperArcBetweenNodes findSuperArcBetweenNodes{
       hierarchicalTreeSuperarcs
     };
+    // Get the number of rounds
+    auto numRoundsArray = hierarchicalTreeDataSet.GetField("NumRounds")
+                            .GetData()
+                            .AsArrayHandle<vtkm::cont::ArrayHandle<vtkm::Id>>();
+    vtkm::Id numRounds = vtkm::cont::ArrayGetValue(0, numRoundsArray);
 
     using vtkm::worklet::scalar_topology::hierarchical_volumetric_branch_decomposer::
       CollapseBranchesWorklet;
-    this->Invoke(CollapseBranchesWorklet{},         // the worklet
+    CollapseBranchesWorklet collapseBranchesWorklet(numRounds);
+    this->Invoke(collapseBranchesWorklet,           // the worklet
                  this->BestUpSupernode,             // input
                  this->BestDownSupernode,           // input
+                 hierarchicalTreeSuperarcs,         // input
                  findRegularByGlobal,               // input ExecutionObject
                  findSuperArcBetweenNodes,          // input ExecutionObject
                  hierarchicalTreeRegular2Supernode, // input

vtkm/filter/scalar_topology/worklet/branch_decomposition/hierarchical_volumetric_branch_decomposer/CollapseBranchesWorklet.h

@@ -65,6 +65,7 @@ public:
   using ControlSignature = void(
     FieldIn bestUpSupernode,
     FieldIn bestDownSupernode,
+    FieldIn superarcs,
     // Execution objects from the hierarchical tree to use the FindRegularByGlobal function
     ExecObject findRegularByGlobal,
     // Execution objects from the hierarchical tree to use the FindSuperArcBetweenNodes, function
@@ -72,12 +73,15 @@ public:
     WholeArrayIn hierarchicalTreeRegular2supernode,
     WholeArrayIn hierarchicalTreeWhichRound,
     WholeArrayInOut branchRoot);
-  using ExecutionSignature = void(InputIndex, _1, _2, _3, _4, _5, _6, _7);
+  using ExecutionSignature = void(InputIndex, _1, _2, _3, _4, _5, _6, _7, _8);
   using InputDomain = _1;
 
   /// Default Constructor
   VTKM_EXEC_CONT
-  CollapseBranchesWorklet() {}
+  CollapseBranchesWorklet(vtkm::Id numRounds)
+    : NumRounds(numRounds)
+  {
+  }
 
   /// operator() of the workelt
   template <typename ExecObjectType1,
@@ -88,6 +92,7 @@ public:
     const vtkm::Id& supernode,                  // iteration index
     const vtkm::Id& bestUpSupernodeId,          // bestUpSupernode[supernode]
     const vtkm::Id& bestDownSupernodeId,        // bestDownSupernode[supernode]
+    const vtkm::Id& superarcsId,                // hierarchicalTree.superarcs[supernode]
     const ExecObjectType1& findRegularByGlobal, // Execution object to call FindRegularByGlobal
     const ExecObjectType2&
       findSuperArcBetweenNodes, // Execution object to call FindSuperArcBetweenNodes
@@ -104,6 +109,18 @@ public:
     //          If it does exist and is an upwards superarc, then the current supernode must have an ascending arc to it, and we're done
     //  Also do the same for the best down, then for each supernode, point the higher numbered at the lower
 
+    // ADDED 19/07/2023
+    // If there are any attachment points left in the hierarchical tree, there is an extra edge case we need to deal with.
+    // It occurs when a supernode is simultaneously the target of an ascending superarc and a descending one
+    // What we do is to test for this here: if we are an attachment point, we omit connecting the best up and down
+    // ADDED 19/07/2023
+    // test for attachment points
+    if ((hierarchicalTreeWhichRoundPortal.Get(supernode) != this->NumRounds) &&
+        (vtkm::worklet::contourtree_augmented::NoSuchElement(superarcsId)))
+    {
+      return;
+    }
+
     // if there is no best up, we're at an upper leaf and will not connect up two superarcs anyway, so we can skip the supernode
     if (vtkm::worklet::contourtree_augmented::NoSuchElement(bestUpSupernodeId))
     {
@@ -233,6 +250,8 @@ public:
     */
   } // operator()()
 
+private:
+  vtkm::Id NumRounds;
 
 }; // CollapseBranchesWorklet
 

vtkm/filter/scalar_topology/worklet/contourtree_augmented/Types.h

@@ -86,12 +86,20 @@ constexpr vtkm::Id IS_REGULAR = static_cast<vtkm::Id>(2);
 constexpr vtkm::Id IS_SADDLE = static_cast<vtkm::Id>(3);
 constexpr vtkm::Id IS_ATTACHMENT = static_cast<vtkm::Id>(4);
 
+// TERMINAL_ELEMENT is primarily used for optimisation of memory access during pointer doubling operations
+// We now need to distinguish between a supernode and superarc when sorting by superarc(node) IDs
+// This only (at present) comes up when processing attachment points, which have null superarcs, so it
+// is reasonable to reuse TERMINAL_ELEMENT for this purpose.  However, we give it a separate macro name with
+// the same value to aid comprehension
+constexpr vtkm::Id TRANSFER_TO_SUPERARC = TERMINAL_ELEMENT;
+
 // clang-format on
 using IdArrayType = vtkm::cont::ArrayHandle<vtkm::Id>;
 
 using EdgePair = vtkm::Pair<vtkm::Id, vtkm::Id>; // here EdgePair.first=low and EdgePair.second=high
 using EdgePairArray = vtkm::cont::ArrayHandle<EdgePair>; // Array of edge pairs
 
+
 // inline functions for retrieving flags or index
 VTKM_EXEC_CONT
 inline bool NoSuchElement(vtkm::Id flaggedIndex)
@@ -143,6 +151,12 @@ inline bool NoFlagsSet(vtkm::Id flaggedIndex)
   return (flaggedIndex & ~INDEX_MASK) == 0;
 } // NoFlagsSet
 
+// Helper function: to check that the TRANSFER_TO_SUPERARC flag is set
+VTKM_EXEC_CONT
+inline bool TransferToSuperarc(vtkm::Id flaggedIndex)
+{ // transferToSuperarc()
+  return ((flaggedIndex & TRANSFER_TO_SUPERARC) != 0);
+} // transferToSuperarc()
 
 // Debug helper function: Assert that an index array has no element with any flags set
 template <typename S>
@@ -225,6 +239,23 @@ inline std::string FlagString(vtkm::Id flaggedIndex)
   return fString;
 } // FlagString()
 
+
+// == comparison operator for edges
+inline bool edgeEqual(const EdgePair& LHS, const EdgePair& RHS)
+{ // operator ==
+
+  if (LHS.first != RHS.first)
+  {
+    return false;
+  }
+  if (LHS.second != RHS.second)
+  {
+    return false;
+  }
+  return true;
+} // operator ==
+
+
 class EdgeDataHeight
 {
 public:

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

@@ -101,6 +101,9 @@
 #include <vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_augmenter/AttachmentIdsEqualComparator.h>
 #include <vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_augmenter/AttachmentSuperparentAndIndexComparator.h>
 #include <vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_augmenter/CopyBaseRegularStructureWorklet.h>
+#include <vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_augmenter/CreateSuperarcsData.h>
+#include <vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_augmenter/CreateSuperarcsSetFirstSupernodePerIterationWorklet.h>
+#include <vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_augmenter/CreateSuperarcsUpdateFirstSupernodePerIterationWorklet.h>
 #include <vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_augmenter/CreateSuperarcsWorklet.h>
 #include <vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_augmenter/FindSuperparentForNecessaryNodesWorklet.h>
 #include <vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_augmenter/HierarchicalAugmenterInOutData.h>
@@ -127,6 +130,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
@@ -173,7 +181,7 @@ public:
 
   /// these are essentially temporary local variables, but are placed here to make the DebugPrint()
   /// more comprehensive. They will be allocated where used
-  /// this one makes a list of attachment Ids and is used in sevral different places, so resize it when done
+  /// this one makes a list of attachment Ids and is used in several different places, so resize it when done
   vtkm::worklet::contourtree_augmented::IdArrayType AttachmentIds;
   /// tracks segments of attachment points by round
   vtkm::worklet::contourtree_augmented::IdArrayType FirstAttachmentPointInRound;
@@ -198,7 +206,12 @@ 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,
+    vtkm::worklet::contourtree_augmented::IdArrayType* volumeArray = NULL,
+    vtkm::Id presimplifyThreshold = 0);
 
   /// routine to prepare the set of attachment points to transfer
   void PrepareOutAttachmentPoints(vtkm::Id round);
@@ -220,7 +233,7 @@ public:
   /// transfer level of superstructure with insertions
   void CopySuperstructure();
   /// reset the super Ids in the hyperstructure to the new values
-  void UpdateHyperstructure();
+  void UpdateHyperstructure(vtkm::Id roundNumber);
   /// copy the remaining base level regular nodes
   void CopyBaseRegularStructure();
 
@@ -249,22 +262,36 @@ 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,
+  vtkm::worklet::contourtree_augmented::IdArrayType* volumeArray,
+  vtkm::Id presimplifyThreshold)
 { // 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
+  // except those under the presimplify threshold. The presimplification is
+  // handled in the IsAttachementPointPredicate
+
   // to do this, we construct an index array with all supernode ID's that satisfy:
   // 1. superparent == NO_SUCH_ELEMENT (i.e. root of interior tree)
   // 2. round < nRounds (except the top level, where 1. indicates the tree root)
-  // initalize AttachementIds
+  // initialize AttachementIds
   {
     vtkm::worklet::contourtree_distributed::hierarchical_augmenter::IsAttachementPointPredicate
-      isAttachementPointPredicate(
-        this->BaseTree->Superarcs, this->BaseTree->WhichRound, this->BaseTree->NumRounds);
+      isAttachementPointPredicate(this->BaseTree->Superarcs,
+                                  this->BaseTree->WhichRound,
+                                  this->BaseTree->NumRounds,
+                                  volumeArray,
+                                  presimplifyThreshold);
     auto tempSupernodeIndex =
       vtkm::cont::ArrayHandleIndex(this->BaseTree->Supernodes.GetNumberOfValues());
     vtkm::cont::Algorithm::CopyIf(
@@ -476,7 +503,11 @@ void HierarchicalAugmenter<FieldType>::ReleaseSwapArrays()
 } // ReleaseSwapArrays()
 
 
-// routine to reconstruct a hierarchical tree using the augmenting supernodes
+// routine to reconstruct a hierarchical tree using the augmenting supernodes.
+// Allowing pre-simplification needs the superstructure and hyperstructure to be done one layer
+// at a time. This means lifting the code from CopySuperstructure up to here CopyHyperstructure()
+// can actually be left the way it is - copying the entire hyperstructure,  as the augmentation
+// process doesn't change it.  It might be sensible to change it anyway, just to make the code better organised
 template <typename FieldType>
 void HierarchicalAugmenter<FieldType>::BuildAugmentedTree()
 { // BuildAugmentedTree()
@@ -484,11 +515,22 @@ void HierarchicalAugmenter<FieldType>::BuildAugmentedTree()
   this->PrepareAugmentedTree();
   // 2.  Copy the hyperstructure, using the old super IDs for now
   this->CopyHyperstructure();
-  // 3.  Copy the superstructure, inserting additional points as we do
-  this->CopySuperstructure();
-  // 4.  Update the hyperstructure to use the new super IDs
-  this->UpdateHyperstructure();
-  // 5.  Copy the remaining regular structure at the bottom level, setting up the regular sort order in the process
+  // 3.	Copy superstructure one round at a time, updating the hyperstructure as well
+  // (needed to permit search for superarcs)
+  //	Loop from the top down:
+  for (vtkm::Id roundNumber = this->BaseTree->NumRounds; roundNumber >= 0; roundNumber--)
+  { // per round
+    // start by retrieving list of old supernodes from the tree (except for attachment points)
+    this->RetrieveOldSupernodes(roundNumber);
+    // since we know the number of attachment points, we can now allocate space for the level
+    // and set up arrays for sorting the supernodes
+    this->ResizeArrays(roundNumber);
+    // now we create the superarcs for the round in the new tree
+    this->CreateSuperarcs(roundNumber);
+    // finally, we update the hyperstructure for the round in the new tree
+    this->UpdateHyperstructure(roundNumber);
+  } // per round
+  // 4.  Copy the remaining regular structure at the bottom level, setting up the regular sort order in the process
   this->CopyBaseRegularStructure();
 } // BuildAugmentedTree()
 
@@ -623,53 +665,91 @@ void HierarchicalAugmenter<FieldType>::CopyHyperstructure()
                                 this->AugmentedTree->FirstHypernodePerIteration[roundNumber]);
   } // per round
 
+  // WARNING 28/05/2023:  Since this resize is for the full hyperstructure, it should be safe to put here.
+  // Unless of course, anything relies on the sizes: but they were 0, so it is unlikely.
+  // A search for hyperarcs.size() & hypernodes.size() in this unit confirmed that nothing uses them.
+  // Nevertheless, set them all to NO_SUCH_ELEMENT out of paranoia
+  // 5. Reset hypernodes, hyperarcs and superchildren using supernode IDs
+  //    The hyperstructure is unchanged, but uses old supernode IDs
+  vtkm::worklet::contourtree_augmented::ResizeVector<vtkm::Id>(
+    this->AugmentedTree->Hypernodes,                      // resize array
+    this->BaseTree->Hypernodes.GetNumberOfValues(),       // new size
+    vtkm::worklet::contourtree_augmented::NO_SUCH_ELEMENT // set all elements to this value
+  );
+  vtkm::worklet::contourtree_augmented::ResizeVector<vtkm::Id>(
+    this->AugmentedTree->Hyperarcs,                       // resize array
+    this->BaseTree->Hyperarcs.GetNumberOfValues(),        // new size
+    vtkm::worklet::contourtree_augmented::NO_SUCH_ELEMENT // set all elements to this value
+  );
+  vtkm::worklet::contourtree_augmented::ResizeVector<vtkm::Id>(
+    this->AugmentedTree->Superchildren,                // resize array
+    this->BaseTree->Superchildren.GetNumberOfValues(), // new size
+    0                                                  // set all elements to this value
+  );
+
 #ifdef DEBUG_PRINT
   VTKM_LOG_S(vtkm::cont::LogLevel::Info, DebugPrint("Hyperstructure Copied", __FILE__, __LINE__));
 #endif
 } // CopyHyperstructure()
 
 
+// WARNING 28/05/2023: deleted and moved up to AugmentTree() routine
 // transfer level of superstructure with insertions
-template <typename FieldType>
-void HierarchicalAugmenter<FieldType>::CopySuperstructure()
-{ // CopySuperstructure()
-
-  //  Loop from the top down:
-  for (vtkm::Id roundNumber = this->BaseTree->NumRounds; roundNumber >= 0; roundNumber--)
-  { // per round
-    // start by retrieving list of old supernodes from the tree (except for attachment points)
-    this->RetrieveOldSupernodes(roundNumber);
-    // since we know the number of attachment points, we can now allocate space for the level
-    // and set up arrays for sorting the supernodes
-    this->ResizeArrays(roundNumber);
-    // now we create the superarcs for the round in the new tree
-    this->CreateSuperarcs(roundNumber);
-  } // per round
-#ifdef DEBUG_PRINT
-  VTKM_LOG_S(vtkm::cont::LogLevel::Info,
-             this->DebugPrint("Superstructure Copied", __FILE__, __LINE__));
-#endif
-} // CopySuperstructure()
+//template <typename FieldType>
+//void HierarchicalAugmenter<FieldType>::CopySuperstructure()
+//{ // CopySuperstructure()
+//
+//  //  Loop from the top down:
+//  for (vtkm::Id roundNumber = this->BaseTree->NumRounds; roundNumber >= 0; roundNumber--)
+//  { // per round
+//    // start by retrieving list of old supernodes from the tree (except for attachment points)
+//    this->RetrieveOldSupernodes(roundNumber);
+//    // since we know the number of attachment points, we can now allocate space for the level
+//    // and set up arrays for sorting the supernodes
+//    this->ResizeArrays(roundNumber);
+//    // now we create the superarcs for the round in the new tree
+//    this->CreateSuperarcs(roundNumber);
+//  } // per round
+//#ifdef DEBUG_PRINT
+//  VTKM_LOG_S(vtkm::cont::LogLevel::Info,
+//             this->DebugPrint("Superstructure Copied", __FILE__, __LINE__));
+//#endif
+//} // CopySuperstructure()
 
 
 // reset the super IDs in the hyperstructure to the new values
 template <typename FieldType>
-void HierarchicalAugmenter<FieldType>::UpdateHyperstructure()
+void HierarchicalAugmenter<FieldType>::UpdateHyperstructure(vtkm::Id roundNumber)
 { // UpdateHyperstructure()
-
-  //  5.  Reset hypernodes, hyperarcs and superchildren using supernode IDs
-  //      The hyperstructure is unchanged, but uses old supernode IDs
+  // now that the superstructure is known, we can find the new supernode IDs for all
+  // of the old hypernodes at this level and update.Wwe want to update the entire round
+  // at once, so we would like to use the firstHypernodePerIteration array.
+  vtkm::Id startIndex =
+    vtkm::cont::ArrayGetValue(0, this->AugmentedTree->FirstHypernodePerIteration[roundNumber]);
+  vtkm::Id stopIndex = vtkm::cont::ArrayGetValue(
+    vtkm::cont::ArrayGetValue(roundNumber, this->AugmentedTree->NumIterations),
+    this->AugmentedTree->FirstHypernodePerIteration[roundNumber]);
+  vtkm::Id selectSize = stopIndex - startIndex;
   {
     vtkm::worklet::contourtree_distributed::hierarchical_augmenter::
       UpdateHyperstructureSetHyperarcsAndNodesWorklet
         updateHyperstructureSetHyperarcsAndNodesWorklet;
-    this->Invoke(
-      updateHyperstructureSetHyperarcsAndNodesWorklet,
-      this->BaseTree->Hypernodes,      // input
-      this->BaseTree->Hyperarcs,       // input
-      this->NewSupernodeIds,           // input
-      this->AugmentedTree->Hypernodes, // output (the array is automatically resized here)
-      this->AugmentedTree->Hyperarcs   // output (the array is automatically resized here)
+    // Create ArrayHandleViews of the subrange of the input and output arrays we need to process
+    auto baseTreeHypernodesView =
+      vtkm::cont::make_ArrayHandleView(this->BaseTree->Hypernodes, startIndex, selectSize);
+    auto baseTreeHyperarcsView =
+      vtkm::cont::make_ArrayHandleView(this->BaseTree->Hyperarcs, startIndex, selectSize);
+    auto augmentedTreeHypernodesView =
+      vtkm::cont::make_ArrayHandleView(this->AugmentedTree->Hypernodes, startIndex, selectSize);
+    auto augmentedTreeHyperarcsView =
+      vtkm::cont::make_ArrayHandleView(this->AugmentedTree->Hyperarcs, startIndex, selectSize);
+    // Invoke the worklet with the subset view of our arrays
+    this->Invoke(updateHyperstructureSetHyperarcsAndNodesWorklet,
+                 baseTreeHypernodesView,      // input
+                 baseTreeHyperarcsView,       // input
+                 this->NewSupernodeIds,       // input
+                 augmentedTreeHypernodesView, // output (the array is automatically resized here)
+                 augmentedTreeHyperarcsView   // output (the array is automatically resized here)
     );
   }
 
@@ -679,10 +759,20 @@ void HierarchicalAugmenter<FieldType>::UpdateHyperstructure()
     vtkm::worklet::contourtree_distributed::hierarchical_augmenter::
       UpdateHyperstructureSetSuperchildrenWorklet updateHyperstructureSetSuperchildrenWorklet(
         this->AugmentedTree->Supernodes.GetNumberOfValues());
-    this->Invoke(
-      updateHyperstructureSetSuperchildrenWorklet,
-      this->AugmentedTree->Hypernodes,   // input
-      this->AugmentedTree->Superchildren // output (the array is automatically resized here)
+    // As above, we need to create views of the relevant subranges of our arrays
+    auto augmentedTreeSuperarcsView =
+      vtkm::cont::make_ArrayHandleView(this->AugmentedTree->Superarcs, startIndex, selectSize);
+    vtkm::Id extraSelectSize =
+      ((startIndex + selectSize) < this->AugmentedTree->Hyperparents.GetNumberOfValues())
+      ? selectSize + 1
+      : selectSize;
+    auto augmentedTreeHyperparentsView = vtkm::cont::make_ArrayHandleView(
+      this->AugmentedTree->Hyperparents, startIndex, extraSelectSize);
+    this->Invoke(updateHyperstructureSetSuperchildrenWorklet,
+                 this->AugmentedTree->Hypernodes,   // input
+                 augmentedTreeSuperarcsView,        // input
+                 augmentedTreeHyperparentsView,     // inpu
+                 this->AugmentedTree->Superchildren // output
     );
   }
 } // UpdateHyperstructure()
@@ -728,12 +818,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 +922,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::
@@ -850,23 +947,22 @@ void HierarchicalAugmenter<FieldType>::RetrieveOldSupernodes(vtkm::Id roundNumbe
   //    TODO PERFORMANCE STATISTICS:
   //      the # of supernodes at each level minus the # of kept supernodes gives us the # of attachment points we lose at this level
   //      in addition to this, the firstAttachmentPointInRound array gives us the # of attachment points we gain at this level
+  //
+  // COMMENT: In adding presimplification, this will have to change.
+  // Previously, it made the hard assumption that all attachment points were transferred & used that to suppress
+  // them. Now it can do that no longer, which gives us a choice.  We can pass in the threshold & volume array
+  // and test here, but that's not ideal as it does the same test in multiple places. Alternatively, we can
+  // look up whether the supernode is already present in the structure, which has an associated search cost.
+  // BUT, we have an array called newSupernodeIDs for the purpose already, so that's how we'll do it.
+
   vtkm::Id supernodeIndexBase =
     vtkm::cont::ArrayGetValue(0, this->BaseTree->FirstSupernodePerIteration[roundNumber]);
   vtkm::cont::ArrayHandleCounting<vtkm::Id> supernodeIdVals(
     supernodeIndexBase, // start
     1,                  // step
     this->BaseTree->NumSupernodesInRound.ReadPortal().Get(roundNumber));
-  // the test for whether to keep it is:
-  // a1. at the top level, keep everything
-  if (!(roundNumber < this->BaseTree->NumRounds))
-  {
-    vtkm::cont::Algorithm::Copy(supernodeIdVals, this->KeptSupernodes);
-  }
-  // a2. at lower levels, keep them if the superarc is NO_SUCH_ELEMENT
-  else
   {
     // Reset this-KeptSupernodes to the right size and initalize with NO_SUCH_ELEMENT.
-    // TODO: Check if a simple free and allocate without initalizing the array is sufficient
     vtkm::cont::Algorithm::Copy(
       // Create const array to copy
       vtkm::cont::ArrayHandleConstant<vtkm::Id>(
@@ -882,7 +978,7 @@ void HierarchicalAugmenter<FieldType>::RetrieveOldSupernodes(vtkm::Id roundNumbe
       supernodeIdVals,
       // Stencil with baseTree->superarcs[supernodeID]
       vtkm::cont::make_ArrayHandleView(
-        this->BaseTree->Superarcs, supernodeIndexBase, this->KeptSupernodes.GetNumberOfValues()),
+        this->NewSupernodeIds, supernodeIndexBase, this->KeptSupernodes.GetNumberOfValues()),
       // And the CopyIf compresses the array to eliminate unnecssary elements
       // save to this->KeptSupernodes
       this->KeptSupernodes,
@@ -908,6 +1004,10 @@ template <typename FieldType>
 void HierarchicalAugmenter<FieldType>::ResizeArrays(vtkm::Id roundNumber)
 { // ResizeArrays()
   // at this point, we know how many supernodes are kept from the same level of the old tree
+  // TODO/WARNING 23/07/2023:
+  // Actually, this is no longer true.  If the same vertex is an attachment point for two adjacent blocks (i.e. it is on the boundary), it is entirely possible
+  // for one block to add it at a higher level than the other.  To preclude this, we will need to edit RetrieveOldSupernodes()
+
   // we can also find out how many supernodes are being inserted, which gives us the correct amount to expand by, saving a double resize() call
   // note that some of these arrays could probably be resized later, but it's cleaner this way
   // also note that if it becomes a problem, we could resize all of the arrays to baseTree->supernodes.size() + # of attachmentPoints as an over-estimate
@@ -1013,6 +1113,9 @@ void HierarchicalAugmenter<FieldType>::ResizeArrays(vtkm::Id roundNumber)
 #endif
 
   // b. Transfer attachment points for level into new supernode array
+  // NOTE: this means the set of attachment points that we have determined by swapping
+  //       need to be inserted onto a superarc at this level. All of them should be from
+  //       lower levels originally, but are being moved up to this level for insertion
   // to copy them in, we use the existing array of attachment point IDs by round
   {
     vtkm::Id firstAttachmentPointInRoundCurrent =
@@ -1061,7 +1164,10 @@ void HierarchicalAugmenter<FieldType>::ResizeArrays(vtkm::Id roundNumber)
                         __LINE__));
 #endif
 
-  // Now we copy in the kept supernodes
+  // Now we copy in the kept supernodes: this used to mean only the non-attachment points
+  // now it includes the attachment points at this level that the simplification removed
+  // so they need to be put back where they were
+  // However, that means that all of them do exist in the base tree, so we can copy from there
   {
     auto oldRegularIdArr =
       vtkm::cont::make_ArrayHandlePermutation(this->KeptSupernodes,        // index
@@ -1135,15 +1241,24 @@ void HierarchicalAugmenter<FieldType>::ResizeArrays(vtkm::Id roundNumber)
       ResizeArraysBuildNewSupernodeIdsWorklet resizeArraysBuildNewSupernodeIdsWorklet(
         numSupernodesAlready);
     auto supernodeIndex = vtkm::cont::ArrayHandleIndex(this->SupernodeSorter.GetNumberOfValues());
-    auto supernodeIdSetPermuted =
-      vtkm::cont::make_ArrayHandlePermutation(this->SupernodeSorter, this->SupernodeIdSet);
+    // auto supernodeIdSetPermuted =
+    //  vtkm::cont::make_ArrayHandlePermutation(this->SupernodeSorter, this->SupernodeIdSet);
+    auto globalRegularIdSetPermuted =
+      vtkm::cont::make_ArrayHandlePermutation(this->SupernodeSorter, this->GlobalRegularIdSet);
+    auto findRegularByGlobal = this->BaseTree->GetFindRegularByGlobal();
     this->Invoke(
       resizeArraysBuildNewSupernodeIdsWorklet,
       supernodeIndex, // input domain. We only need the index because supernodeIdSetPermuted already does the permute
-      supernodeIdSetPermuted, // input input supernodeIDSet permuted by supernodeSorter to allow for FieldIn
+      // supernodeIdSetPermuted, // input input supernodeIDSet permuted by supernodeSorter to allow for FieldIn
+      globalRegularIdSetPermuted,
+      findRegularByGlobal,
+      this->BaseTree->Regular2Supernode,
       this
         ->NewSupernodeIds // output/input (both are necessary since not all values will be overwritten)
     );
+
+    // Add const ExecObjectType1& findRegularByGlobal,
+    // Add baseTree->regular2supernode
   }
 
 #ifdef DEBUG_PRINT
@@ -1161,167 +1276,250 @@ template <typename FieldType>
 void HierarchicalAugmenter<FieldType>::CreateSuperarcs(vtkm::Id roundNumber)
 { // CreateSuperarcs()
   // retrieve the ID number of the first supernode at this level
+#ifdef DEBUG_PRINT
+  VTKM_LOG_S(vtkm::cont::LogLevel::Info,
+             DebugPrint(std::string("Round ") + std::to_string(roundNumber) +
+                          std::string(" Starting CreateSuperarcs()"),
+                        __FILE__,
+                        __LINE__));
+#endif
+
+  vtkm::Id currNumIterations =
+    vtkm::cont::ArrayGetValue(roundNumber, this->AugmentedTree->NumIterations);
   vtkm::Id numSupernodesAlready =
     vtkm::cont::ArrayGetValue(0, this->AugmentedTree->FirstSupernodePerIteration[roundNumber]);
 
-  //  e.  Connect superarcs for the level & set hyperparents & superchildren count, whichRound, whichIteration, super2hypernode
-  { // START scope for e. to delete temporary variables
-    // Note:  The original PPP algorithm performed all operations listed in this block
-    //        in a single parralel for loop. Many of those operations were smart array
-    //        copies. So to simplfy the worklet and to make more effective use of
-    //        VTKm algorithm, a large number of copy operations have been extracted from
-    //        the loop and are performed here via combinations of fancy array handles and
-    //        vtkm::cont::Algorithm::Copy operations.
+  // e.	Connect superarcs for the level & set hyperparents & superchildren count, whichRound, whichIteration, super2hypernode
+  // 24/05/2023: Expansion of comment to help debug.
+  // At this point, we know that all higher rounds are correctly constructed, and that any attachment points that survived simplification
+  // have already been inserted in a higher round.
 
-    // Define the new supernode and regular Id. Both are actually the same here since we are
-    // augmenting the tree here, but for clarity we define them as separate variables.
-    // At all levels above 0, we used to keep regular vertices in case
-    // they are attachment points.  After augmentation, we don't need to.
-    // Instead, at all levels above 0, the regular nodes in each round
-    // are identical to the supernodes. In order to avoid confusion,
-    // we will copy the Id into a separate variable
-    vtkm::cont::ArrayHandleCounting<vtkm::Id> newSupernodeId(
-      numSupernodesAlready,                     // start
-      static_cast<vtkm::Id>(1),                 // step
-      this->SupernodeSorter.GetNumberOfValues() // number of values
-    );
-    auto newRegularId = newSupernodeId;
-    // define the superparentOldSuperId
-    auto permutedSuperparentSet =
-      vtkm::cont::make_ArrayHandlePermutation(this->SupernodeSorter, this->SuperparentSet);
-    auto superparentOldSuperId = vtkm::cont::make_ArrayHandleTransform(
-      permutedSuperparentSet, vtkm::worklet::contourtree_augmented::MaskedIndexFunctor<vtkm::Id>());
+  // The sort should have resulted in the supernodes being segmented along old superarcs.  Most supernodes should be in a segment of length
+  // 1, and should be their own superparent in the sort array.  But we can't readily test that, because other supernodes may also have them
+  // as the superparent.
 
-    // set the supernode's regular Id. Set: augmentedTree->supernodes[newSupernodeID] = newRegularID;
-    auto permutedAugmentedTreeSupernodes =
-      vtkm::cont::make_ArrayHandlePermutation(newSupernodeId, this->AugmentedTree->Supernodes);
-    vtkm::cont::Algorithm::Copy(newRegularId, permutedAugmentedTreeSupernodes);
+  // This loop will principally determine the superarc for each supernode.  For this, the rules break down to:
+  // 1.	If the supernode is the global root, connect it nowhere
+  // 2.	If the supernode is the last in the set of all supernodes in this round, treat it as the end of a segment
+  // 3. If the supernode is the last in a segment by superarc, connect it to the target of its superparent in the old tree, using the new supernode ID
+  // 4.	Otherwise, connect to the new supernode ID of the next supernode in the segment
 
-    // Run the worklet for more complex operations
-    { // START block for CreateSuperarcsWorklet
-      // create the worklet
-      vtkm::worklet::contourtree_distributed::hierarchical_augmenter::CreateSuperarcsWorklet
-        createSuperarcsWorklet(
-          numSupernodesAlready,
-          this->BaseTree->NumRounds,
-          vtkm::cont::ArrayGetValue(roundNumber, this->AugmentedTree->NumIterations),
-          roundNumber,
-          this->AugmentedTree->Supernodes.GetNumberOfValues());
-      // create fancy arrays needed to allow use of FieldIn for worklet parameters
-      auto permutedGlobalRegularIdSet =
-        vtkm::cont::make_ArrayHandlePermutation(this->SupernodeSorter, this->GlobalRegularIdSet);
-      auto augmentedTreeSuperarcsView =
-        vtkm::cont::make_ArrayHandleView(this->AugmentedTree->Superarcs,
-                                         numSupernodesAlready,
-                                         this->SupernodeSorter.GetNumberOfValues());
-      auto augmentedTreeSuper2HypernodeView =
-        vtkm::cont::make_ArrayHandleView(this->AugmentedTree->Super2Hypernode,
-                                         numSupernodesAlready,
-                                         this->SupernodeSorter.GetNumberOfValues());
-      // invoke the worklet
-      this->Invoke(createSuperarcsWorklet,                                       // the worklet
-                   this->SupernodeSorter,                                        // input domain
-                   this->SuperparentSet,                                         // input
-                   this->BaseTree->Superarcs,                                    // input
-                   this->NewSupernodeIds,                                        // input
-                   this->BaseTree->Supernodes,                                   // input
-                   this->BaseTree->RegularNodeGlobalIds,                         // input
-                   permutedGlobalRegularIdSet,                                   // input
-                   this->BaseTree->Super2Hypernode,                              // input
-                   this->BaseTree->WhichIteration,                               // input
-                   augmentedTreeSuperarcsView,                                   // output
-                   this->AugmentedTree->FirstSupernodePerIteration[roundNumber], // input/output
-                   augmentedTreeSuper2HypernodeView                              // output
-      );
-    } // END block for CreateSuperarcsWorklet
+  // In each case, we will need to preserve the ascending / descending flag
 
-    // setting the hyperparent is straightforward since the hyperstructure is preserved
-    // we take the superparent (which is guaranteed to be in the baseTree), find it's hyperparent and use that
-    // Set:  augmentedTree->hyperparents[newSupernodeID] = baseTree->hyperparents[superparentOldSuperID];
-    auto permutedAugmentedTreeHyperparents =
-      vtkm::cont::make_ArrayHandlePermutation(newSupernodeId, this->AugmentedTree->Hyperparents);
-    auto permutedBaseTreeHyperparents =
-      vtkm::cont::make_ArrayHandlePermutation(superparentOldSuperId, this->BaseTree->Hyperparents);
-    vtkm::cont::Algorithm::Copy(permutedBaseTreeHyperparents, permutedAugmentedTreeHyperparents);
+  // We will also have to set the first supernode per iteration - if possible, in a separate loop
 
-    // which round and iteration carry over
-    // Set: augmentedTree->whichRound[newSupernodeID] = baseTree->whichRound[superparentOldSuperID];
-    auto permutedAugmentedTreeWhichRound =
-      vtkm::cont::make_ArrayHandlePermutation(newSupernodeId, this->AugmentedTree->WhichRound);
-    auto permutedBaseTreeWhichRound =
-      vtkm::cont::make_ArrayHandlePermutation(superparentOldSuperId, this->BaseTree->WhichRound);
-    vtkm::cont::Algorithm::Copy(permutedBaseTreeWhichRound, permutedAugmentedTreeWhichRound);
+  // We need this to determine which supernodes are inserted and which are attached (see below)
+  vtkm::Id numInsertedSupernodes =
+    (vtkm::cont::ArrayGetValue(roundNumber + 1, this->FirstAttachmentPointInRound) -
+     vtkm::cont::ArrayGetValue(roundNumber, this->FirstAttachmentPointInRound));
 
-    // Set: augmentedTree->whichIteration[newSupernodeID] = baseTree->whichIteration[superparentOldSuperID];
-    auto permutedAugmentedTreeWhichIteration =
-      vtkm::cont::make_ArrayHandlePermutation(newSupernodeId, this->AugmentedTree->WhichIteration);
-    auto permutedBaseTreeWhichIterationPortal = vtkm::cont::make_ArrayHandlePermutation(
-      superparentOldSuperId, this->BaseTree->WhichIteration);
-    vtkm::cont::Algorithm::Copy(permutedBaseTreeWhichIterationPortal,
-                                permutedAugmentedTreeWhichIteration);
+  { // START Call CreateSuperarcsWorklet (scope to delete temporary variables)
+    // create the worklet
+    vtkm::worklet::contourtree_distributed::hierarchical_augmenter::CreateSuperarcsWorklet<
+      FieldType>
+      createSuperarcsWorklet(
+        numSupernodesAlready, this->BaseTree->NumRounds, numInsertedSupernodes, roundNumber);
 
-    // now we deal with the regular-sized arrays. In the following supernodeSetIndex is simply supernodeSorterPortal.Get(supernode);
-    // copy the global regular Id and data value
-    // Set: augmentedTree->regularNodeGlobalIDs[newRegularID] = globalRegularIDSet[supernodeSetIndex];
-    auto permutedAugmentedTreeRegularNodeGlobalIds = vtkm::cont::make_ArrayHandlePermutation(
-      newRegularId, this->AugmentedTree->RegularNodeGlobalIds);
+    // create fancy arrays needed to allow use of FieldIn for worklet parameters
+    auto permutedSupernodeIdSet =
+      vtkm::cont::make_ArrayHandlePermutation(this->SupernodeSorter, this->SupernodeIdSet);
+    auto permutedBaseTreeSuper2Hypernode = vtkm::cont::make_ArrayHandlePermutation(
+      permutedSupernodeIdSet,
+      this->BaseTree
+        ->Super2Hypernode); // Get this->BaseTree->Super2hypernode[oldSupernodeId]; as FieldIn for the worklet
     auto permutedGlobalRegularIdSet =
       vtkm::cont::make_ArrayHandlePermutation(this->SupernodeSorter, this->GlobalRegularIdSet);
-    vtkm::cont::Algorithm::Copy(permutedGlobalRegularIdSet,
-                                permutedAugmentedTreeRegularNodeGlobalIds);
-    // SetL augmentedTree->dataValues[newRegularID] = dataValueSet[supernodeSetIndex];
-    auto permutedAugmentedTreeDataValues =
-      vtkm::cont::make_ArrayHandlePermutation(newRegularId, this->AugmentedTree->DataValues);
     auto permutedDataValueSet =
       vtkm::cont::make_ArrayHandlePermutation(this->SupernodeSorter, this->DataValueSet);
-    vtkm::cont::Algorithm::Copy(permutedDataValueSet, permutedAugmentedTreeDataValues);
+    auto oldRegularId =
+      vtkm::cont::make_ArrayHandlePermutation(permutedSupernodeIdSet, this->BaseTree->Supernodes);
+    auto oldSuperFrom =
+      vtkm::cont::make_ArrayHandlePermutation(oldRegularId, this->BaseTree->Superparents);
 
-    // the sort order will be dealt with later
-    // since all of these nodes are supernodes, they will be their own superparent, which means that:
-    //  a.  the regular2node can be set immediately
-    //      Set: augmentedTree->regular2supernode[newRegularID] = newSupernodeID;
-    auto permutedAugmentedTreeRegular2Supernode =
-      vtkm::cont::make_ArrayHandlePermutation(newRegularId, this->AugmentedTree->Regular2Supernode);
-    vtkm::cont::Algorithm::Copy(newSupernodeId, permutedAugmentedTreeRegular2Supernode);
+    // Create a view of the range of this->AugmentedTree->Superarcs that will be updated by the worklet
+    // so that we can use FieldOut as the type in the worklet
+    auto augmentedTreeSuperarcsView = vtkm::cont::make_ArrayHandleView(
+      this->AugmentedTree->Superarcs,
+      numSupernodesAlready,                     // start view here
+      this->SupernodeSorter.GetNumberOfValues() // select this many values
+    );
+    auto augmentedTreeHyperparentsView = vtkm::cont::make_ArrayHandleView(
+      this->AugmentedTree->Hyperparents,
+      numSupernodesAlready,                     // start view here
+      this->SupernodeSorter.GetNumberOfValues() // select this many values
+    );
+    auto augmentedTreeSuper2HypernodeView = vtkm::cont::make_ArrayHandleView(
+      this->AugmentedTree->Super2Hypernode,
+      numSupernodesAlready,                     // start view here
+      this->SupernodeSorter.GetNumberOfValues() // select this many values
+    );
+    auto augmentedTreeWhichRoundView = vtkm::cont::make_ArrayHandleView(
+      this->AugmentedTree->WhichRound,
+      numSupernodesAlready,                     // start view here
+      this->SupernodeSorter.GetNumberOfValues() // select this many values
+    );
+    auto augmentedTreeWhichIterationView = vtkm::cont::make_ArrayHandleView(
+      this->AugmentedTree->WhichIteration,
+      numSupernodesAlready,                     // start view here
+      this->SupernodeSorter.GetNumberOfValues() // select this many values
+    );
+    auto augmentedTreeRegularNodeGlobalIdsView = vtkm::cont::make_ArrayHandleView(
+      this->AugmentedTree->RegularNodeGlobalIds,
+      numSupernodesAlready,                     // start view here
+      this->SupernodeSorter.GetNumberOfValues() // select this many values
+    );
+    auto augmentedTreeDataValuesView = vtkm::cont::make_ArrayHandleView(
+      this->AugmentedTree->DataValues,
+      numSupernodesAlready,                     // start view here
+      this->SupernodeSorter.GetNumberOfValues() // select this many values
+    );
+    auto augmentedTreeRegular2SupernodeView = vtkm::cont::make_ArrayHandleView(
+      this->AugmentedTree->Regular2Supernode,
+      numSupernodesAlready,                     // start view here
+      this->SupernodeSorter.GetNumberOfValues() // select this many values
+    );
+    auto augmentedTreeSuperparentsView = vtkm::cont::make_ArrayHandleView(
+      this->AugmentedTree->Superparents,
+      numSupernodesAlready,                     // start view here
+      this->SupernodeSorter.GetNumberOfValues() // select this many values
+    );
 
-    //  b.  as can the superparent
-    //      Set: augmentedTree->superparents[newRegularID] = newSupernodeID;
-    auto permutedAugmentedTreeSuperparents =
-      vtkm::cont::make_ArrayHandlePermutation(newRegularId, this->AugmentedTree->Superparents);
-    vtkm::cont::Algorithm::Copy(newSupernodeId, permutedAugmentedTreeSuperparents);
-  } // END scope for e. to delete temporary variables
+    // Required execution objects to call other functions
+    auto findSuperArcForUnknownNode = this->AugmentedTree->GetFindSuperArcForUnknownNode();
 
-  // We have one last bit of cleanup to do.  If there were attachment points, then the round in which they transfer has been removed
-  // While it is possible to turn this into a null round, it is better to reduce the iteration count by one and resize the arrays
-  // To do this, we access the *LAST* element written and check to see whether it is in the final iteration (according to the base tree)
+    // Execution object used to encapsulate data form the BaseTree to avoid the limit of 20 input parameters per worklet
+    auto createSuperarcsDataExecObj =
+      vtkm::worklet::contourtree_distributed::hierarchical_augmenter::CreateSuperarcsDataExec(
+        this->BaseTree->Hyperparents,
+        this->BaseTree->WhichRound,
+        this->BaseTree->WhichIteration,
+        this->BaseTree->Superarcs,
+        this->BaseTree->Hypernodes,
+        this->SuperparentSet,
+        this->NewSupernodeIds);
+
+    // invoke the worklet
+    this->Invoke(createSuperarcsWorklet, // the worklet
+                 // inputs
+                 this->SupernodeSorter,           // input domain (WholeArrayIn)
+                 permutedSupernodeIdSet,          // input (FieldIn)
+                 permutedBaseTreeSuper2Hypernode, // input (FieldIn)
+                 permutedGlobalRegularIdSet,      // input (FieldIn)
+                 permutedDataValueSet,            // input (FieldIn)
+                 oldSuperFrom,                    // input (FieldIn)
+                 findSuperArcForUnknownNode,      // input (Execution object)
+                 createSuperarcsDataExecObj,      // input (Execution object with BaseTreeData
+                 // Outputs
+                 this->AugmentedTree->Supernodes,       // input/output
+                 augmentedTreeSuperarcsView,            // output
+                 augmentedTreeHyperparentsView,         // output
+                 augmentedTreeSuper2HypernodeView,      // output
+                 augmentedTreeWhichRoundView,           // output
+                 augmentedTreeWhichIterationView,       // output
+                 augmentedTreeRegularNodeGlobalIdsView, // output
+                 augmentedTreeDataValuesView,           // output
+                 augmentedTreeRegular2SupernodeView,    // output
+                 augmentedTreeSuperparentsView          // output
+    );
+
+  } // END Call CreateSuperarcsWorklet
+
+#ifdef DEBUG_PRINT
+  VTKM_LOG_S(vtkm::cont::LogLevel::Info,
+             DebugPrint(std::string("Round ") + std::to_string(roundNumber) +
+                          std::string(" Details Filled in For Supernodes "),
+                        __FILE__,
+                        __LINE__));
+#endif
+
+  // Now, in order to set the first supernode per iteration, we do an additional loop
+  // We are guaranteed that all supernodes at this level are implicitly sorted by iteration, so we test for ends of segments
+  // NOTE that we do this after the previous loop, since we depend on a value that it has set
+  {
+    vtkm::cont::ArrayHandleCounting<vtkm::Id> tempIndex(1, 1, currNumIterations - 1);
+    vtkm::worklet::contourtree_distributed::hierarchical_augmenter::
+      CreateSuperarcsSetFirstSupernodePerIterationWorklet
+        createSuperarcsSetFirstSupernodePerIterationWorklet(numSupernodesAlready);
+    vtkm::cont::ArrayHandleIndex tempSupernodeIndex(this->SupernodeSorter.GetNumberOfValues());
+    this->Invoke(createSuperarcsSetFirstSupernodePerIterationWorklet,
+                 tempSupernodeIndex,
+                 this->AugmentedTree->WhichIteration,                         // input
+                 this->AugmentedTree->FirstSupernodePerIteration[roundNumber] // input/output
+    );
+  }
+
+  // since there's an extra entry in the firstSupernode array as a sentinel, set it
+  vtkm::worklet::contourtree_augmented::IdArraySetValue(
+    currNumIterations,                                           // index
+    this->AugmentedTree->Supernodes.GetNumberOfValues(),         // new value
+    this->AugmentedTree->FirstSupernodePerIteration[roundNumber] // array
+  );
+
+  // The following loop should be safe in parallel since there should never be two zeros in sequence, i.e., the next
+  // entry after a zero will always be valid, regardless of execution order
+  // This was added because in rare cases there are no supernodes transferred in an iteration, for example because there
+  // are no available upper leaves to prune. If this is case, we are guaranteed that there will be available lower leaves
+  // so the next iteration will have a non-zero number.  We had a major bug from this, and it's cropped back up in the
+  // Hierarchical Augmentation, so I'm expanding the comment just in case.
+  {
+    vtkm::cont::ArrayHandleCounting<vtkm::Id> tempIndex(1, 1, currNumIterations - 1);
+    vtkm::worklet::contourtree_distributed::hierarchical_augmenter::
+      CreateSuperarcsUpdateFirstSupernodePerIterationWorklet
+        createSuperarcsUpdateFirstSupernodePerIterationWorklet;
+    this->Invoke(createSuperarcsUpdateFirstSupernodePerIterationWorklet,
+                 tempIndex,                                                   // input index
+                 this->AugmentedTree->FirstSupernodePerIteration[roundNumber] // input/output
+    );
+  }
+
+  // We have one last bit of cleanup to do.  If there were attachment points,
+  // then the round in which they transfer has been removed
+  // While it is possible to turn this into a null round, it is better to
+  // reduce the iteration count by one and resize the arrays
+  // To do this, we access the *LAST* element written and check to see whether
+  // it is in the final iteration (according to the base tree)
   // But there might be *NO* supernodes in the round, so we check first
-  vtkm::Id iterationArraySize =
-    vtkm::cont::ArrayGetValue(roundNumber, this->AugmentedTree->NumIterations);
-  if (iterationArraySize > 0)
+  if (currNumIterations > 0)
   { // at least one iteration
     vtkm::Id lastSupernodeThisLevel = this->AugmentedTree->Supernodes.GetNumberOfValues() - 1;
     vtkm::Id lastIterationThisLevel = vtkm::worklet::contourtree_augmented::MaskedIndex(
       vtkm::cont::ArrayGetValue(lastSupernodeThisLevel, this->AugmentedTree->WhichIteration));
-    // if there were no attachment points, it will be in the last iteration: if there were attachment points, it will be in the previous one
-    if (lastIterationThisLevel < iterationArraySize - 1)
+    // if there were no attachment points, it will be in the last iteration: if there were
+    // attachment points, it will be in the previous one
+    if (lastIterationThisLevel < currNumIterations - 1)
     { // attachment point round was removed
       // decrement the iteration count (still with an extra element as sentinel)
+      vtkm::Id iterationArraySize = currNumIterations;
+      // decrease iterations by 1. I.e,: augmentedTree->nIterations[roundNo]--;
       vtkm::worklet::contourtree_augmented::IdArraySetValue(
-        roundNumber, iterationArraySize - 1, this->AugmentedTree->NumIterations);
+        roundNumber,                       // index
+        currNumIterations - 1,             // new value
+        this->AugmentedTree->NumIterations // array
+      );
       // shrink the supernode array
       this->AugmentedTree->FirstSupernodePerIteration[roundNumber].Allocate(
         iterationArraySize, vtkm::CopyFlag::On); // shrink array but keep values
       vtkm::worklet::contourtree_augmented::IdArraySetValue(
-        iterationArraySize - 1,
-        this->AugmentedTree->Supernodes.GetNumberOfValues(),
-        this->AugmentedTree->FirstSupernodePerIteration[roundNumber]);
+        iterationArraySize - 1,                                      // index
+        this->AugmentedTree->Supernodes.GetNumberOfValues(),         // new value
+        this->AugmentedTree->FirstSupernodePerIteration[roundNumber] // array
+      );
       // for the hypernode array, the last iteration is guaranteed not to have hyperarcs by construction
       // so the last iteration will already have the correct sentinel value, and we just need to shrink the array
-      this->AugmentedTree->FirstHypernodePerIteration[roundNumber].Allocate(
+      this->AugmentedTree->FirstSupernodePerIteration[roundNumber].Allocate(
         iterationArraySize, vtkm::CopyFlag::On); // shrink array but keep values
     }                                            // attachment point round was removed
   }                                              // at least one iteration
 
+
+#ifdef DEBUG_PRINT
+  VTKM_LOG_S(vtkm::cont::LogLevel::Info,
+             DebugPrint(std::string("Round ") + std::to_string(roundNumber) +
+                          std::string(" Superarcs Created "),
+                        __FILE__,
+                        __LINE__));
+#endif
+
   // in the interests of debug, we resize the sorting array to zero here,
   // even though we will re-resize them in the next function
   this->SupernodeSorter.ReleaseResources();

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

@@ -1054,6 +1054,13 @@ void HierarchicalContourTree<FieldType>::AddToVTKMDataSet(vtkm::cont::DataSet& d
   ds.AddField(firstSupernodePerIterationOffsetsField);
   // TODO/FIXME: It seems we may only need the counts for the first iteration, so check, which
   // information we actually need.
+  // Add the number of rounds as an array of length 1
+  vtkm::cont::ArrayHandle<vtkm::Id> tempNumRounds;
+  tempNumRounds.Allocate(1);
+  vtkm::worklet::contourtree_augmented::IdArraySetValue(0, this->NumRounds, tempNumRounds);
+  vtkm::cont::Field numRoundsField(
+    "NumRounds", vtkm::cont::Field::Association::WholeDataSet, tempNumRounds);
+  ds.AddField(numRoundsField);
 }
 
 } // namespace contourtree_distributed

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

@@ -91,7 +91,9 @@
 #include <vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_hyper_sweeper/InitializeIntrinsicVertexCountSubtractLowEndWorklet.h>
 #include <vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_hyper_sweeper/TransferTargetComperator.h>
 #include <vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_hyper_sweeper/TransferWeightsUpdateLHEWorklet.h>
+#include <vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_hyper_sweeper/TransferWeightsUpdateLHEWorkletRound2.h>
 #include <vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_hyper_sweeper/TransferWeightsUpdateRHEWorklet.h>
+#include <vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_hyper_sweeper/TransferWeightsUpdateRHEWorkletRound2.h>
 
 
 namespace vtkm
@@ -530,6 +532,9 @@ void HierarchicalHyperSweeper<SweepValueType, ContourTreeFieldType>::ComputeSupe
   } // scope ComputeSuperarcTransferWeightsWorklet
 
   // 5. Now we need to sort the transfer targets into contiguous segments
+  // TODO / WARNING 11/07/2023
+  // We have now got a flag of ATTACHMENT_POINT_TRANSFER whose effect is to separate out transfers to
+  // the superarc from transfers to the supernode
   {
     // create view of superSortPermute[firstSupernode, lastSupernode) for sorting
     vtkm::cont::ArrayHandleView<vtkm::worklet::contourtree_augmented::IdArrayType>
@@ -578,6 +583,16 @@ void HierarchicalHyperSweeper<SweepValueType, ContourTreeFieldType>::TransferWei
   vtkm::Id firstSupernode,
   vtkm::Id lastSupernode)
 { // TransferWeights()
+  // TODO / WARNING 11/07/2023
+  // This code was originally written on the assumption that the hierarchical tree had been augmented by the attachment points.
+  // As a result, it assumed that no attachment points remained.
+  // It is now being used for partially augmented versions due to pre-simplification, for which the correct treatment is to
+  // transfer not as dependent weight, but as intrinsic weight. Note that this ONLY applies to attachment points: if the
+  // subtree attaches at a proper supernode in the ancestor level, it should still be treated as dependent weight. The logic
+  // behind this is that an attachment point is regular with respect to the superarc along which it inserts.  Adding it as
+  // dependent weight means that it is treated as *OUTSIDE* the superarc in the reverse sweep (or equivalent computation)
+  // Treating it as dependent weight means that both ends of the superarc end up with the correct value.
+
   // 7. Now perform a segmented prefix sum
   vtkm::Id numSupernodesToProcess = lastSupernode - firstSupernode;
   // Same as std::partial_sum(valuePrefixSum.begin() + firstSupernode, valuePrefixSum.begin() + lastSupernode, valuePrefixSum.begin() + firstSupernode);
@@ -602,6 +617,12 @@ void HierarchicalHyperSweeper<SweepValueType, ContourTreeFieldType>::TransferWei
 
   // 7a. and 7b.
   {
+    // TODO / WARNING 11/07/2023
+    // Before dealing with attachment points, we just transferred by segments.  We now have
+    // the possibility of transferring some weight at an attachment point,
+    // and some not.  To avoid write conflicts, we treat this as two passes: one for attachment
+    // points, one for all others. This means duplicating 7a/7b, sadly.
+
     // 7a. Find the RHE of each group and transfer the prefix sum weight
     // Note that we do not compute the transfer weight separately, we add it in place instead
     // Instantiate the worklet
@@ -634,7 +655,7 @@ void HierarchicalHyperSweeper<SweepValueType, ContourTreeFieldType>::TransferWei
     auto valuePrefixSumPreviousValueView = vtkm::cont::make_ArrayHandleView(
       this->ValuePrefixSum, firstSupernode, numSupernodesToProcess - 1);
 
-    // 7b. Now find the LHE of each group and subtract out the prior weight
+    // 7b (non-attachment). Now find the LHE of each group and subtract out the prior weight.
     vtkm::worklet::contourtree_distributed::hierarchical_hyper_sweeper::
       TransferWeightsUpdateLHEWorklet transferWeightsUpdateLHEWorklet;
     this->Invoke(transferWeightsUpdateLHEWorklet,
@@ -643,6 +664,43 @@ void HierarchicalHyperSweeper<SweepValueType, ContourTreeFieldType>::TransferWei
                  valuePrefixSumPreviousValueView,
                  this->DependentValues);
   }
+
+  // 7a (attachment). Find the RHE of each group and transfer the prefix sum weight
+  // Note that we do not compute the transfer weight separately, we add it in place instead
+  {
+    auto supernodeIndex =
+      vtkm::cont::make_ArrayHandleCounting(firstSupernode, vtkm::Id{ 1 }, numSupernodesToProcess);
+    auto valuePrefixSumView = vtkm::cont::make_ArrayHandleView(
+      this->ValuePrefixSum, firstSupernode, numSupernodesToProcess);
+
+    vtkm::worklet::contourtree_distributed::hierarchical_hyper_sweeper::
+      TransferWeightsUpdateRHEWorkletRound2 transferWeightsUpdateRHEWorkletRound2(lastSupernode);
+    // Invoke the worklet
+    this->Invoke(transferWeightsUpdateRHEWorkletRound2, // worklet
+                 supernodeIndex, // input counting array [firstSupernode, lastSupernode)
+                 this->SortedTransferTarget,
+                 valuePrefixSumView, // input view of valuePrefixSum[firstSupernode, lastSupernode)
+                 this->IntrinsicValues,
+                 this->DependentValues);
+  }
+  // 7b (i). Now find the LHE of each group and subtract out the prior weight.
+  {
+    auto sortedTransferTargetView = vtkm::cont::make_ArrayHandleView(
+      this->SortedTransferTarget, firstSupernode + 1, numSupernodesToProcess - 1);
+    auto sortedTransferTargetShiftedView = vtkm::cont::make_ArrayHandleView(
+      this->SortedTransferTarget, firstSupernode, numSupernodesToProcess - 1);
+    auto valuePrefixSumPreviousValueView = vtkm::cont::make_ArrayHandleView(
+      this->ValuePrefixSum, firstSupernode, numSupernodesToProcess - 1);
+
+    vtkm::worklet::contourtree_distributed::hierarchical_hyper_sweeper::
+      TransferWeightsUpdateLHEWorkletRound2 transferWeightsUpdateLHEWorkletRound2;
+    this->Invoke(transferWeightsUpdateLHEWorkletRound2,
+                 sortedTransferTargetView,
+                 sortedTransferTargetShiftedView,
+                 valuePrefixSumPreviousValueView,
+                 this->IntrinsicValues,
+                 this->DependentValues);
+  }
 } // TransferWeights()
 
 

vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_augmenter/CMakeLists.txt

@@ -23,6 +23,9 @@ set(headers
   AttachmentAndSupernodeComparator.h
   ResizeArraysBuildNewSupernodeIdsWorklet.h
   CreateSuperarcsWorklet.h
+  CreateSuperarcsData.h
+  CreateSuperarcsSetFirstSupernodePerIterationWorklet.h
+  CreateSuperarcsUpdateFirstSupernodePerIterationWorklet.h
   HierarchicalAugmenterInOutData.h
 )
 

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

@@ -0,0 +1,169 @@
+//============================================================================
+//  Copyright (c) Kitware, Inc.
+//  All rights reserved.
+//  See LICENSE.txt for details.
+//
+//  This software is distributed WITHOUT ANY WARRANTY; without even
+//  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+//  PURPOSE.  See the above copyright notice for more information.
+//============================================================================
+// Copyright (c) 2018, The Regents of the University of California, through
+// Lawrence Berkeley National Laboratory (subject to receipt of any required approvals
+// from the U.S. Dept. of Energy).  All rights reserved.
+//
+// Redistribution and use in source and binary forms, with or without modification,
+// are permitted provided that the following conditions are met:
+//
+// (1) Redistributions of source code must retain the above copyright notice, this
+//     list of conditions and the following disclaimer.
+//
+// (2) Redistributions in binary form must reproduce the above copyright notice,
+//     this list of conditions and the following disclaimer in the documentation
+//     and/or other materials provided with the distribution.
+//
+// (3) Neither the name of the University of California, Lawrence Berkeley National
+//     Laboratory, U.S. Dept. of Energy nor the names of its contributors may be
+//     used to endorse or promote products derived from this software without
+//     specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+// ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+// WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
+// IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
+// INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
+// BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+// DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
+// OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
+// OF THE POSSIBILITY OF SUCH DAMAGE.
+//
+//=============================================================================
+//
+//  This code is an extension of the algorithm presented in the paper:
+//  Parallel Peak Pruning for Scalable SMP Contour Tree Computation.
+//  Hamish Carr, Gunther Weber, Christopher Sewell, and James Ahrens.
+//  Proceedings of the IEEE Symposium on Large Data Analysis and Visualization
+//  (LDAV), October 2016, Baltimore, Maryland.
+//
+//  The PPP2 algorithm and software were jointly developed by
+//  Hamish Carr (University of Leeds), Gunther H. Weber (LBNL), and
+//  Oliver Ruebel (LBNL)
+//==============================================================================
+
+// This header contains an Execution Object used to pass a arrays to the
+// CreateSuperarcsWorklet to overcome the limitation of 20 input parameters for a worklet
+
+#ifndef vtk_m_worklet_contourtree_distributed_hierarchical_augmenter_create_superarcs_data_h
+#define vtk_m_worklet_contourtree_distributed_hierarchical_augmenter_create_superarcs_data_h
+
+#include <vtkm/filter/scalar_topology/worklet/contourtree_augmented/Types.h>
+
+#include <vtkm/cont/ExecutionObjectBase.h>
+
+namespace vtkm
+{
+namespace worklet
+{
+namespace contourtree_distributed
+{
+namespace hierarchical_augmenter
+{
+
+
+class CreateSuperarcsData
+{
+public:
+  // Sort indicies types
+  using IndicesPortalType = vtkm::worklet::contourtree_augmented::IdArrayType::ReadPortalType;
+
+  VTKM_EXEC_CONT
+  CreateSuperarcsData() {}
+
+  VTKM_CONT
+  CreateSuperarcsData(
+    const vtkm::worklet::contourtree_augmented::IdArrayType& baseTreeHyperparents,
+    const vtkm::worklet::contourtree_augmented::IdArrayType& baseTreeWhichRound,
+    const vtkm::worklet::contourtree_augmented::IdArrayType& baseTreeWhichIteration,
+    const vtkm::worklet::contourtree_augmented::IdArrayType& baseTreeSuperarcs,
+    const vtkm::worklet::contourtree_augmented::IdArrayType& baseTreeHypernodes,
+    const vtkm::worklet::contourtree_augmented::IdArrayType& superparentSet,
+    const vtkm::worklet::contourtree_augmented::IdArrayType& newSupernodeIds,
+    vtkm::cont::DeviceAdapterId device,
+    vtkm::cont::Token& token)
+  {
+    this->BaseTreeHyperparents = baseTreeHyperparents.PrepareForInput(device, token);
+    this->BaseTreeWhichRound = baseTreeWhichRound.PrepareForInput(device, token);
+    this->BaseTreeWhichIteration = baseTreeWhichIteration.PrepareForInput(device, token);
+    this->BaseTreeSuperarcs = baseTreeSuperarcs.PrepareForInput(device, token);
+    this->BaseTreeHypernodes = baseTreeHypernodes.PrepareForInput(device, token);
+    this->SuperparentSet = superparentSet.PrepareForInput(device, token);
+    this->NewSupernodeIds = newSupernodeIds.PrepareForInput(device, token);
+  }
+
+public:
+  IndicesPortalType BaseTreeHyperparents;
+  IndicesPortalType BaseTreeWhichRound;
+  IndicesPortalType BaseTreeWhichIteration;
+  IndicesPortalType BaseTreeSuperarcs;
+  IndicesPortalType BaseTreeHypernodes;
+  IndicesPortalType SuperparentSet;
+  IndicesPortalType NewSupernodeIds;
+};
+
+
+class CreateSuperarcsDataExec : public vtkm::cont::ExecutionObjectBase
+{
+public:
+  VTKM_EXEC_CONT
+  CreateSuperarcsDataExec(
+    const vtkm::worklet::contourtree_augmented::IdArrayType& baseTreeHyperparents,
+    const vtkm::worklet::contourtree_augmented::IdArrayType& baseTreeWhichRound,
+    const vtkm::worklet::contourtree_augmented::IdArrayType& baseTreeWhichIteration,
+    const vtkm::worklet::contourtree_augmented::IdArrayType& baseTreeSuperarcs,
+    const vtkm::worklet::contourtree_augmented::IdArrayType& baseTreeHypernodes,
+    const vtkm::worklet::contourtree_augmented::IdArrayType& superparentSet,
+    const vtkm::worklet::contourtree_augmented::IdArrayType& newSupernodeIds)
+    : BaseTreeHyperparents(baseTreeHyperparents)
+    , BaseTreeWhichRound(baseTreeWhichRound)
+    , BaseTreeWhichIteration(baseTreeWhichIteration)
+    , BaseTreeSuperarcs(baseTreeSuperarcs)
+    , BaseTreeHypernodes(baseTreeHypernodes)
+    , SuperparentSet(superparentSet)
+    , NewSupernodeIds(newSupernodeIds)
+  {
+  }
+
+  VTKM_CONT
+  CreateSuperarcsData PrepareForExecution(vtkm::cont::DeviceAdapterId device,
+                                          vtkm::cont::Token& token) const
+  {
+    return CreateSuperarcsData(BaseTreeHyperparents,
+                               BaseTreeWhichRound,
+                               BaseTreeWhichIteration,
+                               BaseTreeSuperarcs,
+                               BaseTreeHypernodes,
+                               SuperparentSet,
+                               NewSupernodeIds,
+                               device,
+                               token);
+  }
+
+private:
+  // Whole array data used from the BaseTree in CreateSuperarcsWorklet
+  const vtkm::worklet::contourtree_augmented::IdArrayType& BaseTreeHyperparents;
+  const vtkm::worklet::contourtree_augmented::IdArrayType& BaseTreeWhichRound;
+  const vtkm::worklet::contourtree_augmented::IdArrayType& BaseTreeWhichIteration;
+  const vtkm::worklet::contourtree_augmented::IdArrayType& BaseTreeSuperarcs;
+  const vtkm::worklet::contourtree_augmented::IdArrayType& BaseTreeHypernodes;
+  const vtkm::worklet::contourtree_augmented::IdArrayType& SuperparentSet;
+  const vtkm::worklet::contourtree_augmented::IdArrayType& NewSupernodeIds;
+};
+
+
+
+} // namespace hierarchical_augmenter
+} // namespace contourtree_distributed
+} // namespace worklet
+} // namespace vtkm
+
+#endif

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

@@ -0,0 +1,148 @@
+//============================================================================
+//  Copyright (c) Kitware, Inc.
+//  All rights reserved.
+//  See LICENSE.txt for details.
+//
+//  This software is distributed WITHOUT ANY WARRANTY; without even
+//  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+//  PURPOSE.  See the above copyright notice for more information.
+//============================================================================
+// Copyright (c) 2018, The Regents of the University of California, through
+// Lawrence Berkeley National Laboratory (subject to receipt of any required approvals
+// from the U.S. Dept. of Energy).  All rights reserved.
+//
+// Redistribution and use in source and binary forms, with or without modification,
+// are permitted provided that the following conditions are met:
+//
+// (1) Redistributions of source code must retain the above copyright notice, this
+//     list of conditions and the following disclaimer.
+//
+// (2) Redistributions in binary form must reproduce the above copyright notice,
+//     this list of conditions and the following disclaimer in the documentation
+//     and/or other materials provided with the distribution.
+//
+// (3) Neither the name of the University of California, Lawrence Berkeley National
+//     Laboratory, U.S. Dept. of Energy nor the names of its contributors may be
+//     used to endorse or promote products derived from this software without
+//     specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+// ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+// WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
+// IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
+// INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
+// BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+// DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
+// OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
+// OF THE POSSIBILITY OF SUCH DAMAGE.
+//
+//=============================================================================
+//  The PPP2 algorithm and software were jointly developed by
+//  Hamish Carr (University of Leeds), Gunther H. Weber (LBNL), and
+//  Oliver Ruebel (LBNL)
+//==============================================================================
+
+#ifndef vtk_m_worklet_contourtree_distributed_hierarchical_augmenter_create_auperarcs_set_first_supernode_per_iteration_worklet_h
+#define vtk_m_worklet_contourtree_distributed_hierarchical_augmenter_create_auperarcs_set_first_supernode_per_iteration_worklet_h
+
+#include <vtkm/filter/scalar_topology/worklet/contourtree_augmented/Types.h>
+#include <vtkm/worklet/WorkletMapField.h>
+
+namespace vtkm
+{
+namespace worklet
+{
+namespace contourtree_distributed
+{
+namespace hierarchical_augmenter
+{
+
+/// Worklet used in HierarchicalAugmenter::UpdateHyperstructure to set the hyperarcs and hypernodes
+class CreateSuperarcsSetFirstSupernodePerIterationWorklet : public vtkm::worklet::WorkletMapField
+{
+public:
+  /// Control signature for the worklet
+  using ControlSignature = void(FieldIn supernodeIndex,
+                                WholeArrayIn augmentedTreeWhichIteration,
+                                WholeArrayInOut augmentedTreeFirstSupernodePerIteration);
+  using ExecutionSignature = void(_1, _2, _3);
+  using InputDomain = _1;
+
+  // Default Constructor
+  VTKM_EXEC_CONT
+  CreateSuperarcsSetFirstSupernodePerIterationWorklet(vtkm::Id numSupernodesAlready)
+    : NumSupernodesAlready(numSupernodesAlready)
+  {
+  }
+
+  template <typename InFieldPortalType, typename InOutFieldPortalType>
+  VTKM_EXEC void operator()(
+    const vtkm::Id& supernode, // index in supernodeSorter
+    // const vtkm::Id& supernodeSetindex,  // supernodeSorter[supernode]
+    const InFieldPortalType& augmentedTreeWhichIterationPortal,
+    const InOutFieldPortalType& augmentedTreeFirstSupernodePerIterationPortal) const
+  { // operator()()
+    // per supernode in the set
+    // retrieve the index from the sorting index array (Done on input)(NOT USED)
+    // indexType supernodeSetIndex = supernodeSorter[supernode];
+
+    // work out the new supernode ID
+    vtkm::Id newSupernodeId = this->NumSupernodesAlready + supernode;
+
+    // The 0th element sets the first element in the zeroth iteration
+    if (supernode == 0)
+    {
+      augmentedTreeFirstSupernodePerIterationPortal.Set(0, newSupernodeId);
+    }
+    // otherwise, mismatch to the left identifies a new iteration
+    else
+    {
+      if (vtkm::worklet::contourtree_augmented::MaskedIndex(
+            augmentedTreeWhichIterationPortal.Get(newSupernodeId)) !=
+          vtkm::worklet::contourtree_augmented::MaskedIndex(
+            augmentedTreeWhichIterationPortal.Get(newSupernodeId - 1)))
+      { // start of segment
+        augmentedTreeFirstSupernodePerIterationPortal.Set(
+          vtkm::worklet::contourtree_augmented::MaskedIndex(
+            augmentedTreeWhichIterationPortal.Get(newSupernodeId)),
+          newSupernodeId);
+      } // start of segmen
+    }
+
+    /*
+    #pragma omp parallel for
+    for (indexType supernode = 0; supernode < supernodeSorter.size(); supernode++)
+    { // per supernode in the set
+      // retrieve the index from the sorting index array
+      indexType supernodeSetIndex = supernodeSorter[supernode];
+
+      // work out the new supernode ID
+      indexType newSupernodeID = nSupernodesAlready + supernode;
+
+      // The 0th element sets the first element in the zeroth iteration
+      if (supernode == 0)
+        augmentedTree->firstSupernodePerIteration[roundNo][0] = newSupernodeID;
+      // otherwise, mismatch to the left identifies a new iteration
+      else
+      {
+        if (augmentedTree->whichIteration[newSupernodeID] != augmentedTree->whichIteration[newSupernodeID-1])
+          augmentedTree->firstSupernodePerIteration[roundNo][maskedIndex(augmentedTree->whichIteration[newSupernodeID])] = newSupernodeID;
+      }
+    } // per supernode in the set
+    */
+
+  } // operator()()
+
+private:
+  vtkm::Id NumSupernodesAlready;
+
+
+}; // CreateSuperarcsSetFirstSupernodePerIterationWorklet
+
+} // namespace hierarchical_augmenter
+} // namespace contourtree_distributed
+} // namespace worklet
+} // namespace vtkm
+
+#endif

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

@@ -0,0 +1,104 @@
+//============================================================================
+//  Copyright (c) Kitware, Inc.
+//  All rights reserved.
+//  See LICENSE.txt for details.
+//
+//  This software is distributed WITHOUT ANY WARRANTY; without even
+//  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+//  PURPOSE.  See the above copyright notice for more information.
+//============================================================================
+// Copyright (c) 2018, The Regents of the University of California, through
+// Lawrence Berkeley National Laboratory (subject to receipt of any required approvals
+// from the U.S. Dept. of Energy).  All rights reserved.
+//
+// Redistribution and use in source and binary forms, with or without modification,
+// are permitted provided that the following conditions are met:
+//
+// (1) Redistributions of source code must retain the above copyright notice, this
+//     list of conditions and the following disclaimer.
+//
+// (2) Redistributions in binary form must reproduce the above copyright notice,
+//     this list of conditions and the following disclaimer in the documentation
+//     and/or other materials provided with the distribution.
+//
+// (3) Neither the name of the University of California, Lawrence Berkeley National
+//     Laboratory, U.S. Dept. of Energy nor the names of its contributors may be
+//     used to endorse or promote products derived from this software without
+//     specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+// ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+// WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
+// IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
+// INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
+// BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+// DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
+// OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
+// OF THE POSSIBILITY OF SUCH DAMAGE.
+//
+//=============================================================================
+//  The PPP2 algorithm and software were jointly developed by
+//  Hamish Carr (University of Leeds), Gunther H. Weber (LBNL), and
+//  Oliver Ruebel (LBNL)
+//==============================================================================
+
+#ifndef vtk_m_worklet_contourtree_distributed_hierarchical_augmenter_create_auperarcs_update_first_supernode_per_iteration_worklet_h
+#define vtk_m_worklet_contourtree_distributed_hierarchical_augmenter_create_auperarcs_update_first_supernode_per_iteration_worklet_h
+
+#include <vtkm/filter/scalar_topology/worklet/contourtree_augmented/Types.h>
+#include <vtkm/worklet/WorkletMapField.h>
+
+namespace vtkm
+{
+namespace worklet
+{
+namespace contourtree_distributed
+{
+namespace hierarchical_augmenter
+{
+
+/// Worklet used in HierarchicalAugmenter::UpdateHyperstructure to set the hyperarcs and hypernodes
+class CreateSuperarcsUpdateFirstSupernodePerIterationWorklet : public vtkm::worklet::WorkletMapField
+{
+public:
+  /// Control signature for the worklet
+  using ControlSignature = void(FieldIn indexArray,
+                                WholeArrayInOut augmentedTreeFirstSupernodePerIteration);
+  using ExecutionSignature = void(_1, _2);
+  using InputDomain = _1;
+
+  // Default Constructor
+  VTKM_EXEC_CONT
+  CreateSuperarcsUpdateFirstSupernodePerIterationWorklet() {}
+
+  template <typename InOutFieldPortalType>
+  VTKM_EXEC void operator()(
+    const vtkm::Id& iteration,
+    const InOutFieldPortalType& augmentedTreeFirstSupernodePerIterationPortal) const
+  { // operator()()
+    if (augmentedTreeFirstSupernodePerIterationPortal.Get(iteration) == 0)
+    {
+      augmentedTreeFirstSupernodePerIterationPortal.Set(
+        iteration, augmentedTreeFirstSupernodePerIterationPortal.Get(iteration + 1));
+    }
+
+    /*
+    for (indexType iteration = 1; iteration < augmentedTree->nIterations[roundNo]; ++iteration)
+    {
+      if (augmentedTree->firstSupernodePerIteration[roundNo][iteration] == 0)
+      {
+        augmentedTree->firstSupernodePerIteration[roundNo][iteration] =
+          augmentedTree->firstSupernodePerIteration[roundNo][iteration+1];
+      }
+    }
+    */
+  } // operator()()
+};  // CreateSuperarcsUpdateFirstSupernodePerIterationWorklet
+
+} // namespace hierarchical_augmenter
+} // namespace contourtree_distributed
+} // namespace worklet
+} // namespace vtkm
+
+#endif

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

@@ -61,6 +61,7 @@ namespace hierarchical_augmenter
 /// Worklet used to implement the main part of HierarchicalAugmenter::CreateSuperarcs
 /// Connect superarcs for the level & set hyperparents & superchildren count, whichRound,
 /// whichIteration, super2hypernode
+template <typename FieldType>
 class CreateSuperarcsWorklet : public vtkm::worklet::WorkletMapField
 {
 public:
@@ -69,342 +70,559 @@ public:
   /// Control signature for the worklet
   /// @param[in] supernodeSorter input domain. We need access to InputIndex and InputIndex+1,
   ///                           therefore this is a WholeArrayIn transfer.
-  /// @param[in] superparentSet WholeArrayIn because we need access to superparentSet[supernodeSorter[InputIndex]]
-  ///                           and superparentSet[supernodeSorter[InputIndex+1]].
-  /// @param[in] baseTreeSuperarcs WholeArrayIn because we need access to baseTreeSuperarcsPortal.Get(superparentOldSuperId)
-  ///                           While this could be done with fancy array magic, it would require a sequence of multiple
-  ///                           fancy arrays and would likely not be cheaper then computing things in the worklet.
-  /// @param[in] newSupernodeIds WholeArrayIn because we need to access newSupernodeIdsPortal.Get(oldTargetSuperId)
-  ///                           where oldTargetSuperId is the unmasked baseTreeSuperarcsPortal.Get(superparentOldSuperId)
-  /// @param[in] baseTreeSupernodes WholeArrayIn because we need to access baseTreeSupernodesPortal.Get(superparentOldSuperId);
-  /// @param[in] baseTreeRegularNodeGlobalIds WholeArrayIn because we need to access
-  ///                           baseTreeRegularNodeGlobalIdsPortal.Get(superparentOldSuperId);
-  /// @param[in] globalRegularIdSet FieldInd. Permute globalRegularIdSet with supernodeSorter in order to allow this to be a FieldIn.
-  /// @param[in] baseTreeSuper2Hypernode WholeArrayIn because we need to access
-  ///                           baseTreeSuper2HypernodePortal.Get(superparentOldSuperId)
-  /// @param[in] baseTreeWhichIteration WholeArrayIn because we need to access baseTreeWhichIterationPortal.Get(superparentOldSuperId)
-  ///                           and baseTreeWhichIterationPortal.Get(superparentOldSuperId+1)
-  /// @param[in] augmentedTreeSuperarcsView  output view of  this->AugmentedTree->Superarcs with
+  /// @param[in] supernodeIdSetPermuted Field in of supernodeIdSet permuted by the supernodeSorter array to
+  ///                                 allow us to use FieldIn
+  /// @param[in] baseTreeHyperparents whole array input of the BaseTree->Hyperparents
+  /// @param[in] baseTreeSuper2HypernodePermuted  baseTreeSuper2Hypernode permuted by supernodeIdSetPermuted in order
+  ///                          to extract baseTree->super2hypernode[oldSupernodeID];
+  /// @param[in] baseTreeWhichRound
+  /// @param[in] baseTreeWhichIteration
+  /// @param[in] globalRegularIdSetPermuted Field in of globalRegularIdSet permuted by supernodeSorter array to
+  ///                           allow use of FieldIn
+  /// @param[in] dataValueSetPermuted Field in of dataValyeSet permuted by supernodeSorter array to
+  ///                           allow use of FieldIn
+  /// @param[in] baseTreeSuperarcs  BaseTree->Superarcs
+  /// @param[in] oldSuperFrom Permuted baseTree->SuperParents[  baseTree->Supernodes[ supernodeIdSetPermuted  ]  ]
+  /// @param[in] baseTreeHypernodes
+  /// @param[out, in] augmentedTreeSupernodes  this->AugmentedTree->Supernodes array
+  /// @param[out] augmentedTreeSuperarcsView  output view of  this->AugmentedTree->Superarcs with
   ///                           vtkm::cont::make_ArrayHandleView(this->AugmentedTree->Superarcs,
   ///                           numSupernodesAlready, this->SupernodeSorter.GetNumberOfValues()).
   ///                           By using this view allows us to do this one as a FieldOut and it effectively the
   ///                           same as accessing the array at the newSuppernodeId location.
-  /// @param[in] augmentedTreeFirstSupernodePerIteration WholeArrayInOut because we need to update multiple locations.
-  ///                           In is used to preseve original values. Set to augmentedTree->firstSupernodePerIteration[roundNumber].
-  /// @param[in] augmentedTreeSuper2hypernode FieldOut. Output view of this->AugmentedTree->Super2Hypernode
+  /// @param[out] augmentedTreeHyperparentsView  output view of  this->AugmentedTree->Hyperparents with
+  ///                           vtkm::cont::make_ArrayHandleView(this->AugmentedTree->Hyperparents,
+  ///                           numSupernodesAlready, this->SupernodeSorter.GetNumberOfValues()).
+  ///                           By using this view allows us to do this one as a FieldOut and it effectively the
+  ///                           same as accessing the array at the newSuppernodeId location.
+  /// @param[out] augmentedTreeSuper2HypernodeView  output view of  this->AugmentedTree->Super2Hypernode with
   ///                           vtkm::cont::make_ArrayHandleView(this->AugmentedTree->Super2Hypernode,
   ///                           numSupernodesAlready, this->SupernodeSorter.GetNumberOfValues()).
   ///                           By using this view allows us to do this one as a FieldOut and it effectively the
   ///                           same as accessing the array at the newSuppernodeId location.
+  /// @param[out] augmentedTreeWhichRoundView  output view of  this->AugmentedTree->WhichRound with
+  ///                           vtkm::cont::make_ArrayHandleView(this->AugmentedTree->WhichRound,
+  ///                           numSupernodesAlready, this->SupernodeSorter.GetNumberOfValues()).
+  ///                           By using this view allows us to do this one as a FieldOut and it effectively the
+  ///                           same as accessing the array at the newSuppernodeId location.
+  /// @param[out] augmentedTreeWhichIterationView  output view of  this->AugmentedTree->WhichIteration with
+  ///                           vtkm::cont::make_ArrayHandleView(this->AugmentedTree->WhichIteration,
+  ///                           numSupernodesAlready, this->SupernodeSorter.GetNumberOfValues()).
+  ///                           By using this view allows us to do this one as a FieldOut and it effectively the
+  ///                           same as accessing the array at the newSuppernodeId location.
+  /// @param[out] augmentedTreeRegularNodeGlobalIdsView  output view of  this->AugmentedTree->RegularNodeGlobalIds with
+  ///                           vtkm::cont::make_ArrayHandleView(this->AugmentedTree->RegularNodeGlobalIds,
+  ///                           numSupernodesAlready, this->SupernodeSorter.GetNumberOfValues()).
+  ///                           By using this view allows us to do this one as a FieldOut and it effectively the
+  ///                           same as accessing the array at the newRegularId location.
+  /// @param[out] augmentedTreeDataValuesView  output view of  this->AugmentedTree->DataValues with
+  ///                           vtkm::cont::make_ArrayHandleView(this->AugmentedTree->DataValues,
+  ///                           numSupernodesAlready, this->SupernodeSorter.GetNumberOfValues()).
+  ///                           By using this view allows us to do this one as a FieldOut and it effectively the
+  ///                           same as accessing the array at the newRegularId location.
+  /// @param[out] augmentedTreeRegular2SupernodeView  output view of  this->AugmentedTree->Regular2Supernode with
+  ///                           vtkm::cont::make_ArrayHandleView(this->AugmentedTree->Regular2Supernode,
+  ///                           numSupernodesAlready, this->SupernodeSorter.GetNumberOfValues()).
+  ///                           By using this view allows us to do this one as a FieldOut and it effectively the
+  ///                           same as accessing the array at the newRegularId location.
+  /// @param[out] augmentedTreeSuperparentsView  output view of  this->AugmentedTree->Superparents with
+  ///                           vtkm::cont::make_ArrayHandleView(this->AugmentedTree->Superparents,
+  ///                           numSupernodesAlready, this->SupernodeSorter.GetNumberOfValues()).
+  ///                           By using this view allows us to do this one as a FieldOut and it effectively the
+  ///                           same as accessing the array at the newRegularId location.
+
   using ControlSignature = void(
+    // Inputs
     WholeArrayIn supernodeSorter,
-    WholeArrayIn superparentSet,                             // input
-    WholeArrayIn baseTreeSuperarcs,                          // input
-    WholeArrayIn newSupernodeIds,                            // input
-    WholeArrayIn baseTreeSupernodes,                         // input
-    WholeArrayIn baseTreeRegularNodeGlobalIds,               // input
-    FieldIn globalRegularIdSet,                              // input
-    WholeArrayIn baseTreeSuper2Hypernode,                    // input
-    WholeArrayIn baseTreeWhichIteration,                     // input
-    FieldOut augmentedTreeSuperarcsView,                     // output
-    WholeArrayInOut augmentedTreeFirstSupernodePerIteration, // input/output
-    FieldOut augmentedTreeSuper2hypernode                    // ouput
-  );
-  using ExecutionSignature = void(InputIndex, _1, _2, _3, _4, _5, _6, _7, _8, _9, _10, _11, _12);
+    FieldIn supernodeIdSetPermuted,
+    FieldIn baseTreeSuper2HypernodePermuted,
+    FieldIn globalRegularIdSetPermuted,
+    FieldIn dataValueSetPermuted,
+    FieldIn
+      oldSuperFrom, // baseTree->SuperParents[  baseTree->Supernodes[ supernodeIdSetPermuted  ]  ]
+    ExecObject findSuperArcForUnknownNode,
+    ExecObject createSuperarcsData,
+    // Outputs
+    WholeArrayInOut augmentedTreeSupernodes,
+    FieldOut augmentedTreeSuperarcsView,
+    FieldOut augmentedTreeHyperparensView,
+    FieldOut augmentedTreeSuper2Hypernode,
+    FieldOut augmentedTreeWhichRoundView,
+    FieldOut augmentedTreeWhichIterationView,
+    FieldOut augmentedTreeRegularNodeGlobalIdsView,
+    FieldOut augmentedTreeDataValuesView,
+    FieldOut augmentedTreeRegular2SupernodeView,
+    FieldOut augmentedTreeSuperparentsViews);
+  using ExecutionSignature = void(InputIndex,
+                                  _1,
+                                  _2,
+                                  _3,
+                                  _4,
+                                  _5,
+                                  _6,
+                                  _7,
+                                  _8,
+                                  _9,
+                                  _10,
+                                  _11,
+                                  _12,
+                                  _13,
+                                  _14,
+                                  _15,
+                                  _16,
+                                  _17,
+                                  _18);
+  // using ExecutionSignature = void(InputIndex, _1, _2, _3, _4, _5, _6, _7, _8, _9, _10, _11, _12, _13, _14, _15, _16, _17, _18, _19, _20);
   using InputDomain = _1;
 
   /// Default Constructor
   /// @param[in] numSupernodesAlready Set to vtkm::cont::ArrayGetValue(0, this->AugmentedTree->FirstSupernodePerIteration[roundNumber]);
   /// @param[in] baseTreeNumRounds Set to this->BaseTree->NumRounds
-  /// @param[in] augmentedTreeNumIterations Set to  vtkm::cont::ArrayGetValue(roundNumber, this->AugmentedTree->NumIterations);
-  /// @param[in] roundNumber Set the current round
-  /// @param[in] numAugmentedTreeSupernodes Set to augmentedTreeSupernodes this->AugmentedTree->Supernodes.GetNumberOfValues();
+  /// @param[in] numInsertedSupernodes Set to  numInsertedSupernodes
   VTKM_EXEC_CONT
   CreateSuperarcsWorklet(const vtkm::Id& numSupernodesAlready,
                          const vtkm::Id& baseTreeNumRounds,
-                         const vtkm::Id& augmentedTreeNumIterations,
-                         const vtkm::Id& roundNumber,
-                         const vtkm::Id& numAugmentedTreeSupernodes)
+                         const vtkm::Id& numInsertedSupernodes,
+                         const vtkm::Id& roundNo)
     : NumSupernodesAlready(numSupernodesAlready)
     , BaseTreeNumRounds(baseTreeNumRounds)
-    , AugmentedTreeNumIterations(augmentedTreeNumIterations)
-    , RoundNumber(roundNumber)
-    , NumAugmentedTreeSupernodes(numAugmentedTreeSupernodes)
+    , NumInsertedSupernodes(numInsertedSupernodes)
+    , RoundNo(roundNo)
   {
   }
 
   /// operator() of the workelt
-  template <typename InFieldPortalType, typename InOutFieldPortalType>
+  template <typename InFieldPortalType,
+            typename ExecObjType,
+            typename ExecObjectTypeData,
+            typename InOutFieldPortalType>
+  //typename InOutDataFieldPortalType>
+  //typename InOutFieldPortalType,
+  //typename ExecObjectTypeData,
+  //typename ExecObjType>
   VTKM_EXEC void operator()(
+    // Inputs
     const vtkm::Id& supernode, // InputIndex of supernodeSorter
     const InFieldPortalType& supernodeSorterPortal,
-    const InFieldPortalType& superparentSetPortal,
-    const InFieldPortalType& baseTreeSuperarcsPortal,
-    const InFieldPortalType& newSupernodeIdsPortal,
-    const InFieldPortalType& baseTreeSupernodesPortal,
-    const InFieldPortalType& baseTreeRegularNodeGlobalIdsPortal,
-    const vtkm::Id& globalRegularIdSetValue,
-    const InFieldPortalType& baseTreeSuper2HypernodePortal,
-    const InFieldPortalType& baseTreeWhichIterationPortal,
-    vtkm::Id& augmentedTreeSuperarcsValue, // same as augmentedTree->superarcs[newSupernodeId]
-    const InOutFieldPortalType&
-      augmentedTreeFirstSupernodePerIterationPortal, // augmentedTree->firstSupernodePerIteration[roundNumber]
-    vtkm::Id& augmentedTreeSuper2hypernodeValue) const
+    const vtkm::Id& oldSupernodeId, //  supernodeIdSet[supernodeSorterPortal.Get(supernode)]
+    const vtkm::Id&
+      baseTreeSuper2HypernodeOldSupernodeIdValue, // baseTree->super2hypernode[oldSupernodeID];
+    const vtkm::Id&
+      globalRegularIdSetValue, //  globalRegularIdSet[supernodeSorterPortal.Get(supernode)]];
+    const FieldType& dataValueSetValue, //  dataValueSet[supernodeSorterPortal.Get(supernode)]];
+    const vtkm::Id&
+      oldSuperFromValue, // baseTree->SuperParents[  baseTree->Supernodes[ supernodeIdSetPermuted  ]  ]
+    const ExecObjType&
+      findSuperArcForUnknownNode, // Execution object to call FindSuperArcForUnknownNode
+    const ExecObjectTypeData&
+      createSuperarcsData, // Execution object of collect BaseTree data array
+    // Outputs
+    const InOutFieldPortalType& augmentedTreeSupernodesPortal,
+    vtkm::Id& augmentedTreeSuperarcsValue, // set value for AugmentedTree->Auperarcs[newSupernodeId]
+    vtkm::Id&
+      augmentedTreeHyperparentsValue, // set value for AugmentedTree->Hyperparents[newSupernodeId]
+    vtkm::Id&
+      augmentedTreeSuper2HypernodeValue, // set value for AugmentedTree->Super2Hypernode[newSupernodeId]
+    vtkm::Id& augmentedTreeWhichRoundValue,     // AugmentedTree->WhichRound[newSupernodeId]
+    vtkm::Id& augmentedTreeWhichIterationValue, // AugmentedTree->WhichIteration[newSupernodeId]
+    vtkm::Id&
+      augmentedTreeRegularNodeGlobalIdsValue, // AugmentedTree->RegularNodeGlobalIds[newRegularID]
+    FieldType& augmentedTreeDataValuesValue,  // AugmentedTree->DataValues[newRegularID]
+    vtkm::Id& augmentedTreeRegular2SupernodeValue, // AugmentedTree->Regular2Supernode[newRegularID]
+    vtkm::Id& augmentedTreeSuperparentsValue       // AugmentedTree->Superparents[newRegularID]
+  ) const
   {
     // per supernode in the set
     // retrieve the index from the sorting index array
     vtkm::Id supernodeSetIndex = supernodeSorterPortal.Get(supernode);
 
-    // work out the new supernode Id. We have this defined on the outside as a fancy array handle,
-    // however, using the fancy handle here would not really make a performance differnce and
-    // computing it here is more readable
+    // work out the new supernode ID
     vtkm::Id newSupernodeId = this->NumSupernodesAlready + supernode;
 
-    // NOTE: The newRegularId is no longer needed here since all parts
-    //       that used it in the worklet have been moved outside
-    // vtkm::Id newRegularId = newSupernodeId;
+    // and the old supernode ID
+    //  vtkm::Id oldSupernodeId = supernodeIDSet[supernodeSetIndex];  Extracted on call and provides as input
 
-    // NOTE: This part has been moved out of the worklet and is performed using standard vtkm copy constructs
-    // // setting the supernode's regular Id is now trivial
-    // augmentedTreeSupernodesPortal.Set(newSupernodeId, newRegularId);
+    //	At all levels above 0, we used to keep regular vertices in case they are attachment points.
+    //  After augmentation, we don't need to.
+    //	Instead, at all levels above 0, the regular nodes in each round are identical to the supernodes
+    //	In order to avoid confusion, we will copy the ID into a separate variable
+    vtkm::Id newRegularId = newSupernodeId;
 
-    // retrieve the ascending flag from the superparent
-    vtkm::Id superparentSetVal = superparentSetPortal.Get(supernodeSetIndex);
-    // get the ascending flag from the parent
-    bool superarcAscends = vtkm::worklet::contourtree_augmented::IsAscending(superparentSetVal);
-    // strip the ascending flag from the superparent.
-    vtkm::Id superparentOldSuperId =
-      vtkm::worklet::contourtree_augmented::MaskedIndex(superparentSetVal);
+    // setting the supernode's regular ID is now trivial
+    augmentedTreeSupernodesPortal.Set(newSupernodeId, newRegularId);
 
-    // setting the superarc is done the usual way.  Our sort routine has ended up
-    // with the supernodes arranged in either ascending or descending order
-    // inwards along the parent superarc (as expressed by the superparent Id).
-    // Each superarc except the last in the segment points to the next one:
-    // the last one points to the target of the original superarc.
-    // first test to see if we're the last in the array
-    if (supernode == supernodeSorterPortal.GetNumberOfValues() - 1)
-    { // last in the array
-      // special case for root of entire tree at end of top level
-      if (RoundNumber == this->BaseTreeNumRounds)
+    // retrieve the old superID of the superparent.  This is slightly tricky, as we have four classes of supernodes:
+    // 1. the root of the entire tree
+    // 2. attachment points not being inserted. In this case, the supernode ID is stored in the superparentSet
+    //    array, not the superparent for insertion purposes
+    // 3. attachment points being inserted.  In this case, the superparent is stored in the superparentSet array
+    // 4. "ordinary" supernodes, where the superparent is the same as the supernode ID anyway
+    //
+    // Note that an attachment point gets inserted into a parent superarc.  But the attachment point itself has
+    // a NULL superarc, because it's only a virtual insertion.
+    // This means that such an attachment superarc cannot be the superparent of any other attachment point
+    // It is therefore reasonable to deal with 1. & 2 separately. 3. & 4. then combine together
+
+    // first we test for the root of the tree
+    if ((this->RoundNo == BaseTreeNumRounds) &&
+        (supernode == supernodeSorterPortal.GetNumberOfValues() - 1))
+    { // root of the tree
+      // note that oldSupernodeID is guaranteed not to be NO_SUCH_ELEMENT, as the root is in every tree
+      // set the super arrays
+      augmentedTreeSuperarcsValue = vtkm::worklet::contourtree_augmented::NO_SUCH_ELEMENT;
+      // hyperstructure carries over, so we use the same hyperparent as before
+      augmentedTreeHyperparentsValue = createSuperarcsData.BaseTreeHyperparents.Get(oldSupernodeId);
+      // and set the hypernode ID
+      augmentedTreeSuper2HypernodeValue = baseTreeSuper2HypernodeOldSupernodeIdValue;
+      // and the round and iteration
+      augmentedTreeWhichRoundValue = createSuperarcsData.BaseTreeWhichRound.Get(oldSupernodeId);
+      augmentedTreeWhichIterationValue =
+        createSuperarcsData.BaseTreeWhichIteration.Get(oldSupernodeId);
+      // and set the relevant regular arrays
+      augmentedTreeRegularNodeGlobalIdsValue = globalRegularIdSetValue;
+      augmentedTreeDataValuesValue = dataValueSetValue;
+      // for the root, these always point to itself
+      augmentedTreeRegular2SupernodeValue = newSupernodeId;
+      augmentedTreeSuperparentsValue = newSupernodeId;
+    } // root of the tree
+    // now deal with unsimplified attachment points, which we can identify because they were in the "kept" batch, not the "inserted" batch,
+    // and this is given away by the index into the set of supernodes to be added
+    // and the fact that the superarc is NO_SUCH_ELEMENT
+    else if ((supernodeSetIndex >= this->NumInsertedSupernodes) &&
+             (vtkm::worklet::contourtree_augmented::NoSuchElement(
+               createSuperarcsData.BaseTreeSuperarcs.Get(oldSupernodeId))))
+    { // preserved attachment point
+      // note that oldSupernodeID is guaranteed not to be NO_SUCH_ELEMENT, as the supernode came from this block originally
+      // set the superarc to NO_SUCH_ELEMENT, as before
+      augmentedTreeSuperarcsValue = vtkm::worklet::contourtree_augmented::NO_SUCH_ELEMENT;
+      // hyperstructure carries over, so we use the same hyperparent as before
+      augmentedTreeHyperparentsValue = createSuperarcsData.BaseTreeHyperparents.Get(oldSupernodeId);
+      // attachment points are never hypernodes anyway, so set it directly
+      augmentedTreeSuper2HypernodeValue = vtkm::worklet::contourtree_augmented::NO_SUCH_ELEMENT;
+      // and the round and iteration
+      augmentedTreeWhichRoundValue = createSuperarcsData.BaseTreeWhichRound.Get(oldSupernodeId);
+      augmentedTreeWhichIterationValue =
+        createSuperarcsData.BaseTreeWhichIteration.Get(oldSupernodeId);
+      // and set the relevant regular arrays
+      augmentedTreeRegularNodeGlobalIdsValue = globalRegularIdSetValue;
+      augmentedTreeDataValuesValue = dataValueSetValue;
+      // for a preserved attachment point, this always points to itself
+      augmentedTreeRegular2SupernodeValue = newSupernodeId;
+      // the superparent is the tricky one, as the old one may have been broken up by insertions at a higher level
+
+      // Here, what we need to do is a search in the augmented tree to find which superarc to attach to.  This is necessary
+      // because the old superarc it attached to may have been broken up.
+      // We are guaranteed that there is one, and that it only uses the higher levels of the augmented tree,
+      // so the fact that we are partially constructed doesn't get in the way.  To do this, we need supernodes
+      // known to be in the higher level that are above and below the supernode.
+      // Since the point was an attachment point in the base tree, that means that there is a higher round superarc
+      // it inserts into.  Moreover, the algorithm ALWAYS inserts a supernode at or above its original round, so
+      // we can guarantee that both ends of the parent are in the higher levels.  Which means we only need to work
+      // out which end is higher.
+
+      // oldSuperFrom is provided as input and extracted as FieldIn on call. oldRegularId is not needed here
+      //      indexType oldRegularID = baseTree->supernodes[oldSupernodeID];
+      //      indexType oldSuperFrom = baseTree->superparents[oldRegularID];
+      //      indexType oldSuperTo = baseTree->superarcs[oldSuperFrom];
+      vtkm::Id oldSuperToValue = createSuperarcsData.BaseTreeSuperarcs.Get(oldSuperFromValue);
+
+      // retrieve the ascending flag
+      bool ascendingSuperarc = vtkm::worklet::contourtree_augmented::IsAscending(oldSuperToValue);
+      // and mask out the flags
+      vtkm::Id oldSuperToMaskedIndex =
+        vtkm::worklet::contourtree_augmented::MaskedIndex(oldSuperToValue);
+
+      // since we haven't set up the regular search array yet, we can't use that
+      // instead, we know that the two supernodes must be in the new tree, so we retrieve their new super IDs
+      // and convert them to regular
+
+      // retrieve their new super IDs
+      vtkm::Id newSuperFrom = createSuperarcsData.NewSupernodeIds.Get(oldSuperFromValue);
+      vtkm::Id newSuperTo = createSuperarcsData.NewSupernodeIds.Get(oldSuperToMaskedIndex);
+
+      // convert to regular IDs (which is what the FindSuperArcForUnknownNode() routine assumes)
+      vtkm::Id newRegularFrom = augmentedTreeSupernodesPortal.Get(newSuperFrom);
+      vtkm::Id newRegularTo = augmentedTreeSupernodesPortal.Get(newSuperTo);
+
+      // the new superparent after the search
+      vtkm::Id newSuperparentId = vtkm::worklet::contourtree_augmented::NO_SUCH_ELEMENT;
+
+      // depending on the ascending flag
+      if (ascendingSuperarc)
       {
-        augmentedTreeSuperarcsValue = vtkm::worklet::contourtree_augmented::NO_SUCH_ELEMENT;
+        newSuperparentId = findSuperArcForUnknownNode.FindSuperArcForUnknownNode(
+          globalRegularIdSetValue, dataValueSetValue, newRegularTo, newRegularFrom);
       }
       else
-      { // not the tree root
-        // retrieve the target of the superarc from the base tree (masking to strip out the ascending flag)
-        vtkm::Id oldTargetSuperId = vtkm::worklet::contourtree_augmented::MaskedIndex(
-          baseTreeSuperarcsPortal.Get(superparentOldSuperId));
-        // convert to a new supernode Id
-        vtkm::Id newTargetSuperId = newSupernodeIdsPortal.Get(oldTargetSuperId);
-        // add the ascending flag back in and store in the array
-        augmentedTreeSuperarcsValue = newTargetSuperId |
-          (superarcAscends ? vtkm::worklet::contourtree_augmented::IS_ASCENDING : 0x00);
-      } // not the tree root
-      // since there's an extra entry in the firstSupernode array as a sentinel, set it
-      augmentedTreeFirstSupernodePerIterationPortal.Set(this->AugmentedTreeNumIterations,
-                                                        NumAugmentedTreeSupernodes);
-    } // last in the array
-    else if (superparentOldSuperId !=
-             vtkm::worklet::contourtree_augmented::MaskedIndex(
-               superparentSetPortal.Get(supernodeSorterPortal.Get(supernode + 1))))
-    { // last in the segment
-      // retrieve the target of the superarc from the base tree (masking to strip out the ascending flag)
-      vtkm::Id oldTargetSuperId = vtkm::worklet::contourtree_augmented::MaskedIndex(
-        baseTreeSuperarcsPortal.Get(superparentOldSuperId));
-      // convert to a new supernode Id
-      vtkm::Id newTargetSuperId = newSupernodeIdsPortal.Get(oldTargetSuperId);
-      // add the ascending flag back in and store in the array
-      augmentedTreeSuperarcsValue = newTargetSuperId |
-        (superarcAscends ? vtkm::worklet::contourtree_augmented::IS_ASCENDING : 0x00);
+      {
+        newSuperparentId = findSuperArcForUnknownNode.FindSuperArcForUnknownNode(
+          globalRegularIdSetValue, dataValueSetValue, newRegularFrom, newRegularTo);
+      }
 
-      // since we're the last in the segment, we check to see if we are at the end of an iteration
-      vtkm::Id iterationNumber = vtkm::worklet::contourtree_augmented::MaskedIndex(
-        baseTreeWhichIterationPortal.Get(superparentOldSuperId));
-      vtkm::Id iterationNumberOfNext = vtkm::worklet::contourtree_augmented::MaskedIndex(
-        baseTreeWhichIterationPortal.Get(superparentOldSuperId + 1));
+      // attachment points use the superparent to store the superarc they insert onto
+      augmentedTreeSuperparentsValue = newSuperparentId;
 
-      if (iterationNumber != iterationNumberOfNext)
-      { // boundary of iterations
-        // If so, we set the "firstSupernodePerIteration" for the next
-        augmentedTreeFirstSupernodePerIterationPortal.Set(iterationNumberOfNext,
-                                                          newSupernodeId + 1);
-      } // boundary of iterations
-    }   // last in the segment
+    } // preserved attachment point
     else
-    { // not last in the segment
-      // the target is always the next one, so just store it with the ascending flag
-      augmentedTreeSuperarcsValue = (newSupernodeId + 1) |
-        (superarcAscends ? vtkm::worklet::contourtree_augmented::IS_ASCENDING : 0x00);
-    } // not last in the segment
+    { // raised attachment point or "ordinary" supernodes
+      // Since all of the superparents must be in the base tree, we can now retrieve the target
+      vtkm::Id superparentOldSuperId = vtkm::worklet::contourtree_augmented::MaskedIndex(
+        createSuperarcsData.SuperparentSet.Get(supernodeSetIndex));
+      vtkm::Id oldTargetSuperId = createSuperarcsData.BaseTreeSuperarcs.Get(superparentOldSuperId);
+      // and break it into a target and flags
+      bool ascendingSuperarc = vtkm::worklet::contourtree_augmented::IsAscending(oldTargetSuperId);
+      // NOTE: if the target was NO_SUCH_ELEMENT, this will hold 0
+      oldTargetSuperId = vtkm::worklet::contourtree_augmented::MaskedIndex(oldTargetSuperId);
 
-    // set the first supernode in the first iteration to the beginning of the round
-    augmentedTreeFirstSupernodePerIterationPortal.Set(0, this->NumSupernodesAlready);
+      // and another boolean for whether we are the last element in a segment
+      bool isLastInSegment = false;
 
-
-    // NOTE: This part has been moved out of the worklet and is performed using standard vtkm copy constructs
-    // // setting the hyperparent is straightforward since the hyperstructure is preserved
-    // // we take the superparent (which is guaranteed to be in the baseTree), find it's hyperparent and use that
-    // augmentedTreeHyperparentsPortal.Set(newSupernodeId, baseTreeHyperparentsPortal.Get(superparentOldSuperId));
-
-    // NOTE: This part could potentially be made a separate worklet but it does not seem necessary
-    // similarly, the super2hypernode should carry over, but it's harder to test because of the attachment points which
-    // do not have valid old supernode Ids.  Instead, we check their superparent's regular global Id against them: if it
-    // matches, then it must be the start of the superarc, in which case it does have an old Id, and we can then use the
-    // existing hypernode Id
-    vtkm::Id superparentOldRegularId = baseTreeSupernodesPortal.Get(superparentOldSuperId);
-    vtkm::Id superparentGlobalId = baseTreeRegularNodeGlobalIdsPortal.Get(superparentOldRegularId);
-    // Here: globalRegularIdSetValue is the same as globalRegularIdSetPortal.Get(supernodeSetIndex)
-    if (superparentGlobalId == globalRegularIdSetValue)
-    {
-      // augmentedTreeSuper2hypernodePortal.Set(newSupernodeId, baseTreeSuper2HypernodePortal.Get(superparentOldSuperId));
-      augmentedTreeSuper2hypernodeValue = baseTreeSuper2HypernodePortal.Get(superparentOldSuperId);
-    }
-    else
-    {
-      // augmentedTreeSuper2hypernodePortal.Set(newSupernodeId, vtkm::worklet::contourtree_augmented::NO_SUCH_ELEMENT);
-      augmentedTreeSuper2hypernodeValue = vtkm::worklet::contourtree_augmented::NO_SUCH_ELEMENT;
-    }
-
-    // NOTE: This part has been moved out of the worklet and is performed using standard vtkm copy constructs
-    // // which round and iteration carry over
-    // augmentedTreeWhichRoundPortal.Set(newSupernodeId, baseTreeWhichRoundPortal.Get(superparentOldSuperId));
-    // augmentedTreeWhichIterationPortal.Set(newSupernodeId, baseTreeWhichIterationPortal.Get(superparentOldSuperId));
-
-    // now we deal with the regular-sized arrays
-
-    // NOTE: This part has been moved out of the worklet and is performed using standard vtkm copy constructs
-    // // copy the global regular Id and data value
-    // augmentedTreeRegularNodeGlobalIdsPortal.Set(newRegularId, globalRegularIdSetPortal.Get(supernodeSetIndex));
-    // augmentedTreeDataValuesPortal.Set(newRegularId, dataValueSetPortal.Get(supernodeSetIndex));
-
-    // NOTE: This part has been moved out of the worklet and is performed using standard vtkm copy constructs
-    // // the sort order will be dealt with later
-    // // since all of these nodes are supernodes, they will be their own superparent, which means that:
-    // //  a.  the regular2node can be set immediately
-    // augmentedTreeRegular2SupernodePortal.Set(newRegularId, newSupernodeId);
-    // //  b.  as can the superparent
-    // augmentedTreeSuperparentsPortal.Set(newRegularId, newSupernodeId);
-
-    // In serial this worklet implements the following operation
-    /*
-    for (vtkm::Id supernode = 0; supernode < supernodeSorter.size(); supernode++)
-    { // per supernode in the set
-      // retrieve the index from the sorting index array
-      vtkm::Id supernodeSetIndex = supernodeSorter[supernode];
-
-      // work out the new supernode ID
-      vtkm::Id newSupernodeID = numSupernodesAlready + supernode;
-
-      //  At all levels above 0, we used to keep regular vertices in case they are attachment points.  After augmentation, we don't need to.
-      //  Instead, at all levels above 0, the regular nodes in each round are identical to the supernodes
-      //  In order to avoid confusion, we will copy the ID into a separate variable
-      vtkm::Id newRegularID = newSupernodeID;
-
-      // setting the supernode's regular ID is now trivial
-      augmentedTree->supernodes      [newSupernodeID] = newRegularID;
-
-      // retrieve the ascending flag from the superparent
-      bool superarcAscends = isAscending(superparentSet[supernodeSetIndex]);
-
-      // strip the ascending flag from the superparent
-      vtkm::Id superparentOldSuperID = maskedIndex(superparentSet[supernodeSetIndex]);
+      // end of the entire array counts as last in segment
+      if (supernode == supernodeSorterPortal.GetNumberOfValues() - 1)
+      {
+        isLastInSegment = true;
+      }
+      // otherwise, check for a mismatch in the sorting superparent which indicates the end of a segment
+      else if (vtkm::worklet::contourtree_augmented::MaskedIndex(
+                 createSuperarcsData.SuperparentSet.Get(supernodeSetIndex)) !=
+               vtkm::worklet::contourtree_augmented::MaskedIndex(
+                 createSuperarcsData.SuperparentSet.Get(supernodeSorterPortal.Get(supernode + 1))))
+      {
+        isLastInSegment = true;
+      }
 
       // setting the superarc is done the usual way.  Our sort routine has ended up with the supernodes arranged in either ascending or descending order
       // inwards along the parent superarc (as expressed by the superparent ID).  Each superarc except the last in the segment points to the next one:
       // the last one points to the target of the original superarc.
-      // first test to see if we're the last in the array
-      if (supernode == supernodeSorter.size() - 1)
-      { // last in the array
-        // special case for root of entire tree at end of top level
-        if (roundNumber == baseTree->nRounds)
-        {
-          augmentedTree->superarcs[newSupernodeID] = NO_SUCH_ELEMENT;
-        }
-        else
-        { // not the tree root
-          // retrieve the target of the superarc from the base tree (masking to strip out the ascending flag)
-          vtkm::Id oldTargetSuperID = maskedIndex(baseTree->superarcs[superparentOldSuperID]);
-          // convert to a new supernode ID
-          vtkm::Id newTargetSuperID = newSupernodeIDs[oldTargetSuperID];
-          // add the ascending flag back in and store in the array
-          augmentedTree->superarcs[newSupernodeID] = newTargetSuperID | (superarcAscends ? IS_ASCENDING : 0x00);
-        } // not the tree root
-        // since there's an extra entry in the firstSupernode array as a sentinel, set it
-        augmentedTree->firstSupernodePerIteration[roundNumber][augmentedTree->nIterations[roundNumber]] = augmentedTree->supernodes.size();
-      } // last in the array
-      else if (superparentOldSuperID != maskedIndex(superparentSet[supernodeSorter[supernode+1]]))
-      { // last in the segment
-        // retrieve the target of the superarc from the base tree (masking to strip out the ascending flag)
-        vtkm::Id oldTargetSuperID = maskedIndex(baseTree->superarcs[superparentOldSuperID]);
-        // convert to a new supernode ID
-        vtkm::Id newTargetSuperID = newSupernodeIDs[oldTargetSuperID];
-        // add the ascending flag back in and store in the array
-        augmentedTree->superarcs[newSupernodeID] = newTargetSuperID | (superarcAscends ? IS_ASCENDING : 0x00);
-
-        // since we're the last in the segment, we check to see if we are at the end of an iteration
-        vtkm::Id iterationNumber     = maskedIndex(baseTree->whichIteration[superparentOldSuperID]);
-        vtkm::Id iterationNumberOfNext = maskedIndex(baseTree->whichIteration[superparentOldSuperID + 1]);
-
-        if (iterationNumber != iterationNumberOfNext)
-        { // boundary of iterations
-          // If so, we set the "firstSupernodePerIteration" for the next
-          augmentedTree->firstSupernodePerIteration[roundNumber][iterationNumberOfNext] = newSupernodeID + 1;
-        } // boundary of iterations
-      } // last in the segment
+      if (isLastInSegment)
+      { // last in segment
+        // we take the old target of the superarc (in old supernode IDs) and convert it to a new supernode ID
+        augmentedTreeSuperarcsValue = createSuperarcsData.NewSupernodeIds.Get(oldTargetSuperId) |
+          (ascendingSuperarc ? vtkm::worklet::contourtree_augmented::IS_ASCENDING : 0x00);
+      } // last in segment
       else
-      { // not last in the segment
+      { // not last in segment
         // the target is always the next one, so just store it with the ascending flag
-        augmentedTree->superarcs[newSupernodeID] = (newSupernodeID+1) | (superarcAscends ? IS_ASCENDING : 0x00);
-      } // not last in the segment
+        augmentedTreeSuperarcsValue = (newSupernodeId + 1) |
+          (ascendingSuperarc ? vtkm::worklet::contourtree_augmented::IS_ASCENDING : 0x00);
+      } // not last in segment
 
-      // set the first supernode in the first iteration to the beginning of the round
-      augmentedTree->firstSupernodePerIteration[roundNumber][0] = numSupernodesAlready;
+      // first we identify the hyperarc on which the superarc sits
+      // this will be visible in the old base tree, since hyperstructure carries over
+      vtkm::Id oldHyperparent = createSuperarcsData.BaseTreeHyperparents.Get(superparentOldSuperId);
 
-      // setting the hyperparent is straightforward since the hyperstructure is preserved
-      // we take the superparent (which is guaranteed to be in the baseTree), find it's hyperparent and use that
-      augmentedTree->hyperparents      [newSupernodeID] = baseTree->hyperparents      [superparentOldSuperID];
+      // hyperstructure carries over, so we use the same hyperparent as the superparent
+      augmentedTreeHyperparentsValue = oldHyperparent;
 
-      // similarly, the super2hypernode should carry over, but it's harder to test because of the attachment points which
-      // do not have valid old supernode IDs.  Instead, we check their superparent's regular global ID against them: if it
-      // matches, then it must be the start of the superarc, in which case it does have an old ID, and we can then use the
-      // existing hypernode ID
-      vtkm::Id superparentOldRegularID = baseTree->supernodes[superparentOldSuperID];
-      vtkm::Id superparentGlobalID = baseTree->regularNodeGlobalIDs[superparentOldRegularID];
-      if (superparentGlobalID == globalRegularIDSet[supernodeSetIndex])
+      // retrieve the hyperparent's old supernode ID & convert to a new one, then test it
+      if (createSuperarcsData.NewSupernodeIds.Get(
+            createSuperarcsData.BaseTreeHypernodes.Get(oldHyperparent)) == newSupernodeId)
       {
-        augmentedTree->super2hypernode  [newSupernodeID]   = baseTree->super2hypernode[superparentOldSuperID];
+        augmentedTreeSuper2HypernodeValue = oldHyperparent;
       }
       else
       {
-        augmentedTree->super2hypernode  [newSupernodeID]   = NO_SUCH_ELEMENT;
+        augmentedTreeSuper2HypernodeValue = vtkm::worklet::contourtree_augmented::NO_SUCH_ELEMENT;
       }
 
-      // which round and iteration carry over
-      augmentedTree->whichRound      [newSupernodeID]   = baseTree->whichRound[superparentOldSuperID];
-      augmentedTree->whichIteration    [newSupernodeID]   = baseTree->whichIteration[superparentOldSuperID];
+      // round and iteration are set from the superparent, since we are raising to its level
+      augmentedTreeWhichRoundValue =
+        createSuperarcsData.BaseTreeWhichRound.Get(superparentOldSuperId);
+      augmentedTreeWhichIterationValue =
+        createSuperarcsData.BaseTreeWhichIteration.Get(superparentOldSuperId);
+      // and set the relevant regular arrays
+      augmentedTreeRegularNodeGlobalIdsValue = globalRegularIdSetValue;
+      augmentedTreeDataValuesValue = dataValueSetValue;
+      // for all supernodes, this points to itself
+      augmentedTreeRegular2SupernodeValue = newSupernodeId;
+      // and since we're inserted, so does this
+      augmentedTreeSuperparentsValue = newSupernodeId;
+    } // raised attachment point or "ordinary" supernodes
 
-      // now we deal with the regular-sized arrays
+    /*  Original PPP2 code
+    for (indexType supernode = 0; supernode < supernodeSorter.size(); supernode++)
+        { // per supernode in the set
+        // retrieve the index from the sorting index array
+        indexType supernodeSetIndex = supernodeSorter[supernode];
 
-      // copy the global regular ID and data value
-      augmentedTree->regularNodeGlobalIDs  [newRegularID]    = globalRegularIDSet[supernodeSetIndex];
-      augmentedTree->dataValues      [newRegularID]    = dataValueSet[supernodeSetIndex];
+        // work out the new supernode ID
+        indexType newSupernodeID = nSupernodesAlready + supernode;
 
-      // the sort order will be dealt with later
-      // since all of these nodes are supernodes, they will be their own superparent, which means that:
-      //  a.  the regular2node can be set immediately
-      augmentedTree->regular2supernode  [newRegularID]     = newSupernodeID;
-      //  b.  as can the superparent
-      augmentedTree->superparents      [newRegularID]     = newSupernodeID;
-    } // per supernode in the set
+        // and the old supernode ID
+        // NB: May be NO_SUCH_ELEMENT if not already in the base tree
+        indexType oldSupernodeID = supernodeIDSet[supernodeSetIndex];
+
+        //	At all levels above 0, we used to keep regular vertices in case they are attachment points.  After augmentation, we don't need to.
+        //	Instead, at all levels above 0, the regular nodes in each round are identical to the supernodes
+        //	In order to avoid confusion, we will copy the ID into a separate variable
+        indexType newRegularID = newSupernodeID;
+
+        // setting the supernode's regular ID is now trivial
+        augmentedTree->supernodes			[newSupernodeID] = newRegularID;
+
+        // retrieve the old superID of the superparent.  This is slightly tricky, as we have four classes of supernodes:
+        // 1. the root of the entire tree
+        // 2. attachment points not being inserted. In this case, the supernode ID is stored in the superparentSet array, not the superparent for insertion purposes
+        // 3. attachment points being inserted.  In this case, the superparent is stored in the superparentSet array
+        // 4. "ordinary" supernodes, where the superparent is the same as the supernode ID anyway
+        //
+        // Note that an attachment point gets inserted into a parent superarc.  But the attachment point itself has a NULL superarc, because it's only a virtual insertion
+        // This means that such an attachment superarc cannot be the superparent of any other attachment point
+        // It is therefore reasonable to deal with 1. & 2 separately. 3. & 4. then combine together
+
+        // first we test for the root of the tree
+        if ((roundNo == baseTree->nRounds) && (supernode == supernodeSorter.size() - 1))
+            { // root of the tree
+            // note that oldSupernodeID is guaranteed not to be NO_SUCH_ELEMENT, as the root is in every tree
+            // set the super arrays
+            augmentedTree->superarcs			[newSupernodeID]	= NO_SUCH_ELEMENT;
+            // hyperstructure carries over, so we use the same hyperparent as before
+            augmentedTree->hyperparents			[newSupernodeID]	= baseTree->hyperparents[oldSupernodeID];
+            // and set the hypernode ID
+            augmentedTree->super2hypernode		[newSupernodeID]	= baseTree->super2hypernode[oldSupernodeID];
+            // and the round and iteration
+            augmentedTree->whichRound			[newSupernodeID]	= baseTree->whichRound[oldSupernodeID];
+            augmentedTree->whichIteration		[newSupernodeID]	= baseTree->whichIteration[oldSupernodeID];
+            // and set the relevant regular arrays
+            augmentedTree->regularNodeGlobalIDs	[newRegularID]		= globalRegularIDSet[supernodeSetIndex];
+            augmentedTree->dataValues			[newRegularID]		= dataValueSet[supernodeSetIndex];
+            // for the root, these always point to itself
+            augmentedTree->regular2supernode	[newRegularID]		= newSupernodeID;
+            augmentedTree->superparents			[newRegularID]		= newSupernodeID;
+            } // root of the tree
+        // now deal with unsimplified attachment points, which we can identify because they were in the "kept" batch, not the "inserted" batch,
+        // and this is given away by the index into the set of supernodes to be added
+        // and the fact that the superarc is NO_SUCH_ELEMENT
+        else if	((supernodeSetIndex >= nInsertedSupernodes) && (noSuchElement(baseTree->superarcs[oldSupernodeID])))
+            { // preserved attachment point
+            // note that oldSupernodeID is guaranteed not to be NO_SUCH_ELEMENT, as the supernode came from this block originally
+            // set the superarc to NO_SUCH_ELEMENT, as before
+            augmentedTree->superarcs			[newSupernodeID]	= NO_SUCH_ELEMENT;
+            // hyperstructure carries over, so we use the same hyperparent as before
+            augmentedTree->hyperparents			[newSupernodeID]	= baseTree->hyperparents[oldSupernodeID];
+            // attachment points are never hypernodes anyway, so set it directly
+            augmentedTree->super2hypernode		[newSupernodeID]	= NO_SUCH_ELEMENT;
+            // and the round and iteration
+            augmentedTree->whichRound			[newSupernodeID]	= baseTree->whichRound[oldSupernodeID];
+            augmentedTree->whichIteration		[newSupernodeID]	= baseTree->whichIteration[oldSupernodeID];
+            // and set the relevant regular arrays
+            augmentedTree->regularNodeGlobalIDs	[newRegularID]		= globalRegularIDSet[supernodeSetIndex];
+            augmentedTree->dataValues			[newRegularID]		= dataValueSet[supernodeSetIndex];
+            // for a preserved attachment point, this always points to itself
+            augmentedTree->regular2supernode	[newRegularID]		= newSupernodeID;
+            // the superparent is the tricky one, as the old one may have been broken up by insertions at a higher level
+
+            // Here, what we need to do is a search in the augmented tree to find which superarc to attach to.  This is necessary
+            // because the old superarc it attached to may have been broken up.
+            // We are guaranteed that there is one, and that it only uses the higher levels of the augmented tree,
+            // so the fact that we are partially constructed doesn't get in the way.  To do this, we need supernodes
+            // known to be in the higher level that are above and below the supernode.
+            // Since the point was an attachment point in the base tree, that means that there is a higher round superarc
+            // it inserts into.  Moreover, the algorithm ALWAYS inserts a supernode at or above its original round, so
+            // we can guarantee that both ends of the parent are in the higher levels.  Which means we only need to work
+            // out which end is higher.
+            indexType oldRegularID = baseTree->supernodes[oldSupernodeID];
+            indexType oldSuperFrom = baseTree->superparents[oldRegularID];
+            indexType oldSuperTo = baseTree->superarcs[oldSuperFrom];
+            // retrieve the ascending flag
+            bool ascendingSuperarc = isAscending(oldSuperTo);
+            // and mask out the flags
+            oldSuperTo = maskedIndex(oldSuperTo);
+
+            // since we haven't set up the regular search array yet, we can't use that
+            // instead, we know that the two supernodes must be in the new tree, so we retrieve their new super IDs
+            // and convert them to regular
+
+            // retrieve their new super IDs
+            indexType newSuperFrom = newSupernodeIDs[oldSuperFrom];
+            indexType newSuperTo = newSupernodeIDs[oldSuperTo];
+
+            // convert to regular IDs (which is what the FindSuperArcForUnknownNode() routine assumes)
+            indexType newRegularFrom = augmentedTree->supernodes[newSuperFrom];
+            indexType newRegularTo = augmentedTree->supernodes[newSuperTo];
+
+            // the new superparent after the search
+            indexType newSuperparentID = NO_SUCH_ELEMENT;
+
+            // depending on the ascending flag
+            if (ascendingSuperarc)
+                newSuperparentID = augmentedTree->FindSuperArcForUnknownNode(globalRegularIDSet[supernodeSetIndex],dataValueSet[supernodeSetIndex], newRegularTo, newRegularFrom);
+            else
+                newSuperparentID = augmentedTree->FindSuperArcForUnknownNode(globalRegularIDSet[supernodeSetIndex],dataValueSet[supernodeSetIndex], newRegularFrom, newRegularTo);
+
+            // attachment points use the superparent to store the superarc they insert onto
+            augmentedTree->superparents			[newRegularID] 		= newSuperparentID;
+
+            } // preserved attachment point
+        else
+            { // raised attachment point or "ordinary" supernodes
+            // Since all of the superparents must be in the base tree, we can now retrieve the target
+            indexType superparentOldSuperID = maskedIndex(superparentSet[supernodeSetIndex]);
+            indexType oldTargetSuperID = baseTree->superarcs[superparentOldSuperID];
+            // and break it into a target and flags
+            bool ascendingSuperarc = isAscending(oldTargetSuperID);
+            // NOTE: if the target was NO_SUCH_ELEMENT, this will hold 0
+            oldTargetSuperID = maskedIndex(oldTargetSuperID);
+
+            // and another boolean for whether we are the last element in a segment
+            bool isLastInSegment = false;
+
+            // end of the entire array counts as last in segment
+            if (supernode == supernodeSorter.size() - 1)
+                isLastInSegment = true;
+            // otherwise, check for a mismatch in the sorting superparent which indicates the end of a segment
+            else if (maskedIndex(superparentSet[supernodeSetIndex]) != maskedIndex(superparentSet[supernodeSorter[supernode+1]]))
+                isLastInSegment = true;
+
+            // setting the superarc is done the usual way.  Our sort routine has ended up with the supernodes arranged in either ascending or descending order
+            // inwards along the parent superarc (as expressed by the superparent ID).  Each superarc except the last in the segment points to the next one:
+            // the last one points to the target of the original superarc.
+            if (isLastInSegment)
+                { // last in segment
+                // we take the old target of the superarc (in old supernode IDs) and convert it to a new supernode ID
+                augmentedTree->superarcs[newSupernodeID] = newSupernodeIDs[oldTargetSuperID] | (ascendingSuperarc ? IS_ASCENDING : 0x00);
+                } // last in segment
+            else
+                { // not last in segment
+                // the target is always the next one, so just store it with the ascending flag
+                augmentedTree->superarcs[newSupernodeID] = (newSupernodeID+1) | (ascendingSuperarc ? IS_ASCENDING : 0x00);
+                } // not last in segment
+
+            // first we identify the hyperarc on which the superarc sits
+            // this will be visible in the old base tree, since hyperstructure carries over
+            indexType oldHyperparent = baseTree->hyperparents[superparentOldSuperID];
+
+            // hyperstructure carries over, so we use the same hyperparent as the superparent
+            augmentedTree->hyperparents			[newSupernodeID]	= oldHyperparent;
+            // retrieve the hyperparent's old supernode ID & convert to a new one, then test it
+            if (newSupernodeIDs[baseTree->hypernodes[oldHyperparent]] == newSupernodeID)
+                augmentedTree->super2hypernode	[newSupernodeID]	= oldHyperparent;
+            else
+                augmentedTree->super2hypernode	[newSupernodeID]	= NO_SUCH_ELEMENT;
+
+            // round and iteration are set from the superparent, since we are raising to its level
+            augmentedTree->whichRound			[newSupernodeID]	= baseTree->whichRound[superparentOldSuperID];
+            augmentedTree->whichIteration		[newSupernodeID]	= baseTree->whichIteration[superparentOldSuperID];
+            // and set the relevant regular arrays
+            augmentedTree->regularNodeGlobalIDs	[newRegularID]		= globalRegularIDSet[supernodeSetIndex];
+            augmentedTree->dataValues			[newRegularID]		= dataValueSet[supernodeSetIndex];
+            // for all supernodes, this points to itself
+            augmentedTree->regular2supernode	[newRegularID]		= newSupernodeID;
+            // and since we're inserted, so does this
+            augmentedTree->superparents			[newRegularID]		= newSupernodeID;
+            } // raised attachment point or "ordinary" supernodes
+
+        } // per supernode in the set
     */
+
+
   } // operator()()
 
 private:
   const vtkm::Id NumSupernodesAlready;
   const vtkm::Id BaseTreeNumRounds;
-  const vtkm::Id AugmentedTreeNumIterations;
-  const vtkm::Id RoundNumber;
-  const vtkm::Id NumAugmentedTreeSupernodes;
+  const vtkm::Id NumInsertedSupernodes;
+  const vtkm::Id RoundNo;
 
 }; // CreateSuperarcsWorklet
 

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

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

@@ -81,27 +81,50 @@ public:
     const vtkm::worklet::contourtree_augmented::IdArrayType& superarcs,
     const vtkm::worklet::contourtree_augmented::IdArrayType& whichRound,
     const vtkm::Id numRounds,
+    vtkm::worklet::contourtree_augmented::IdArrayType* volumeArray,
+    vtkm::Id presimplifyThreshold,
     vtkm::cont::DeviceAdapterId device,
     vtkm::cont::Token& token)
     : SuperarcsPortal(superarcs.PrepareForInput(device, token))
     , WhichRoundPortal(whichRound.PrepareForInput(device, token))
     , NumRounds(numRounds)
+    , PresimplifyThreshold(presimplifyThreshold)
   { // constructor
+    this->Presimplify = ((volumeArray != NULL) && (presimplifyThreshold > 0));
+    // If we presimplify then store the volumeArray. Otherwise we don't need to volume array and we
+    // set it to another portal, just to make sure the variable is being initalized with something
+    this->VolumeArrayPortal =
+      this->Presimplify ? volumeArray->PrepareForInput(device, token) : this->WhichRoundPortal;
   } // constructor
 
   // () operator - gets called to do comparison
   VTKM_EXEC
   bool operator()(const vtkm::Id& supernode) const
   { // operator()
-    return (
-      vtkm::worklet::contourtree_augmented::NoSuchElement(this->SuperarcsPortal.Get(supernode)) &&
-      (this->WhichRoundPortal.Get(supernode) < this->NumRounds));
+    // an attachment point is defined by having no superarc (NO_SUCH_ELEMENT) and not being in
+    // the final round (where this indicates the global root)
+    bool predicate =
+      (vtkm::worklet::contourtree_augmented::NoSuchElement(this->SuperarcsPortal.Get(supernode)) &&
+       (this->WhichRoundPortal.Get(supernode) < this->NumRounds));
+    // if we pass this check then we need to also check that the supernode passes the pre-simplification threshold
+    if (predicate && this->Presimplify)
+    {
+      // suppress if it's volume is at or below the threshold
+      if (this->VolumeArrayPortal.Get(supernode) <= this->PresimplifyThreshold)
+      {                    // below threshold
+        predicate = false; // do not keep attachement point below the simplification threshold
+      }                    // below threshold
+    }
+    return predicate;
   } // operator()
 
 private:
   IdPortalType SuperarcsPortal;
   IdPortalType WhichRoundPortal;
   const vtkm::Id NumRounds;
+  bool Presimplify;
+  IdPortalType VolumeArrayPortal;
+  vtkm::Id PresimplifyThreshold;
 
 
 }; // IsAttachementPointPredicateImpl
@@ -113,24 +136,35 @@ public:
   VTKM_CONT
   IsAttachementPointPredicate(const vtkm::worklet::contourtree_augmented::IdArrayType& superarcs,
                               const vtkm::worklet::contourtree_augmented::IdArrayType& whichRound,
-                              const vtkm::Id numRounds)
+                              const vtkm::Id numRounds,
+                              vtkm::worklet::contourtree_augmented::IdArrayType* volumeArray = NULL,
+                              vtkm::Id presimplifyThreshold = 0)
     : Superarcs(superarcs)
     , WhichRound(whichRound)
     , NumRounds(numRounds)
+    , VolumeArray(volumeArray)
+    , PresimplifyThreshold(presimplifyThreshold)
   {
   }
 
   VTKM_CONT IsAttachementPointPredicateImpl PrepareForExecution(vtkm::cont::DeviceAdapterId device,
                                                                 vtkm::cont::Token& token) const
   {
-    return IsAttachementPointPredicateImpl(
-      this->Superarcs, this->WhichRound, this->NumRounds, device, token);
+    return IsAttachementPointPredicateImpl(this->Superarcs,
+                                           this->WhichRound,
+                                           this->NumRounds,
+                                           this->VolumeArray,
+                                           this->PresimplifyThreshold,
+                                           device,
+                                           token);
   }
 
 private:
   vtkm::worklet::contourtree_augmented::IdArrayType Superarcs;
   vtkm::worklet::contourtree_augmented::IdArrayType WhichRound;
   const vtkm::Id NumRounds;
+  vtkm::worklet::contourtree_augmented::IdArrayType* VolumeArray;
+  vtkm::Id PresimplifyThreshold;
 }; // IsAttachementPointPredicate
 
 } // namespace hierarchical_augmenter

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

@@ -65,12 +65,16 @@ public:
   /// Control signature for the worklet
   using ControlSignature = void(
     FieldIn supernodeIndex, // input domain ArrayHandleIndex(SupernodeSorter.GetNumberOfValues())
+    // FieldIn supernodeIdSetPermuted, // input supernodeIDSet permuted by supernodeSorter to allow for FieldIn
     FieldIn
-      supernodeIdSetPermuted, // input supernodeIDSet permuted by supernodeSorter to allow for FieldIn
+      globalRegularIdSet, // input globalRegularIdSet permuted by supernodeSorter to allow for FieldIn
+    ExecObject findRegularByGlobal,
+    WholeArrayIn baseTreeRegular2Supernode,
     WholeArrayInOut
       newSupernodeIds // output/input (both are necessary since not all valyes will be overwritten)
   );
-  using ExecutionSignature = void(_1, _2, _3);
+
+  using ExecutionSignature = void(_1, _2, _3, _4, _5);
   using InputDomain = _1;
 
   // Default Constructor
@@ -80,10 +84,13 @@ public:
   {
   }
 
-  template <typename InOutFieldPortalType>
+  template <typename InOutFieldPortalType, typename InFieldPortalType, typename ExecObjectType>
   VTKM_EXEC void operator()(
-    const vtkm::Id& supernode,      // InputIndex of supernodeSorter
-    const vtkm::Id& oldSupernodeId, // same as supernodeIDSet[supernodeSetIndex];
+    const vtkm::Id& supernode, // InputIndex of supernodeSorter
+    // const vtkm::Id& oldSupernodeId, // same as supernodeIDSet[supernodeSetIndex];
+    const vtkm::Id& globalRegularIdSetValue, // same as globalRegularIDSet[supernodeSetIndex]
+    const ExecObjectType& findRegularByGlobal,
+    const InFieldPortalType& baseTreeRegular2SupernodePortal,
     const InOutFieldPortalType& newSupernodeIdsPortal) const
   {
     // per supernode
@@ -96,6 +103,17 @@ public:
     // that if it came from another block it will be set to NO_SUCH_ELEMENT
     // vtkm::Id oldSupernodeId set on input since we use ArrayHandlePermutation to
     // shuffle supernodeIDSet by supernodeSorter;
+    // TODO/WARNING:  Logic error in that comment for presimplified trees, but not for the original version.  See RetrieveOldSupernodes() for why.
+
+    // TODO/WARNING:  We substitute a search in the old hierarchical tree for the supernode.  If it is present, then we fill in it's entry in the
+
+    // newSupernodeIDs array. If not, we're happy.
+    vtkm::Id oldRegularId = findRegularByGlobal.FindRegularByGlobal(globalRegularIdSetValue);
+    vtkm::Id oldSupernodeId = vtkm::worklet::contourtree_augmented::NO_SUCH_ELEMENT;
+    if (!vtkm::worklet::contourtree_augmented::NoSuchElement(oldRegularId))
+    {
+      oldSupernodeId = baseTreeRegular2SupernodePortal.Get(oldRegularId);
+    }
 
     // and write to the lookup array
     if (!vtkm::worklet::contourtree_augmented::NoSuchElement(oldSupernodeId))

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

@@ -66,10 +66,12 @@ class UpdateHyperstructureSetSuperchildrenWorklet : public vtkm::worklet::Workle
 public:
   /// Control signature for the worklet
   using ControlSignature = void(
-    WholeArrayIn augmentedTreeHypernodes, // input (we need both this and the next value)
-    FieldOut augmentedTreeSuperchildren   // output
+    WholeArrayIn augmentedTreeHypernodes,      // input (we need both this and the next value)
+    FieldIn augmentedTreeSuperarcs,            // input
+    WholeArrayIn augmentedTreeHyperparents,    // input
+    WholeArrayInOut augmentedTreeSuperchildren // output
   );
-  using ExecutionSignature = void(InputIndex, _1, _2);
+  using ExecutionSignature = void(InputIndex, _1, _2, _3, _4);
   using InputDomain = _1;
 
   // Default Constructor
@@ -80,46 +82,58 @@ public:
   }
 
 
-  template <typename InFieldPortalType>
-  VTKM_EXEC void operator()(
-    const vtkm::Id& hypernode,
-    const InFieldPortalType& augmentedTreeHypernodesPortal,
-    vtkm::Id&
-      augmentedTreeSuperchildrenValue // same as augmentedTree->superchildren[InputIndex] = ...
-  ) const
+  template <typename InFieldPortalType1, typename InFieldPortalType2, typename OutFieldPortalType>
+  VTKM_EXEC void operator()(const vtkm::Id& supernode,
+                            const InFieldPortalType1& augmentedTreeHypernodesPortal,
+                            const vtkm::Id& augmentedTreeSuperarcsValue,
+                            const InFieldPortalType2& augmentedTreeHyperparentsPortal,
+                            const OutFieldPortalType& augmentedTreeSuperchildrenPortal) const
   {
-    // per hypernode
-    // retrieve the new superId
-    vtkm::Id superId = augmentedTreeHypernodesPortal.Get(hypernode);
-    // and the next one over
-    vtkm::Id nextSuperId;
-    if (hypernode == augmentedTreeHypernodesPortal.GetNumberOfValues() - 1)
+    // per supernode
+    // attachment points have NULL superarcs and are skipped
+    if (vtkm::worklet::contourtree_augmented::NoSuchElement(augmentedTreeSuperarcsValue))
     {
-      nextSuperId = this->AugmentedTreeNumSupernodes;
+      return;
     }
-    else
-    {
-      nextSuperId = augmentedTreeHypernodesPortal.Get(hypernode + 1);
-    }
-    // the difference is the number of superchildren
-    augmentedTreeSuperchildrenValue = nextSuperId - superId;
+    // we are now guaranteed to have a valid hyperparent
+    vtkm::Id hyperparent = augmentedTreeHyperparentsPortal.Get(supernode);
+    vtkm::Id hyperparentSuperId = augmentedTreeHypernodesPortal.Get(hyperparent);
 
+    // we could be at the end of the array, so test explicitly
+    if (supernode == this->AugmentedTreeNumSupernodes - 1)
+    {
+      // this means that we are the end of the segment and can subtract the hyperparent's super ID to get the number of superchildren
+      augmentedTreeSuperchildrenPortal.Set(hyperparent,
+                                           this->AugmentedTreeNumSupernodes - hyperparentSuperId);
+    }
+    // otherwise, if our hyperparent is different from our neighbor's, we are the end of the segment
+    else if (hyperparent != augmentedTreeHyperparentsPortal.Get(supernode + 1))
+    {
+      // again, subtract to get the number
+      augmentedTreeSuperchildrenPortal.Set(
+        hyperparent, this->AugmentedTreeNumSupernodes + 1 - hyperparentSuperId);
+    }
+    // per supernode
     // In serial this worklet implements the following operation
     /*
-      for (vtkm::Id hypernode = 0; hypernode < augmentedTree->hypernodes.size(); hypernode++)
-      { // per hypernode
-        // retrieve the new super ID
-        vtkm::Id superID = augmentedTree->hypernodes[hypernode];
-        // and the next one over
-        vtkm::Id nextSuperID;
-        if (hypernode == augmentedTree->hypernodes.size() - 1)
-          nextSuperID = augmentedTree->supernodes.size();
-        else
-          nextSuperID = augmentedTree->hypernodes[hypernode+1];
-        // the difference is the number of superchildren
-        augmentedTree->superchildren[hypernode] = nextSuperID - superID;
-      } // per hypernode
+      for (indexType supernode = augmentedTree->firstSupernodePerIteration[roundNo][0]; supernode < augmentedTree->firstSupernodePerIteration[roundNo][augmentedTree->nIterations[roundNo]]; supernode++)
+      { // per supernode
+        // attachment points have NULL superarcs and are skipped
+       if (noSuchElement(augmentedTree->superarcs[supernode]))
+          continue;
+       // we are now guaranteed to have a valid hyperparent
+       indexType hyperparent = augmentedTree->hyperparents[supernode];
+       indexType hyperparentSuperID = augmentedTree->hypernodes[hyperparent];
 
+       // we could be at the end of the array, so test explicitly
+      if (supernode == augmentedTree->supernodes.size() - 1)
+        // this means that we are the end of the segment and can subtract the hyperparent's super ID to get the number of superchildren
+        augmentedTree->superchildren[hyperparent] = augmentedTree->supernodes.size() - hyperparentSuperID;
+      // otherwise, if our hyperparent is different from our neighbor's, we are the end of the segment
+      else if (hyperparent != augmentedTree->hyperparents[supernode + 1])
+        // again, subtract to get the number
+        augmentedTree->superchildren[hyperparent] = supernode + 1 - hyperparentSuperID;
+      } // per supernode
     */
   } // operator()()
 

vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_hyper_sweeper/CMakeLists.txt

@@ -17,6 +17,8 @@ set(headers
   TransferTargetComperator.h
   TransferWeightsUpdateRHEWorklet.h
   TransferWeightsUpdateLHEWorklet.h
+  TransferWeightsUpdateRHEWorkletRound2.h
+  TransferWeightsUpdateLHEWorkletRound2.h
 )
 
 vtkm_declare_headers(${headers})

vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_hyper_sweeper/ComputeSuperarcTransferWeightsWorklet.h

@@ -116,7 +116,8 @@ public:
       else
       { // attachment point
         // set the transfer target
-        transferTarget = hierarchicalTreeSuperparentsPortal.Get(supernodeRegularId);
+        transferTarget = hierarchicalTreeSuperparentsPortal.Get(supernodeRegularId) |
+          vtkm::worklet::contourtree_augmented::TRANSFER_TO_SUPERARC;
       } // attachment point
     }   // null superarc
     else

vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_hyper_sweeper/TransferWeightsUpdateLHEWorklet.h

@@ -74,10 +74,11 @@ public:
                                 FieldIn sortedTransferTargetShiftedView,
                                 FieldIn valuePrefixSumShiftedView,
                                 WholeArrayInOut dependentValuesPortal);
-  using ExecutionSgnature = void(_1, _2, _3, _4);
+  using ExecutionSignature = void(InputIndex, _1, _2, _3, _4);
 
   template <typename InOutPortalType>
-  VTKM_EXEC void operator()(const vtkm::Id& sortedTransferTargetValue,
+  VTKM_EXEC void operator()(const vtkm::Id& supernode,
+                            const vtkm::Id& sortedTransferTargetValue,
                             const vtkm::Id& sortedTransferTargetPreviousValue,
                             const vtkm::Id& valuePrefixSumPreviousValue,
                             InOutPortalType& dependentValuesPortal) const
@@ -88,30 +89,33 @@ public:
     {
       return;
     }
-    if (sortedTransferTargetValue != sortedTransferTargetPreviousValue)
+
+    // we need to separate out the flag for attachment points
+    bool superarcTransfer =
+      vtkm::worklet::contourtree_augmented::TransferToSuperarc(sortedTransferTargetValue);
+    vtkm::Id superarcOrNodeId =
+      vtkm::worklet::contourtree_augmented::MaskedIndex(sortedTransferTargetValue);
+
+    // ignore the transfers for attachment points
+    if (superarcTransfer)
     {
-      auto originalValue = dependentValuesPortal.Get(sortedTransferTargetValue);
+      return;
+    }
+
+    // the LHE at 0 is special - it subtracts zero.  In practice, since NO_SUCH_ELEMENT will sort low, this will never
+    // occur, but let's keep the logic strict
+    auto originalValue = dependentValuesPortal.Get(superarcOrNodeId);
+    if (supernode ==
+        0) // i.e., supernode == firstSupernode but we view the arrays when we call the worklet
+    {      // LHE 0
+      dependentValuesPortal.Set(superarcOrNodeId, originalValue - 0);
+    }
+    else if (sortedTransferTargetValue != sortedTransferTargetPreviousValue)
+    { // LHE not 0
+
       dependentValuesPortal.Set(sortedTransferTargetValue,
                                 originalValue - valuePrefixSumPreviousValue);
     }
-
-
-    // In serial this worklet implements the following operation
-    /*
-    for (vtkm::Id supernode = firstSupernode + 1; supernode < lastSupernode; supernode++)
-    { // per supernode
-      // ignore any that point at NO_SUCH_ELEMENT
-      if (noSuchElement(sortedTransferTarget[supernode]))
-        continue;
-
-      // the LHE at 0 is special - it subtracts zero.  In practice, since NO_SUCH_ELEMENT will sort low, this will never
-      // occur, but let's keep the logic strict
-      if (sortedTransferTarget[supernode] != sortedTransferTarget[supernode-1])
-      { // LHE not 0
-        dependentValues[sortedTransferTarget[supernode]] -= valuePrefixSum[supernode-1];
-      } // LHE not 0
-    } // per supernode
-    */
   } // operator()()
 };  // TransferWeightsUpdateLHEWorklet
 

vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_hyper_sweeper/TransferWeightsUpdateLHEWorkletRound2.h

@@ -0,0 +1,132 @@
+//============================================================================
+//  Copyright (c) Kitware, Inc.
+//  All rights reserved.
+//  See LICENSE.txt for details.
+//
+//  This software is distributed WITHOUT ANY WARRANTY; without even
+//  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+//  PURPOSE.  See the above copyright notice for more information.
+//============================================================================
+// Copyright (c) 2018, The Regents of the University of California, through
+// Lawrence Berkeley National Laboratory (subject to receipt of any required approvals
+// from the U.S. Dept. of Energy).  All rights reserved.
+//
+// Redistribution and use in source and binary forms, with or without modification,
+// are permitted provided that the following conditions are met:
+//
+// (1) Redistributions of source code must retain the above copyright notice, this
+//     list of conditions and the following disclaimer.
+//
+// (2) Redistributions in binary form must reproduce the above copyright notice,
+//     this list of conditions and the following disclaimer in the documentation
+//     and/or other materials provided with the distribution.
+//
+// (3) Neither the name of the University of California, Lawrence Berkeley National
+//     Laboratory, U.S. Dept. of Energy nor the names of its contributors may be
+//     used to endorse or promote products derived from this software without
+//     specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+// ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+// WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
+// IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
+// INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
+// BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+// DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
+// OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
+// OF THE POSSIBILITY OF SUCH DAMAGE.
+//
+//=============================================================================
+//
+//  This code is an extension of the algorithm presented in the paper:
+//  Parallel Peak Pruning for Scalable SMP Contour Tree Computation.
+//  Hamish Carr, Gunther Weber, Christopher Sewell, and James Ahrens.
+//  Proceedings of the IEEE Symposium on Large Data Analysis and Visualization
+//  (LDAV), October 2016, Baltimore, Maryland.
+//
+//  The PPP2 algorithm and software were jointly developed by
+//  Hamish Carr (University of Leeds), Gunther H. Weber (LBNL), and
+//  Oliver Ruebel (LBNL)
+//==============================================================================
+
+#ifndef vtk_m_worklet_contourtree_distributed_hierarchical_hyper_sweeper_transfer_weights_update_lhe_worklet_round2_h
+#define vtk_m_worklet_contourtree_distributed_hierarchical_hyper_sweeper_transfer_weights_update_lhe_worklet_round2_h
+
+#include <vtkm/filter/scalar_topology/worklet/contourtree_augmented/Types.h>
+#include <vtkm/worklet/WorkletMapField.h>
+
+namespace vtkm
+{
+namespace worklet
+{
+namespace contourtree_distributed
+{
+namespace hierarchical_hyper_sweeper
+{
+
+/// Worklet used in HierarchicalHyperSweeper.TransferWeights(...) to implement
+/// step 7b. Now find the LHE of each group and subtract out the prior weight
+class TransferWeightsUpdateLHEWorkletRound2 : public vtkm::worklet::WorkletMapField
+{
+public:
+  using ControlSignature = void(FieldIn sortedTransferTargetPortal,
+                                FieldIn sortedTransferTargetShiftedView,
+                                FieldIn valuePrefixSumShiftedView,
+                                WholeArrayInOut intrinsicValues,
+                                WholeArrayInOut dependentValues);
+  using ExecutionSignature = void(InputIndex, _1, _2, _3, _4, _5);
+
+  template <typename InOutPortalType>
+  VTKM_EXEC void operator()(const vtkm::Id& supernode,
+                            const vtkm::Id& sortedTransferTargetValue,
+                            const vtkm::Id& sortedTransferTargetPreviousValue,
+                            const vtkm::Id& valuePrefixSumPreviousValue,
+                            InOutPortalType& intrinsicValuesPortal,
+                            InOutPortalType& dependentValuesPortal) const
+  {
+    // per supernode
+    // ignore any that point at NO_SUCH_ELEMENT
+    if (vtkm::worklet::contourtree_augmented::NoSuchElement(sortedTransferTargetValue))
+    {
+      return;
+    }
+
+    // we need to separate out the flag for attachment points
+    bool superarcTransfer =
+      vtkm::worklet::contourtree_augmented::TransferToSuperarc(sortedTransferTargetValue);
+    vtkm::Id superarcOrNodeId =
+      vtkm::worklet::contourtree_augmented::MaskedIndex(sortedTransferTargetValue);
+
+    // ignore the transfers for attachment points
+    if (superarcTransfer)
+    {
+      return;
+    }
+
+    // the LHE at 0 is special - it subtracts zero.  In practice, since NO_SUCH_ELEMENT will sort low, this will never
+    // occur, but let's keep the logic strict
+    auto originalIntrinsicValue = dependentValuesPortal.Get(superarcOrNodeId);
+    auto originalDependentValue = dependentValuesPortal.Get(superarcOrNodeId);
+    if (supernode ==
+        0) // i.e., supernode == firstSupernode but we view the arrays when we call the worklet
+    {      // LHE 0
+      intrinsicValuesPortal.Set(superarcOrNodeId, originalIntrinsicValue - 0);
+      dependentValuesPortal.Set(superarcOrNodeId, originalDependentValue - 0);
+    }
+    else if (sortedTransferTargetValue != sortedTransferTargetPreviousValue)
+    { // LHE not 0
+      intrinsicValuesPortal.Set(sortedTransferTargetValue,
+                                originalIntrinsicValue - valuePrefixSumPreviousValue);
+      dependentValuesPortal.Set(sortedTransferTargetValue,
+                                originalDependentValue - valuePrefixSumPreviousValue);
+    }
+  } // operator()()
+};  // TransferWeightsUpdateLHEWorklet
+
+} // namespace hierarchical_hyper_sweeper
+} // namespace contourtree_distributed
+} // namespace worklet
+} // namespace vtkm
+
+#endif

vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_hyper_sweeper/TransferWeightsUpdateRHEWorklet.h

@@ -100,8 +100,18 @@ public:
       if ((supernode == this->LastSupernode - 1) ||
           (transferTarget != sortedTransferTargetPortal.Get(supernode + 1)))
       { // RHE of segment
-        auto originalValue = dependentValuesPortal.Get(transferTarget);
-        dependentValuesPortal.Set(transferTarget, originalValue + valuePrefixSum);
+        // we need to separate out the flag for attachment points
+        bool superarcTransfer =
+          vtkm::worklet::contourtree_augmented::TransferToSuperarc(transferTarget);
+        vtkm::Id superarcOrNodeId =
+          vtkm::worklet::contourtree_augmented::MaskedIndex(transferTarget);
+        // we ignore attachment points
+        if (superarcTransfer)
+        {
+          return;
+        }
+        auto originalValue = dependentValuesPortal.Get(superarcOrNodeId);
+        dependentValuesPortal.Set(superarcOrNodeId, originalValue + valuePrefixSum);
       } // RHE of segment
     }
 
@@ -113,10 +123,19 @@ public:
         if (noSuchElement(sortedTransferTarget[supernode]))
           continue;
 
-        // the RHE of each segment transfers its weight (including all irrelevant prefixes)
         if ((supernode == lastSupernode - 1) || (sortedTransferTarget[supernode] != sortedTransferTarget[supernode+1]))
         { // RHE of segment
-          dependentValues[sortedTransferTarget[supernode]] += valuePrefixSum[supernode];
+            // TODO / WARNING 11/07/2023
+            // we need to separate out the flag for attachment points
+            bool superarcTransfer = transferToSuperarc(sortedTransferTarget[supernode]);
+            indexType superarcOrNodeID = maskedIndex(sortedTransferTarget[supernode]);
+
+            // we ignore attachment points
+            if (superarcTransfer)
+               continue;
+
+            // transfer as dependent weight
+            dependentValues[superarcOrNodeID] += valuePrefixSum[supernode];
         } // RHE of segment
       } // per supernode
     */

vtkm/filter/scalar_topology/worklet/contourtree_distributed/hierarchical_hyper_sweeper/TransferWeightsUpdateRHEWorkletRound2.h

@@ -0,0 +1,161 @@
+//============================================================================
+//  Copyright (c) Kitware, Inc.
+//  All rights reserved.
+//  See LICENSE.txt for details.
+//
+//  This software is distributed WITHOUT ANY WARRANTY; without even
+//  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+//  PURPOSE.  See the above copyright notice for more information.
+//============================================================================
+// Copyright (c) 2018, The Regents of the University of California, through
+// Lawrence Berkeley National Laboratory (subject to receipt of any required approvals
+// from the U.S. Dept. of Energy).  All rights reserved.
+//
+// Redistribution and use in source and binary forms, with or without modification,
+// are permitted provided that the following conditions are met:
+//
+// (1) Redistributions of source code must retain the above copyright notice, this
+//     list of conditions and the following disclaimer.
+//
+// (2) Redistributions in binary form must reproduce the above copyright notice,
+//     this list of conditions and the following disclaimer in the documentation
+//     and/or other materials provided with the distribution.
+//
+// (3) Neither the name of the University of California, Lawrence Berkeley National
+//     Laboratory, U.S. Dept. of Energy nor the names of its contributors may be
+//     used to endorse or promote products derived from this software without
+//     specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+// ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+// WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
+// IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
+// INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
+// BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+// DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
+// OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
+// OF THE POSSIBILITY OF SUCH DAMAGE.
+//
+//=============================================================================
+//
+//  This code is an extension of the algorithm presented in the paper:
+//  Parallel Peak Pruning for Scalable SMP Contour Tree Computation.
+//  Hamish Carr, Gunther Weber, Christopher Sewell, and James Ahrens.
+//  Proceedings of the IEEE Symposium on Large Data Analysis and Visualization
+//  (LDAV), October 2016, Baltimore, Maryland.
+//
+//  The PPP2 algorithm and software were jointly developed by
+//  Hamish Carr (University of Leeds), Gunther H. Weber (LBNL), and
+//  Oliver Ruebel (LBNL)
+//==============================================================================
+
+#ifndef vtk_m_worklet_contourtree_distributed_hierarchical_hyper_sweeper_transfer_weights_update_rhe_worklet_round2_h
+#define vtk_m_worklet_contourtree_distributed_hierarchical_hyper_sweeper_transfer_weights_update_rhe_worklet_round2_h
+
+#include <vtkm/filter/scalar_topology/worklet/contourtree_augmented/Types.h>
+#include <vtkm/worklet/WorkletMapField.h>
+
+namespace vtkm
+{
+namespace worklet
+{
+namespace contourtree_distributed
+{
+namespace hierarchical_hyper_sweeper
+{
+
+/// Worklet used in HierarchicalHyperSweeper.TransferWeights(...) to implement
+/// step 7a. Find the RHE of each group and transfer the prefix sum weight.
+/// Note that we do not compute the transfer weight separately, we add it in place instead
+class TransferWeightsUpdateRHEWorkletRound2 : public vtkm::worklet::WorkletMapField
+{
+public:
+  using ControlSignature =
+    void(FieldIn supernodeIndex, // input counting array [firstSupernode, lastSupernode)
+         WholeArrayIn sortedTransferTarget,
+         FieldIn valuePrefixSumView, // input view of valuePrefixSum[firstSupernode, lastSupernode)
+         WholeArrayInOut intrinsicValuesPortal,
+         WholeArrayInOut dependentValuesPortal);
+  using ExecutionSignature = void(_1, _2, _3, _4, _5);
+
+  // Default Constructor
+  VTKM_EXEC_CONT
+  TransferWeightsUpdateRHEWorkletRound2(const vtkm::Id& lastSupernode)
+    : LastSupernode(lastSupernode)
+  {
+  }
+
+  template <typename InPortalType, typename OutPortalType>
+  VTKM_EXEC void operator()(const vtkm::Id& supernode,
+                            const InPortalType& sortedTransferTargetPortal,
+                            const vtkm::Id& valuePrefixSum, // same as valuePrefixSum[supernode]
+                            OutPortalType& intrinsicValuesPortal,
+                            OutPortalType& dependentValuesPortal) const
+  {
+    // per supernode
+    // ignore any that point at NO_SUCH_ELEMENT
+    vtkm::Id transferTarget = sortedTransferTargetPortal.Get(supernode);
+    if (!vtkm::worklet::contourtree_augmented::NoSuchElement(transferTarget))
+    {
+      // the RHE of each segment transfers its weight (including all irrelevant prefixes)
+      if ((supernode == this->LastSupernode - 1) ||
+          (transferTarget != sortedTransferTargetPortal.Get(supernode + 1)))
+      { // RHE of segment
+        // we need to separate out the flag for attachment points
+        bool superarcTransfer =
+          vtkm::worklet::contourtree_augmented::TransferToSuperarc(transferTarget);
+        vtkm::Id superarcOrNodeId =
+          vtkm::worklet::contourtree_augmented::MaskedIndex(transferTarget);
+        // we ignore attachment points
+        if (superarcTransfer)
+        {
+          return;
+        }
+        // we modify both intrinsic and dependent values
+        // we modify both intrinsic and dependent values
+        auto originalIntrinsicValue = intrinsicValuesPortal.Get(superarcOrNodeId);
+        intrinsicValuesPortal.Set(superarcOrNodeId, originalIntrinsicValue + valuePrefixSum);
+        auto originalDependentValue = dependentValuesPortal.Get(superarcOrNodeId);
+        dependentValuesPortal.Set(superarcOrNodeId, originalDependentValue + valuePrefixSum);
+      } // RHE of segment
+    }
+
+    // In serial this worklet implements the following operation
+    /*
+    for (indexType supernode = firstSupernode; supernode < lastSupernode; supernode++)
+    { // per supernode
+      // ignore any that point at NO_SUCH_ELEMENT
+      if (noSuchElement(sortedTransferTarget[supernode]))
+        continue;
+
+      // the RHE of each segment transfers its weight (including all irrelevant prefixes)
+      if ((supernode == lastSupernode - 1) || (sortedTransferTarget[supernode] != sortedTransferTarget[supernode+1]))
+      { // RHE of segment
+        // we need to separate out the flag for attachment points
+        bool superarcTransfer = transferToSuperarc(sortedTransferTarget[supernode]);
+        indexType superarcOrNodeID = maskedIndex(sortedTransferTarget[supernode]);
+        // ignore the transfers for non-attachment points
+        if (!superarcTransfer)
+          continue;
+
+        // we modify both intrinsic and dependent values
+        intrinsicValues[superarcOrNodeID] += valuePrefixSum[supernode];
+        dependentValues[superarcOrNodeID] += valuePrefixSum[supernode];
+      } // RHE of segment
+
+    } // per supernode
+    */
+  } // operator()()
+
+private:
+  const vtkm::Id LastSupernode;
+
+}; // TransferWeightsUpdateRHEWorklet
+
+} // namespace hierarchical_hyper_sweeper
+} // namespace contourtree_distributed
+} // namespace worklet
+} // namespace vtkm
+
+#endif