Merge branch 'add/presimplification' into 'master'

Draft: Add pre-simplification to address scaling performance issues

See merge request vtk/vtk-m!3227

.gitlab/ci/ubuntu1804.yml

@@ -261,6 +261,7 @@ test:ubuntu1804_cuda_perftest:
   variables:
     TEST_INCLUSIONS: "PerformanceTest"
     VTKm_PERF_REMOTE_URL: "https://vbolea:$VTKM_BENCH_RECORDS_TOKEN@gitlab.kitware.com/vbolea/vtk-m-benchmark-records.git"
+    VTKm_PERF_BENCH_DEVICE: "cuda"
     VTKm_PERF_ALPHA: "0.05"
     VTKm_PERF_REPETITIONS: "10"
     VTKm_PERF_DIST: "t"

docs/users-guide/acknowledgements.rst

@@ -90,5 +90,8 @@ Sandia National Laboratories is a multimission laboratory managed and operated b
 This research was supported by the Exascale Computing Project (17-SC-20-SC), a joint project of the U.S.
 Department of Energy's Office of Science and National Nuclear Security Administration, responsible for delivering a capable exascale ecosystem, including software, applications, and hardware technology, to support the nation's exascale computing imperative.
 
-This material is based upon work supported by the U.S.
-Department of Energy, Office of Science, Office of Advanced Scientific Computing Research, Scientific Discovery through Advanced Computing (SciDAC) program.
+This work was supported in part by the U.S. Department of Energy (DOE) RAPIDS SciDAC project under contract number DE-AC05-00OR22725 and by the Exascale Computing Project (17-SC-20-SC), a collaborative effort of the U.S. Department of Energy Office of Science and the National Nuclear Security Administration.
+
+This research used resources of the Oak Ridge Leadership Computing Facility at the Oak Ridge National Laboratory, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC05-00OR22725.
+
+This research used resources of the Argonne Leadership Computing Facility, a U.S. Department of Energy (DOE) Office of Science user facility at Argonne National Laboratory and is based on research supported by the U.S. DOE Office of Science-Advanced Scientific Computing Research Program, under Contract No. DE-AC02-06CH11357.

docs/users-guide/advanced-types.rst

@@ -636,9 +636,7 @@ This is most typically used with C++ run-time type information to convert a run-
 
 .. doxygenfunction:: vtkm::ListForEach(Functor &&f, vtkm::List<Ts...>, Args&&... args)
 
-The following example shows a rudimentary version of converting a dynamically-typed array to a statically-typed array similar to what is done in |VTKm| classes like :class:`vtkm::cont::UnknownArrayHandle` (which is documented in Chapter~\ref{chap:UnknownArrayHandle}).
-
-.. todo:: Fix ``UnknownArrayHandle`` chapter reference above.
+The following example shows a rudimentary version of converting a dynamically-typed array to a statically-typed array similar to what is done in |VTKm| classes like :class:`vtkm::cont::UnknownArrayHandle`, which is documented in :chapref:`unknown-array-handle:Unknown Array Handles`.
 
 .. load-example:: ListForEach
    :file: GuideExampleLists.cxx

docs/users-guide/basic-filter-impl.rst

@@ -72,9 +72,7 @@ The filter implementation can get the appropriate field to operate on using the
 
 One of the challenges with writing filters is determining the actual types the algorithm is operating on.
 The :class:`vtkm::cont::Field` object pulled from the input :class:`vtkm::cont::DataSet` contains a :class:`vtkm::cont::ArrayHandle` (see :chapref:`basic-array-handles:Basic Array Handles`), but you do not know what the template parameters of the :class:`vtkm::cont::ArrayHandle` are.
-There are numerous ways to extract an array of an unknown type out of a :class:`vtkm::cont::ArrayHandle` (many of which will be explored later in Chapter \ref{chap:UnknownArrayHandle}), but the :class:`vtkm::filter::Filter` contains some convenience functions to simplify this.
-
-.. todo:: Fix above reference to unknown array handle chapter.
+There are numerous ways to extract an array of an unknown type out of a :class:`vtkm::cont::ArrayHandle`, many of which will be explored later in :chapref:`unknown-array-handle:Unknown Array Handles`, but the :class:`vtkm::filter::Filter` contains some convenience functions to simplify this.
 
 In particular, this filter operates specifically on scalar fields.
 For this purpose, :class:`vtkm::filter::Filter` provides the :func:`vtkm::filter::Filter::CastAndCallScalarField` helper method.

docs/users-guide/conf.py

@@ -20,7 +20,7 @@ project = "The VTK-m User's Guide"
 copyright = 'Kitware Inc., National Technology & Engineering Solutions of Sandia LLC, UT-Battelle LLC, Los Alamos National Security LLC'
 author = 'Kenneth Moreland'
 version = '@VTKm_VERSION_FULL@'
-release = '@VTKm_VERSION_FULL@'
+release = '@VTKm_VERSION_MAJOR@.@VTKm_VERSION_MINOR@'
 
 # -- General configuration ---------------------------------------------------
 # https://www.sphinx-doc.org/en/master/usage/configuration.html#general-configuration
@@ -76,8 +76,8 @@ today_fmt = '%B %d, %Y'
 rst_prolog = '''
 .. |VTKm| replace:: VTK‑m
 .. |Veclike| replace:: ``Vec``-like
-.. |report-year| replace:: 2023
-.. |report-number| replace:: ORNL/TM-2023/3182
+.. |report-year| replace:: 2024
+.. |report-number| replace:: ORNL/TM-2024/3443
 '''
 
 breathe_projects = { 'vtkm': '@doxygen_xml_output_dir@' }

docs/users-guide/examples/GuideExampleUnknownArrayHandle.cxx

@@ -345,7 +345,7 @@ void TryPrintArrayContents()
   ////
   //// BEGIN-EXAMPLE CastAndCallWithFloatFallback
   ////
-  uncertainArray.CastAndCall(PrintArrayContentsFunctor{});
+  uncertainArray.CastAndCallWithFloatFallback(PrintArrayContentsFunctor{});
   ////
   //// END-EXAMPLE CastAndCallWithFloatFallback
   ////

docs/users-guide/extended-filter-impl.rst

@@ -88,7 +88,7 @@ That calls a lambda function that invokes a worklet to create the output field.
 .. didyouknow::
    The filter implemented in :numref:`ex:FilterFieldImpl` is limited to only find the magnitude of :class:`vtkm::Vec`'s with 3 components.
    It may be the case you wish to implement a filter that operates on :class:`vtkm::Vec`'s of multiple sizes (or perhaps even any size).
-   Chapter \ref{chap:UnknownArrayHandle} discusses how you can use the :class:`vtkm::cont::UnknownArrayHandle` contained in the :class:`vtkm::cont::Field` to more expressively decide what types to check for.
+   :chapref:`unknown-array-handle:Unknown Array Handles` discusses how you can use the :class:`vtkm::cont::UnknownArrayHandle` contained in the :class:`vtkm::cont::Field` to more expressively decide what types to check for.
 
 .. doxygenfunction:: vtkm::filter::Filter::CastAndCallVariableVecField(const vtkm::cont::UnknownArrayHandle&, Functor&&, Args&&...) const
 .. doxygenfunction:: vtkm::filter::Filter::CastAndCallVariableVecField(const vtkm::cont::Field&, Functor&&, Args&&...) const

docs/users-guide/fancy-array-handles.rst

@@ -2,6 +2,8 @@
 Fancy Array Handles
 ==============================
 
+.. todo:: Document :class:`vtkm::cont::ArrayHandleMultiplexer`.
+
 .. index::
    double: array handle; fancy
 

docs/users-guide/index.rst

@@ -23,9 +23,9 @@ David Pugmire,
 Nick Thompson,
 Allison Vacanti,
 Abhishek Yenpure,
-and the |VTKm| community
+and the |VTKm| community.
 
-Moreland, K. (|report-year|). *The VTK-m User's Guide*, Tech report |report-number|, Oak Ridge National Laboratory.
+Moreland, K. (|report-year|). *The VTK-m User's Guide*, version |release|, Tech report |report-number|, Oak Ridge National Laboratory.
 
 .. centered:: Join the VTK-m Community at http://m.vtk.org.
 

docs/users-guide/memory-layout.rst

@@ -83,8 +83,6 @@ Additionally, you can use its constructors or the :func:`vtkm::cont::make_ArrayH
 Strided Arrays
 --------------------
 
-.. todo:: Should this be moved to the chapter/section on transformed arrays?
-
 .. index::
    double: array handle; stride
    double: array handle; offset
@@ -163,3 +161,14 @@ This is convenient for operations that want to operate on arrays with an unknown
 .. load-example:: GetRuntimeVec
    :file: GuideExampleArrayHandleRuntimeVec.cxx
    :caption: Using :class:`vtkm::cont::ArrayHandleRuntimeVec` to get an array regardless of the size of the contained :class:`vtkm::Vec` values.
+
+
+---------------------------------------------
+Recombined Vec Arrays of Strided Components
+---------------------------------------------
+
+|VTKm| contains a special array, :class:`vtkm::cont::ArrayHandleRecombineVec`, to combine component arrays represented in :class:`vtkm::cont::ArrayHandleStride` together to form `Vec` values.
+:class:`vtkm::cont::ArrayHandleRecombineVec` is similar to :class:`vtkm::cont::ArrayHandleSOA` (see :secref:`memory-layout:Structure of Arrays`) except that (1) it holds stride arrays for its components instead of basic arrays and that (2) the number of components can be specified at runtime.
+:class:`vtkm::cont::ArrayHandleRecombineVec` is mainly provided for the implementation of extracting arrays out of a :class:`vtkm::cont::UnknownArrayHandle` (see :secref:`unknown-array-handle:Extracting All Components`).
+
+.. doxygenclass:: vtkm::cont::ArrayHandleRecombineVec

docs/users-guide/part-developing.rst

@@ -9,3 +9,4 @@ Developing Algorithms
    basic-array-handles.rst
    simple-worklets.rst
    basic-filter-impl.rst
+   unknown-array-handle.rst

docs/users-guide/unknown-array-handle.rst

@@ -0,0 +1,457 @@
+==============================
+Unknown Array Handles
+==============================
+
+.. index::
+   single: unknown array handle
+   single: array handle; unknown
+
+The :class:`vtkm::cont::ArrayHandle` class uses templating to make very efficient and type-safe access to data.
+However, it is sometimes inconvenient or impossible to specify the element type and storage at run-time.
+The :class:`vtkm::cont::UnknownArrayHandle` class provides a mechanism to manage arrays of data with unspecified types.
+
+:class:`vtkm::cont::UnknownArrayHandle` holds a reference to an array.
+Unlike :class:`vtkm::cont::ArrayHandle`, :class:`vtkm::cont::UnknownArrayHandle` is *not* templated.
+Instead, it uses C++ run-type type information to store the array without type and cast it when appropriate.
+
+.. doxygenclass:: vtkm::cont::UnknownArrayHandle
+
+.. index:: unknown array handle; construct
+
+An :class:`vtkm::cont::UnknownArrayHandle` can be established by constructing it with or assigning it to an :class:`vtkm::cont::ArrayHandle`.
+The following example demonstrates how an :class:`vtkm::cont::UnknownArrayHandle` might be used to load an array whose type is not known until run-time.
+
+.. load-example:: CreateUnknownArrayHandle
+   :file: GuideExampleUnknownArrayHandle.cxx
+   :caption: Creating an :class:`vtkm::cont::UnknownArrayHandle`.
+
+It is possible to construct a :class:`vtkm::cont::UnknownArrayHandle` that does not point to any :class:`vtkm::cont::ArrayHandle`.
+In this case, the :class:`vtkm::cont::UnknownArrayHandle` is considered not "valid."
+Validity can be tested with the :func:`vtkm::cont::UnknownArrayHandle::IsValid` method.
+
+.. doxygenfunction:: vtkm::cont::UnknownArrayHandle::IsValid
+
+Most of the following operations on :class:`vtkm::cont::UnknownArrayHandle` will fail by throwing an exception if it is not valid.
+Note that it is also possible for a :class:`vtkm::cont::UnknownArrayHandle` to contain an empty :class:`vtkm::cont::ArrayHandle`.
+A :class:`vtkm::cont::UnknownArrayHandle` that contains a :class:`vtkm::cont::ArrayHandle` but has no memory allocated is still considered valid.
+
+Some basic, human-readable information can be retrieved using the :func:`vtkm::cont::UnknownArrayHandle::PrintSummary` method.
+It will print the type and size of the array along with some or all of the values.
+
+.. doxygenfunction:: vtkm::cont::UnknownArrayHandle::PrintSummary
+
+
+------------------------------
+Allocation
+------------------------------
+
+.. index:: unknown array handle; allocation
+
+Data pointed to by an :class:`vtkm::cont::UnknownArrayHandle` is not directly accessible.
+However, it is still possible to do some type-agnostic manipulation of the array allocations.
+
+First, it is always possible to call :func:`vtkm::cont::UnknownArrayHandle::GetNumberOfValues` to retrieve the current size of the array.
+It is also possible to call :func:`vtkm::cont::UnknownArrayHandle::Allocate` to change the size of an unknown array.
+:class:`vtkm::cont::UnknownArrayHandle`'s :func:`vtkm::cont::UnknownArrayHandle::Allocate` works exactly the same as the :func:`vtkm::cont::ArrayHandle::Allocate` in the basic :class:`vtkm::cont::ArrayHandle`.
+
+.. doxygenfunction:: vtkm::cont::UnknownArrayHandle::GetNumberOfValues
+.. doxygenfunction:: vtkm::cont::UnknownArrayHandle::Allocate(vtkm::Id, vtkm::CopyFlag, vtkm::cont::Token&) const
+.. doxygenfunction:: vtkm::cont::UnknownArrayHandle::Allocate(vtkm::Id, vtkm::CopyFlag) const
+
+.. load-example:: UnknownArrayHandleResize
+   :file: GuideExampleUnknownArrayHandle.cxx
+   :caption: Checking the size of a :class:`vtkm::cont::ArrayHandle` and resizing it.
+
+It is often the case where you have an :class:`vtkm::cont::UnknownArrayHandle` as the input to an operation and you want to generate an output of the same type.
+To handle this case, use the :func:`vtkm::cont::UnknownArrayHandle::NewInstance` method to create a new array of the same type (without having to determine the type).
+
+.. doxygenfunction:: vtkm::cont::UnknownArrayHandle::NewInstance
+
+.. load-example:: NonTypeUnknownArrayHandleNewInstance
+   :file: GuideExampleUnknownArrayHandle.cxx
+   :caption: Creating a new instance of an unknown array handle.
+
+That said, there are many special array handles described in :chapref:`memory-layout:Memory Layout of Array Handles` and :chapref:`fancy-array-handles:Fancy Array Handles` that either cannot be directly constructed or cannot be used as outputs.
+Thus, if you do not know the storage of the array, the similar array returned by :func:`vtkm::cont::UnknownArrayHandle::NewInstance` could be infeasible for use as an output.
+Thus, :class:`vtkm::cont::UnknownArrayHandle` also contains the :func:`vtkm::cont::UnknownArrayHandle::NewInstanceBasic` method to create a new array with the same value type but using the basic array storage, which can always be resized and written to.
+
+.. doxygenfunction:: vtkm::cont::UnknownArrayHandle::NewInstanceBasic
+
+.. load-example:: UnknownArrayHandleBasicInstance
+   :file: GuideExampleUnknownArrayHandle.cxx
+   :caption: Creating a new basic instance of an unknown array handle.
+
+It is sometimes the case that you need a new array of a similar type, but that type has to hold floating point values.
+For example, if you had an operation that computed a discrete cosine transform on an array, the result would be very inaccurate if stored as integers.
+In this case, you would actually want to store the result in an array of floating point values.
+For this case, you can use the :func:`vtkm::cont::UnknownArrayHandle::NewInstanceFloatBasic` to create a new basic :class:`vtkm::cont::ArrayHandle` with the component type changed to :type:`vtkm::FloatDefault`.
+For example, if the :class:`vtkm::cont::UnknownArrayHandle` stores an :class:`vtkm::cont::ArrayHandle` of type :type:`vtkm::Id`, :func:`vtkm::cont::UnknownArrayHandle::NewInstanceFloatBasic` will create an :class:`vtkm::cont::ArrayHandle` of type :type:`vtkm::FloatDefault`.
+If the :class:`vtkm::cont::UnknownArrayHandle` stores an :class:`vtkm::cont::ArrayHandle` of type :type:`vtkm::Id3`, :func:`vtkm::cont::UnknownArrayHandle::NewInstanceFloatBasic` will create an :class:`vtkm::cont::ArrayHandle` of type :type:`vtkm::Vec3f`.
+
+.. doxygenfunction:: vtkm::cont::UnknownArrayHandle::NewInstanceFloatBasic
+
+.. load-example:: UnknownArrayHandleFloatInstance
+   :file: GuideExampleUnknownArrayHandle.cxx
+   :caption: Creating a new array instance with floating point values.
+
+Finally, it may be the case where you are finished using a :class:`vtkm::cont::UnknownArrayHandle`.
+If you want to free up memory on the device, which may have limited memory, you can do so with :func:`vtkm::cont::UnknownArrayHandle::ReleaseResourcesExecution`, which will free any memory on the device but preserve the data on the host.
+If the data will never be used again, all memory can be freed with :func:`vtkm::cont::UnknownArrayHandle::ReleaseResources`
+
+.. doxygenfunction:: vtkm::cont::UnknownArrayHandle::ReleaseResourcesExecution
+.. doxygenfunction:: vtkm::cont::UnknownArrayHandle::ReleaseResources
+
+
+------------------------------
+Casting to Known Types
+------------------------------
+
+.. index::
+   single: unknown array handle; cast
+   single: unknown array handle; as array handle
+
+Data pointed to by an :class:`vtkm::cont::UnknownArrayHandle` is not directly
+accessible.
+To access the data, you need to retrieve the data as an :class:`vtkm::cont::ArrayHandle`.
+If you happen to know (or can guess) the type, you can use the :func:`vtkm::cont::UnknownArrayHandle::AsArrayHandle` method to retrieve the array as a specific type.
+
+.. doxygenfunction:: vtkm::cont::UnknownArrayHandle::AsArrayHandle(vtkm::cont::ArrayHandle<T, S>&) const
+.. doxygenfunction:: vtkm::cont::UnknownArrayHandle::AsArrayHandle() const
+
+.. load-example:: UnknownArrayHandleAsArrayHandle1
+   :file: GuideExampleUnknownArrayHandle.cxx
+   :caption: Retrieving an array of a known type from :class:`vtkm::cont::UnknownArrayHandle`.
+
+:func:`vtkm::cont::UnknownArrayHandle::AsArrayHandle` actually has two forms.
+The first form, shown in the previous example, has no arguments and returns the :class:`vtkm::cont::ArrayHandle`.
+This form requires you to specify the type of array as a template parameter.
+The alternate form has you pass a reference to a concrete :class:`vtkm::cont::ArrayHandle` as an argument as shown in the following example.
+This form can imply the template parameter from the argument.
+
+.. load-example:: UnknownArrayHandleAsArrayHandle2
+   :file: GuideExampleUnknownArrayHandle.cxx
+   :caption: Alternate form for retrieving an array of a known type from :class:`vtkm::cont::UnknownArrayHandle`.
+
+:func:`vtkm::cont::UnknownArrayHandle::AsArrayHandle` treats :class:`vtkm::cont::ArrayHandleCast` and :class:`vtkm::cont::ArrayHandleMultiplexer` special.
+If the special :class:`vtkm::cont::ArrayHandle` can hold the actual array stored, then :func:`vtkm::cont::UnknownArrayHandle::AsArrayHandle` will return successfully.
+In the following example, :func:`vtkm::cont::UnknownArrayHandle::AsArrayHandle` returns an array of type :type:`vtkm::Float32` as an :class:`vtkm::cont::ArrayHandleCast` that converts the values to :type:`vtkm::Float64`.
+
+.. load-example:: UnknownArrayHandleAsCastArray
+   :file: GuideExampleUnknownArrayHandle.cxx
+   :caption: Getting a cast array handle from an :class:`vtkm::cont::ArrayHandleCast`.
+
+.. didyouknow::
+   The inverse retrieval works as well.
+   If you create an :class:`vtkm::cont::UnknownArrayHandle` with an :class:`vtkm::cont::ArrayHandleCast` or :class:`vtkm::cont::ArrayHandleMultiplexer`, you can get the underlying array with :func:`vtkm::cont::UnknownArrayHandle::AsArrayHandle`.
+   These relationships also work recursively (e.g. an array placed in a cast array that is placed in a multiplexer).
+
+.. index:: unknown array handle; query type
+
+If the :class:`vtkm::cont::UnknownArrayHandle` cannot store its array in the type given to :func:`vtkm::cont::UnknownArrayHandle::AsArrayHandle`, it will throw an exception.
+Thus, you should not use :func:`vtkm::cont::UnknownArrayHandle::AsArrayHandle` with types that you are not sure about.
+Use the :func:`vtkm::cont::UnknownArrayHandle::CanConvert` method to determine if a given :class:`vtkm::cont::ArrayHandle` type will work with :func:`vtkm::cont::UnknownArrayHandle::AsArrayHandle`.
+
+.. doxygenfunction:: vtkm::cont::UnknownArrayHandle::CanConvert
+
+.. load-example:: UnknownArrayHandleCanConvert
+   :file: GuideExampleUnknownArrayHandle.cxx
+   :caption: Querying whether a given :class:`vtkm::cont::ArrayHandle` can be retrieved from a :class:`vtkm::cont::UnknownArrayHandle`.
+
+By design, :func:`vtkm::cont::UnknownArrayHandle::CanConvert` will return true for types that are not actually stored in the :class:`vtkm::cont::UnknownArrayHandle` but can be retrieved.
+If you need to know specifically what type is stored in the :class:`vtkm::cont::UnknownArrayHandle`, you can use the :func:`vtkm::cont::UnknownArrayHandle::IsType` method instead.
+
+.. doxygenfunction:: vtkm::cont::UnknownArrayHandle::IsType
+
+If you need to query either the value type or the storage, you can use :func:`vtkm::cont::UnknownArrayHandle::IsValueType` and :func:`vtkm::cont::UnknownArrayHandle::IsStorageType`, respectively.
+:class:`vtkm::cont::UnknownArrayHandle` also provides :func:`vtkm::cont::UnknownArrayHandle::GetValueTypeName`, :func:`vtkm::cont::UnknownArrayHandle::GetStorageTypeName`, and :func:`vtkm::cont::UnknownArrayHandle::GetArrayTypeName` for debugging purposes.
+
+.. doxygenfunction:: vtkm::cont::UnknownArrayHandle::IsValueType
+.. doxygenfunction:: vtkm::cont::UnknownArrayHandle::IsStorageType
+.. doxygenfunction:: vtkm::cont::UnknownArrayHandle::GetValueTypeName
+.. doxygenfunction:: vtkm::cont::UnknownArrayHandle::GetStorageTypeName
+.. doxygenfunction:: vtkm::cont::UnknownArrayHandle::GetArrayTypeName
+
+.. commonerrors::
+   :func:`vtkm::cont::UnknownArrayHandle::CanConvert` is almost always safer to use than :func:`vtkm::cont::UnknownArrayHandle::IsType` or its similar methods.
+   Even though :func:`vtkm::cont::UnknownArrayHandle::IsType` reflects the actual array type, :func:`vtkm::cont::UnknownArrayHandle::CanConvert` better describes how :class:`vtkm::cont::UnknownArrayHandle` will behave.
+
+If you do not know the exact type of the array contained in an :class:`vtkm::cont::UnknownArrayHandle`, a brute force method to get the data out is to copy it to an array of a known type.
+This can be done with the :func:`vtkm::cont::UnknownArrayHandle::DeepCopyFrom` method, which will copy the contents of a target array into an existing array of a (potentially) different type.
+
+.. doxygenfunction:: vtkm::cont::UnknownArrayHandle::DeepCopyFrom(const vtkm::cont::UnknownArrayHandle&)
+.. doxygenfunction:: vtkm::cont::UnknownArrayHandle::DeepCopyFrom(const vtkm::cont::UnknownArrayHandle&) const
+
+.. load-example:: UnknownArrayHandleDeepCopy
+   :file: GuideExampleUnknownArrayHandle.cxx
+   :caption: Deep copy arrays of unknown types.
+
+It is often the case that you have good reason to believe that an array is of an expected type, but you have no way to be sure.
+To simplify code, the most rational thing to do is to get the array as the expected type if that is indeed what it is, or to copy it to an array of that type otherwise.
+The :func:`vtkm::cont::UnknownArrayHandle::CopyShallowIfPossible` does just that.
+
+.. doxygenfunction:: vtkm::cont::UnknownArrayHandle::CopyShallowIfPossible(const vtkm::cont::UnknownArrayHandle&)
+.. doxygenfunction:: vtkm::cont::UnknownArrayHandle::CopyShallowIfPossible(const vtkm::cont::UnknownArrayHandle&) const
+
+.. load-example:: UnknownArrayHandleShallowCopy
+   :file: GuideExampleUnknownArrayHandle.cxx
+   :caption: Using :func:`vtkm::cont::UnknownArrayHandle::CopyShallowIfPossible` to get an unknown array as a particular type.
+
+.. didyouknow::
+   The :class:`vtkm::cont::UnknownArrayHandle` copy methods behave similarly to the :func:`vtkm::cont::ArrayCopy` functions.
+
+
+----------------------------------------
+Casting to a List of Potential Types
+----------------------------------------
+
+.. index:: unknown array handle; cast
+
+Using :func:`vtkm::cont::UnknownArrayHandle::AsArrayHandle` is fine as long as the correct types are known, but often times they are not.
+For this use case :class:`vtkm::cont::UnknownArrayHandle` has a method named :func:`vtkm::cont::UnknownArrayHandle::CastAndCallForTypes` that attempts to cast the array to some set of types.
+
+The :func:`vtkm::cont::UnknownArrayHandle::CastAndCallForTypes` method accepts a functor to run on the appropriately cast array.
+The functor must have an overloaded const parentheses operator that accepts an :class:`vtkm::cont::ArrayHandle` of the appropriate type.
+You also have to specify two template parameters that specify a :class:`vtkm::List` of value types to try and a :class:`vtkm::List` of storage types to try, respectively.
+The macros :c:macro:`VTKM_DEFAULT_TYPE_LIST` and :c:macro:`VTKM_DEFAULT_STORAGE_LIST` are often used when nothing more specific is known.
+
+.. doxygenfunction:: vtkm::cont::UnknownArrayHandle::CastAndCallForTypes
+
+.. load-example:: UsingCastAndCallForTypes
+   :file: GuideExampleUnknownArrayHandle.cxx
+   :caption: Operating on an :class:`vtkm::cont::UnknownArrayHandle` with :func:`vtkm::cont::UnknownArrayHandle::CastAndCallForTypes`.
+
+.. didyouknow::
+   The first (required) argument to :func:`vtkm::cont::UnknownArrayHandle::CastAndCallForTypes` is the functor to call with the array.
+   You can supply any number of optional arguments after that.
+   Those arguments will be passed directly to the functor.
+   This makes it easy to pass state to the functor.
+
+.. didyouknow::
+   When an :class:`vtkm::cont::UnknownArrayHandle` is used in place of an :class:`vtkm::cont::ArrayHandle` as an argument to a worklet invocation, it will internally use :func:`vtkm::cont::UnknownArrayHandle::CastAndCallForTypes` to attempt to call the worklet with an :class:`vtkm::cont::ArrayHandle` of the correct type.
+
+:class:`vtkm::cont::UnknownArrayHandle` has a simple subclass named :class:`vtkm::cont::UncertainArrayHandle` for use when you can narrow the array to a finite set of types.
+:class:`vtkm::cont::UncertainArrayHandle` has two template parameters that must be specified: a :class:`vtkm::List` of value types and a :class:`vtkm::List` of storage types.
+
+.. doxygenclass:: vtkm::cont::UncertainArrayHandle
+
+:class:`vtkm::cont::UncertainArrayHandle` has a method named :func:`vtkm::cont::UncertainArrayHandle::CastAndCall` that behaves the same as :func:`vtkm::cont::UnknownArrayHandle::CastAndCallForTypes` except that you do not have to specify the types to try.
+Instead, the types are taken from the template parameters of the :class:`vtkm::cont::UncertainArrayHandle` itself.
+
+.. doxygenfunction:: vtkm::cont::UncertainArrayHandle::CastAndCall
+
+.. load-example:: UncertainArrayHandle
+   :file: GuideExampleUnknownArrayHandle.cxx
+   :caption: Using :class:`vtkm::cont::UncertainArrayHandle` to cast and call a functor.
+
+.. didyouknow::
+   Like with :class:`vtkm::cont::UnknownArrayHandle`, if an :class:`vtkm::cont::UncertainArrayHandle` is used in a worklet invocation, it will internally use :func:`vtkm::cont::UncertainArrayHandle::CastAndCall`.
+   This provides a convenient way to specify what array types the invoker should try.
+
+Both :class:`vtkm::cont::UnknownArrayHandle` and :class:`vtkm::cont::UncertainArrayHandle` provide a method named :func:`vtkm::cont::UnknownArrayHandle::ResetTypes` to redefine the types to try.
+:func:`vtkm::cont::UncertainArrayHandle::ResetTypes` has two template parameters that are the :class:`vtkm::List`'s of value and storage types.
+:func:`vtkm::cont::UnknownArrayHandle::ResetTypes` returns a new :class:`vtkm::cont::UncertainArrayHandle` with the given types.
+This is a convenient way to pass these types to functions.
+
+.. doxygenfunction:: vtkm::cont::UnknownArrayHandle::ResetTypes
+
+:class:`vtkm::cont::UncertainArrayHandle` additionally has methods named :func:`vtkm::cont::UncertainArrayHandle::ResetValueTypes` and :func:`vtkm::cont::UncertainArrayHandle::ResetStorageTypes` to reset the value types and storage types, respectively, without modifying the other.
+
+.. doxygenfunction:: vtkm::cont::UncertainArrayHandle::ResetValueTypes
+.. doxygenfunction:: vtkm::cont::UncertainArrayHandle::ResetStorageTypes
+
+.. load-example:: UnknownArrayResetTypes
+   :file: GuideExampleUnknownArrayHandle.cxx
+   :caption: Resetting the types of an :class:`vtkm::cont::UnknownArrayHandle`.
+
+.. commonerrors::
+   Because it returns an :class:`vtkm::cont::UncertainArrayHandle`, you need to include :file:`vtkm/cont/UncertainArrayHandle.h` if you use :func:`vtkm::cont::UnknownArrayHandle::ResetTypes`.
+   This is true even if you do not directly use the returned object.
+
+
+------------------------------
+Accessing Truly Unknown Arrays
+------------------------------
+
+So far in :secref:`unknown-array-handle:Casting to Known Types` and :secref:`unknown-array-handle:Casting to a List of Potential Types` we explored how to access the data in an :class:`vtkm::cont::UnknownArrayHandle` when you actually know the array type or can narrow down the array type to some finite number of candidates.
+But what happens if you cannot practically narrow down the types in the :class:`vtkm::cont::UnknownArrayHandle`?
+For this case, :class:`vtkm::cont::UnknownArrayHandle` provides mechanisms for extracting data knowing little or nothing about the types.
+
+Cast with Floating Point Fallback
+========================================
+
+.. index:: unknown array handle; fallback
+
+The problem with :func:`vtkm::cont::UnknownArrayHandle::CastAndCallForTypes` and :func:`vtkm::cont::UncertainArrayHandle::CastAndCall` is that you can only list a finite amount of value types and storage types to try.
+If you encounter an :class:`vtkm::cont::UnknownArrayHandle` containing a different :class:`vtkm::cont::ArrayHandle` type, the cast and call will simply fail.
+Since the compiler must create a code path for each possible :class:`vtkm::cont::ArrayHandle` type, it may not even be feasible to list all known types.
+
+:func:`vtkm::cont::UnknownArrayHandle::CastAndCallForTypesWithFloatFallback` works around this problem by providing a fallback in case the contained :class:`vtkm::cont::ArrayHandle` does not match any of the types tried.
+If none of the types match, then :func:`vtkm::cont::UnknownArrayHandle::CastAndCallForTypesWithFloatFallback` will copy the data to a :class:`vtkm::cont::ArrayHandle` with :type:`vtkm::FloatDefault` values (or some compatible :class:`vtkm::Vec` with :type:`vtkm::FloatDefault` components) and basic storage.
+It will then attempt to match again with this copied array.
+
+.. doxygenfunction:: vtkm::cont::UnknownArrayHandle::CastAndCallForTypesWithFloatFallback
+
+.. load-example:: CastAndCallForTypesWithFloatFallback
+   :file: GuideExampleUnknownArrayHandle.cxx
+   :caption: Cast and call a functor from an :class:`vtkm::cont::UnknownArrayHandle` with a float fallback.
+
+In this case, we do not have to list every possible type because the array will be copied to a known type if nothing matches.
+Note that when using :func:`vtkm::cont::UnknownArrayHandle::CastAndCallForTypesWithFloatFallback`, you still need to include an appropriate type based on :type:`vtkm::FloatDefault` in the value type list and :class:`vtkm::cont::StorageTagBasic` in the storage list so that the copied array can match.
+
+:class:`vtkm::cont::UncertainArrayHandle` has a matching method named :func:`vtkm::cont::UncertainArrayHandle::CastAndCallWithFloatFallback` that does the same operation using the types specified in the :class:`vtkm::cont::UncertainArrayHandle`.
+
+.. doxygenfunction:: vtkm::cont::UncertainArrayHandle::CastAndCallWithFloatFallback
+
+.. load-example:: CastAndCallWithFloatFallback
+   :file: GuideExampleUnknownArrayHandle.cxx
+   :caption: Cast and call a functor from an :class:`vtkm::cont::UncertainArrayHandle` with a float fallback.
+
+Extracting Components
+==============================
+
+Using a floating point fallback allows you to use arrays of unknown types in most circumstances, but it does have a few drawbacks.
+First, and most obvious, is that you may not operate on the data in its native format.
+If you want to preserve the integer format of data, this may not be the method.
+Second, the fallback requires a copy of the data.
+If :func:`vtkm::cont::UnknownArrayHandle::CastAndCallForTypesWithFloatFallback` does not match the type of the array, it copies the array to a new type that (hopefully) can be matched.
+Third, :func:`vtkm::cont::UnknownArrayHandle::CastAndCallForTypesWithFloatFallback` still needs to match the number of components in each array value.
+If the contained :class:`vtkm::cont::ArrayHandle` contains values that are :class:`vtkm::Vec`'s of length 2, then the data will be copied to an array of :type:`vtkm::Vec2f`'s.
+If :type:`vtkm::Vec2f` is not included in the types to try, the cast and call will still fail.
+
+.. index:: unknown array handle; extract component
+
+A way to get around these problems is to extract a single component from the array.
+You can use the :func:`vtkm::cont::UnknownArrayHandle::ExtractComponent` method to return an :class:`vtkm::cont::ArrayHandle` with the values for a given component for each value in the array.
+The type of the returned :class:`vtkm::cont::ArrayHandle` will be the same regardless of the actual array type stored in the :class:`vtkm::cont::UnknownArrayHandle`.
+
+.. doxygenfunction:: vtkm::cont::UnknownArrayHandle::ExtractComponent
+
+:func:`vtkm::cont::UnknownArrayHandle::ExtractComponent` must be given a template argument for the base component type.
+The following example extracts the first component of all :class:`vtkm::Vec` values in an :class:`vtkm::cont::UnknownArrayHandle` assuming that the component is of type :type:`vtkm::FloatDefault` (:exlineref:`ex:UnknownArrayExtractComponent:Call`).
+
+.. load-example:: UnknownArrayExtractComponent
+   :file: GuideExampleUnknownArrayHandle.cxx
+   :caption: Extracting the first component of every value in an :class:`vtkm::cont::UnknownArrayHandle`.
+
+The code in :numref:`ex:UnknownArrayExtractComponent` works with any array with values based on the default floating point type.
+If the :class:`vtkm::cont::UnknownArrayHandle` has an array containing :type:`vtkm::FloatDefault`, then the returned array has all the same values.
+If the :class:`vtkm::cont::UnknownArrayHandle` contains values of type :type:`vtkm::Vec3f`, then each value in the returned array will be the first component of this array.
+
+If the :class:`vtkm::cont::UnknownArrayHandle` really contains an array with incompatible value types (such as ``vtkm::cont::ArrayHandle<vtkm::Id>``), then an :class:`vtkm::cont::ErrorBadType` will be thrown.
+To check if the :class:`vtkm::cont::UnknownArrayHandle` contains an array of a compatible type, use the :func:`vtkm::cont::UnknownArrayHandle::IsBaseComponentType` method to check the component type being used as the template argument to :func:`vtkm::cont::UnknownArrayHandle::ExtractComponent`.
+
+.. doxygenfunction:: vtkm::cont::UnknownArrayHandle::IsBaseComponentType
+
+.. load-example:: UnknownArrayBaseComponentType
+   :file: GuideExampleUnknownArrayHandle.cxx
+   :caption: Checking the base component type in an :class:`vtkm::cont::UnknownArrayHandle`.
+
+it is also possible to get a name for the base component type (mostly for debugging purposes) with :func:`vtkm::cont::UnknownArrayHandle::GetBaseComponentTypeName`.
+
+.. doxygenfunction:: vtkm::cont::UnknownArrayHandle::GetBaseComponentTypeName
+
+You will often need to query the number of components that can be extracted from the array.
+This can be queried with :func:`vtkm::cont::UnknownArrayHandle::GetNumberOfComponentsFlat`.
+
+.. doxygenfunction:: vtkm::cont::UnknownArrayHandle::GetNumberOfComponentsFlat
+
+This section started with the motivation of getting data from an :class:`vtkm::cont::UnknownArrayHandle` without knowing anything about the type, yet :func:`vtkm::cont::UnknownArrayHandle::ExtractComponent` still requires a type parameter.
+However, by limiting the type needed to the base component type, you only need to check the base C types (standard integers and floating points) available in C++.
+You do not need to know whether these components are arranged in :class:`vtkm::Vec`'s or the size of the :class:`vtkm::Vec`.
+A general implementation of an algorithm might have to deal with scalars as well as :class:`vtkm::Vec`'s of size 2, 3, and 4.
+If we consider operations on tensors, :class:`vtkm::Vec`'s of size 6 and 9 can be common as well.
+But when using :func:`vtkm::cont::UnknownArrayHandle::ExtractComponent`, a single condition can handle any potential :class:`vtkm::Vec` size.
+
+Another advantage of :func:`vtkm::cont::UnknownArrayHandle::ExtractComponent` is that the type of storage does not need to be specified.
+:func:`vtkm::cont::UnknownArrayHandle::ExtractComponent` works with any type of :class:`vtkm::cont::ArrayHandle` storage (with some caveats).
+So, :numref:`ex:UnknownArrayExtractComponent` works equally as well with :class:`vtkm::cont::ArrayHandleBasic`, :class:`vtkm::cont::ArrayHandleSOA`, :class:`vtkm::cont::ArrayHandleUniformPointCoordinates`, :class:`vtkm::cont::ArrayHandleCartesianProduct`, and many others.
+Trying to capture all reasonable types of arrays could easily require hundreds of conditions, all of which and more can be captured with :func:`vtkm::cont::UnknownArrayHandle::ExtractComponent` and the roughly 12 basic C data types.
+In practice, you often only really have to worry about floating point components, which further reduces the cases down to (usually) 2.
+
+:func:`vtkm::cont::UnknownArrayHandle::ExtractComponent` works by returning an :class:`vtkm::cont::ArrayHandleStride`.
+This is a special :class:`vtkm::cont::ArrayHandle` that can access data buffers by skipping values at regular intervals.
+This allows it to access data packed in different ways such as :class:`vtkm::cont::ArrayHandleBasic`, :class:`vtkm::cont::ArrayHandleSOA`, and many others.
+That said, :class:`vtkm::cont::ArrayHandleStride` is not magic, so if it cannot directly access memory, some or all of it may be copied.
+If you are attempting to use the array from :func:`vtkm::cont::UnknownArrayHandle::ExtractComponent` as an output array, pass :enumerator:`vtkm::CopyFlag::Off` as a second argument.
+This will ensure that data are not copied so that any data written will go to the original array (or throw an exception if this cannot be done).
+
+.. commonerrors::
+   Although :func:`vtkm::cont::UnknownArrayHandle::ExtractComponent` will technically work with any :class:`vtkm::cont::ArrayHandle` (of simple :class:`vtkm::Vec` types), it may require a very inefficient memory copy.
+   Pay attention if :func:`vtkm::cont::UnknownArrayHandle::ExtractComponent` issues a warning about an inefficient memory copy.
+   This is likely a serious performance issue, and the data should be retrieved in a different way (or better yet stored in a different way).
+
+Extracting All Components
+==============================
+
+:numref:`ex:UnknownArrayExtractComponent` accesses the first component of each :class:`vtkm::Vec` in an array.
+But in practice you usually want to operate on all components stored in the array.
+A simple solution is to iterate over each component.
+
+.. load-example:: UnknownArrayExtractComponentsMultiple
+   :file: GuideExampleUnknownArrayHandle.cxx
+   :caption: Extracting each component from an :class:`vtkm::cont::UnknownArrayHandle`.
+
+To ensure that the type of the extracted component is a basic C type, the :class:`vtkm::Vec` values are "flattened."
+That is, they are treated as if they are a single level :class:`vtkm::Vec`.
+For example, if you have a value type of ``vtkm::Vec<vtkm::Id3, 2>``, :func:`vtkm::cont::UnknownArrayHandle::ExtractComponent` treats this type as ``vtkm::Vec<vtkm::Id, 6>``.
+This allows you to extract the components as type :type:`vtkm::Id` rather than having a special case for :type:`vtkm::Id3`.
+
+Although iterating over components works fine, it can be inconvenient.
+An alternate mechanism is to use :func:`vtkm::cont::UnknownArrayHandle::ExtractArrayFromComponents` to get all the components at once.
+:func:`vtkm::cont::UnknownArrayHandle::ExtractArrayFromComponents` works like :func:`vtkm::cont::UnknownArrayHandle::ExtractComponent` except that instead of returning an :class:`vtkm::cont::ArrayHandleStride`, it returns a special :class:`vtkm::cont::ArrayHandleRecombineVec` that behaves like an :class:`vtkm::cont::ArrayHandle` to reference all component arrays at once.
+
+.. doxygenfunction:: vtkm::cont::UnknownArrayHandle::ExtractArrayFromComponents
+
+.. load-example:: UnknownArrayExtractArrayFromComponents
+   :file: GuideExampleUnknownArrayHandle.cxx
+   :caption: Extracting all components from an :class:`vtkm::cont::UnknownArrayHandle` at once.
+
+.. commonerrors::
+   Although it has the same interface as other :class:`vtkm::cont::ArrayHandle`'s, :class:`vtkm::cont::ArrayHandleRecombineVec` has a special value type that breaks some conventions.
+   For example, when used in a worklet, the value type passed from this array to the worklet cannot be replicated.
+   That is, you cannot create a temporary stack value of the same type.
+
+Because you still need to specify a base component type, you will likely still need to check several types to safely extract data from an :class:`vtkm::cont::UnknownArrayHandle` by component.
+To do this automatically, you can use the :func:`vtkm::cont::UnknownArrayHandle::CastAndCallWithExtractedArray`.
+This method behaves similarly to :func:`vtkm::cont::UncertainArrayHandle::CastAndCall` except that it internally uses :func:`vtkm::cont::UnknownArrayHandle::ExtractArrayFromComponents`.
+
+.. doxygenfunction:: vtkm::cont::UnknownArrayHandle::CastAndCallWithExtractedArray
+
+.. load-example:: UnknownArrayCallWithExtractedArray
+   :file: GuideExampleUnknownArrayHandle.cxx
+   :caption: Calling a functor for nearly any type of array stored in an :class:`vtkm::cont::UnknownArrayHandle`.
+
+
+------------------------------
+Mutability
+------------------------------
+
+.. index:: unknown array handle; const
+
+One subtle feature of :class:`vtkm::cont::UnknownArrayHandle` is that the class is, in principle, a pointer to an array pointer.
+This means that the data in an :class:`vtkm::cont::UnknownArrayHandle` is always mutable even if the class is declared ``const``.
+The upshot is that you can pass output arrays as constant :class:`vtkm::cont::UnknownArrayHandle` references.
+
+.. load-example:: UnknownArrayConstOutput
+   :file: GuideExampleUnknownArrayHandle.cxx
+   :caption: Using a ``const`` :class:`vtkm::cont::UnknownArrayHandle` for a function output.
+
+Although it seems strange, there is a good reason to allow an output :class:`vtkm::cont::UnknownArrayHandle` to be ``const``.
+It allows a typed :class:`vtkm::cont::ArrayHandle` to be used as the argument to the function.
+In this case, the compiler will automatically convert the :class:`vtkm::cont::ArrayHandle` to a :class:`vtkm::cont::UnknownArrayHandle`.
+When C++ creates objects like this, they can only be passed a constant reference, an Rvalue reference, or by value.
+So, declaring the output parameter as ``const`` :class:`vtkm::cont::UnknownArrayHandle` allows it to be used for code like this.
+
+.. load-example:: UseUnknownArrayConstOutput
+   :file: GuideExampleUnknownArrayHandle.cxx
+   :caption: Passing an :class:`vtkm::cont::ArrayHandle` as an output :class:`vtkm::cont::UnknownArrayHandle`.
+
+Of course, you could also declare the output by value instead of by reference, but this has the same semantics with extra internal pointer management.
+
+.. didyouknow::
+   When possible, it is better to pass a :class:`vtkm::cont::UnknownArrayHandle` as a constant reference (or by value) rather than a mutable reference, even if the array contents are going to be modified.
+   This allows the function to support automatic conversion of an output :class:`vtkm::cont::ArrayHandle`.
+
+So if a constant :class:`vtkm::cont::UnknownArrayHandle` can have its contents modified, what is the difference between a constant reference and a non-constant reference?
+The difference is that the constant reference can change the array's content, but not the array itself.
+If you want to do operations like doing a shallow copy or changing the underlying type of the array, a non-constant reference is needed.

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/ArrayHandleRecombineVec.h

@@ -587,19 +587,19 @@ public:
 
 } // namespace internal
 
-/// \brief A grouping of `ArrayHandleStride`s into an `ArrayHandle` of `Vec`s.
+/// @brief A grouping of `ArrayHandleStride`s into an `ArrayHandle` of `vtkm::Vec`s.
 ///
 /// The main intention of `ArrayHandleStride` is to pull out a component of an
-/// `ArrayHandle` without knowing there `ArrayHandle`'s storage or `Vec` shape.
+/// `ArrayHandle` without knowing there `ArrayHandle`'s storage or `vtkm::Vec` shape.
 /// However, usually you want to do an operation on all the components together.
 /// `ArrayHandleRecombineVec` implements the functionality to easily take a
 /// group of extracted components and treat them as a single `ArrayHandle` of
-/// `Vec` values.
+/// `vtkm::Vec` values.
 ///
 /// Note that caution should be used with `ArrayHandleRecombineVec` because the
-/// size of the `Vec` values is not known at compile time. Thus, the value
+/// size of the `vtkm::Vec` values is not known at compile time. Thus, the value
 /// type of this array is forced to a special `RecombineVec` class that can cause
-/// surprises if treated as a `Vec`. In particular, the static `NUM_COMPONENTS`
+/// surprises if treated as a `vtkm::Vec`. In particular, the static `NUM_COMPONENTS`
 /// expression does not exist. Furthermore, new variables of type `RecombineVec`
 /// cannot be created. This means that simple operators like `+` will not work
 /// because they require an intermediate object to be created. (Equal operators
@@ -618,17 +618,33 @@ public:
     (vtkm::cont::ArrayHandle<internal::detail::RecombinedValueType<ComponentType>,
                              vtkm::cont::internal::StorageTagRecombineVec>));
 
+  /// @brief Return the number of components in each value of the array.
+  ///
+  /// This is also equal to the number of component arrays referenced by this
+  /// fancy array.
+  ///
+  /// `ArrayHandleRecombineVec` always stores flat Vec values. As such, this number
+  /// of components is the same as the number of base components.
   vtkm::IdComponent GetNumberOfComponents() const
   {
     return StorageType::GetNumberOfComponents(this->GetBuffers());
   }
 
+  /// @brief Get the array storing the values for a particular component.
+  ///
+  /// The returned array is a `vtkm::cont::ArrayHandleStride`. It is possible
+  /// that the returned arrays from different components reference the same area
+  /// of physical memory (usually referencing values interleaved with each other).
   vtkm::cont::ArrayHandleStride<ComponentType> GetComponentArray(
     vtkm::IdComponent componentIndex) const
   {
     return StorageType::ArrayForComponent(this->GetBuffers(), componentIndex);
   }
 
+  /// @brief Add a component array.
+  ///
+  /// `AppendComponentArray()` provides an easy way to build an `ArrayHandleRecombineVec`
+  /// by iteratively adding the component arrays.
   void AppendComponentArray(
     const vtkm::cont::ArrayHandle<ComponentType, vtkm::cont::StorageTagStride>& array)
   {

vtkm/cont/UncertainArrayHandle.h

@@ -109,10 +109,10 @@ public:
 
   /// \brief Call a functor using the underlying array type with a float cast fallback.
   ///
-  /// `CastAndCallWithFloatFallback` attempts to cast the held array to a specific value type,
+  /// `CastAndCallWithFloatFallback()` attempts to cast the held array to a specific value type,
   /// and then calls the given functor with the cast array. If the underlying array
   /// does not match any of the requested array types, the array is copied to a new
-  /// `ArrayHandleBasic` with `FloatDefault` components in its value and attempts to
+  /// `ArrayHandleBasic` with `vtkm::FloatDefault` components in its value and attempts to
   /// cast to those types.
   ///
   template <typename Functor, typename... Args>

vtkm/cont/UnknownArrayHandle.h

@@ -397,7 +397,7 @@ inline std::shared_ptr<UnknownAHContainer> MakeUnknownAHContainerFunctor::operat
 template <typename ValueTypeList, typename StorageTypeList>
 class UncertainArrayHandle;
 
-/// \brief An ArrayHandle of an unknown value type and storage.
+/// @brief An ArrayHandle of an unknown value type and storage.
 ///
 /// `UnknownArrayHandle` holds an `ArrayHandle` object using runtime polymorphism
 /// to manage different value and storage types rather than compile-time templates.
@@ -408,14 +408,14 @@ class UncertainArrayHandle;
 /// types.
 ///
 /// To interface between the runtime polymorphism and the templated algorithms
-/// in VTK-m, `UnknownArrayHandle` contains a method named `CastAndCallForTypes`
+/// in VTK-m, `UnknownArrayHandle` contains a method named `CastAndCallForTypes()`
 /// that determines the correct type from some known list of value types and
 /// storage. This mechanism is used internally by VTK-m's worklet invocation
 /// mechanism to determine the type when running algorithms.
 ///
 /// If the `UnknownArrayHandle` is used in a context where the possible array
 /// types can be whittled down to a finite list (or you have to), you can
-/// specify lists of value types and storage using the `ResetTypesAndStorage`
+/// specify lists of value types and storage using the `ResetTypesAndStorage()`
 /// method. This will convert this object to an `UncertainArrayHandle` of the
 /// given types. In cases where a finite set of types need to specified but
 /// there is no known subset, `VTKM_DEFAULT_TYPE_LIST` and
@@ -444,7 +444,7 @@ public:
   {
   }
 
-  /// \brief Returns whether an array is stored in this `UnknownArrayHandle`.
+  /// @brief Returns whether an array is stored in this `UnknownArrayHandle`.
   ///
   /// If the `UnknownArrayHandle` is constructed without an `ArrayHandle`, it
   /// will not have an underlying type, and therefore the operations will be
@@ -452,7 +452,7 @@ public:
   /// `ArrayHandle`.
   VTKM_CONT bool IsValid() const;
 
-  /// \brief Create a new array of the same type as this array.
+  /// @brief Create a new array of the same type as this array.
   ///
   /// This method creates a new array that is the same type as this one and
   /// returns a new `UnknownArrayHandle` for it. This method is convenient when
@@ -460,7 +460,7 @@ public:
   ///
   VTKM_CONT UnknownArrayHandle NewInstance() const;
 
-  /// \brief Create a new `ArrayHandleBasic` with the same `ValueType` as this array.
+  /// @brief Create a new `ArrayHandleBasic` with the same `ValueType` as this array.
   ///
   /// This method creates a new `ArrayHandleBasic` that has the same `ValueType` as the
   /// array held by this one and returns a new `UnknownArrayHandle` for it. This method
@@ -469,7 +469,7 @@ public:
   ///
   VTKM_CONT UnknownArrayHandle NewInstanceBasic() const;
 
-  /// \brief Create a new `ArrayHandleBasic` with the base component of `FloatDefault`
+  /// @brief Create a new `ArrayHandleBasic` with the base component of `vtkm::FloatDefault`
   ///
   /// This method creates a new `ArrayHandleBasic` that has a `ValueType` that is similar
   /// to the array held by this one except that the base component type is replaced with
@@ -487,22 +487,22 @@ public:
   ///
   VTKM_CONT UnknownArrayHandle NewInstanceFloatBasic() const;
 
-  /// \brief Returns the name of the value type stored in the array.
+  /// @brief Returns the name of the value type stored in the array.
   ///
   /// Returns an empty string if no array is stored.
   VTKM_CONT std::string GetValueTypeName() const;
 
-  /// \brief Returns the name of the base component of the value type stored in the array.
+  /// @brief Returns the name of the base component of the value type stored in the array.
   ///
   /// Returns an empty string if no array is stored.
   VTKM_CONT std::string GetBaseComponentTypeName() const;
 
-  /// \brief Returns the name of the storage tag for the array.
+  /// @brief Returns the name of the storage tag for the array.
   ///
   /// Returns an empty string if no array is stored.
   VTKM_CONT std::string GetStorageTypeName() const;
 
-  /// \brief Returns a string representation of the underlying data type.
+  /// @brief Returns a string representation of the underlying data type.
   ///
   /// The returned string will be of the form `vtkm::cont::ArrayHandle<T, S>` rather than the name
   /// of an actual subclass. If no array is stored, an empty string is returned.
@@ -525,7 +525,7 @@ public:
     return this->IsStorageTypeImpl(typeid(StorageType));
   }
 
-  /// \brief Returns true if this array's `ValueType` has the provided base component type.
+  /// @brief Returns true if this array's `ValueType` has the provided base component type.
   ///
   /// The base component type is the recursive component type of any `Vec`-like object. So
   /// if the array's `ValueType` is `vtkm::Vec<vtkm::Float32, 3>`, then the base component
@@ -542,17 +542,17 @@ public:
     return this->IsBaseComponentTypeImpl(detail::UnknownAHComponentInfo::Make<BaseComponentType>());
   }
 
-  /// \brief Returns true if this array matches the ArrayHandleType template argument.
+  /// @brief Returns true if this array matches the ArrayHandleType template argument.
   ///
   /// Note that `UnknownArrayHandle` has some special handling for `ArrayHandleCast` and
   /// `ArrayHandleMultiplexer`. If you stored an array of one of these types into an
-  /// `UnknownArrayHandle`, the type of the underlying array will change and `IsType`
+  /// `UnknownArrayHandle`, the type of the underlying array will change and `IsType()`
   /// will fail. However, you can still get the array back out as that type using
   /// `AsArrayHandle`.
   ///
-  /// Use the `CanConvert` method instead to determine if the `UnknownArrayHandle`
+  /// Use the `CanConvert()` method instead to determine if the `UnknownArrayHandle`
   /// contains an array that "matches" the array of a given type. Under most
-  /// circumstances, you should prefer `CanConvert` over `IsType`.
+  /// circumstances, you should prefer `CanConvert()` over `IsType()`.
   ///
   template <typename ArrayHandleType>
   VTKM_CONT bool IsType() const
@@ -562,7 +562,7 @@ public:
             this->IsStorageType<typename ArrayHandleType::StorageTag>());
   }
 
-  /// \brief Assigns potential value and storage types.
+  /// @brief Assigns potential value and storage types.
   ///
   /// Calling this method will return an `UncertainArrayHandle` with the provided
   /// value and storage type lists. The returned object will hold the same
@@ -575,11 +575,11 @@ public:
     NewValueTypeList = NewValueTypeList{},
     NewStorageTypeList = NewStorageTypeList{}) const;
 
-  /// \brief Returns the number of values in the array.
+  /// @brief Returns the number of values in the array.
   ///
   VTKM_CONT vtkm::Id GetNumberOfValues() const;
 
-  /// \brief Returns the number of components for each value in the array.
+  /// @brief Returns the number of components for each value in the array.
   ///
   /// If the array holds `vtkm::Vec` objects, this will return the number of components
   /// in each value. If the array holds a basic C type (such as `float`), this will return 1.
@@ -588,7 +588,7 @@ public:
   ///
   VTKM_CONT vtkm::IdComponent GetNumberOfComponents() const;
 
-  /// \brief Returns the total number of components for each value in the array.
+  /// @brief Returns the total number of components for each value in the array.
   ///
   /// If the array holds `vtkm::Vec` objects, this will return the total number of components
   /// in each value assuming the object is flattened out to one level of `Vec` objects.
@@ -606,21 +606,19 @@ public:
   ///
   VTKM_CONT vtkm::IdComponent GetNumberOfComponentsFlat() const;
 
-  /// \brief Reallocate the data in the array.
+  /// @brief Reallocate the data in the array.
   ///
-  /// The allocation works the same as the `Allocate` method of `vtkm::cont::ArrayHandle`.
-  ///
-  /// @{
+  /// The allocation works the same as the `Allocate()` method of `vtkm::cont::ArrayHandle`.
   VTKM_CONT void Allocate(vtkm::Id numValues,
                           vtkm::CopyFlag preserve,
                           vtkm::cont::Token& token) const;
+  /// @copydoc Allocate
   VTKM_CONT void Allocate(vtkm::Id numValues, vtkm::CopyFlag preserve = vtkm::CopyFlag::Off) const;
-  /// @}
 
-  /// \brief Determine if the contained array can be passed to the given array type.
+  /// @brief Determine if the contained array can be passed to the given array type.
   ///
-  /// This method will return true if calling `AsArrayHandle` of the given type will
-  /// succeed. The result is similar to `IsType`, and if `IsType` returns true, then
+  /// This method will return true if calling `AsArrayHandle()` of the given type will
+  /// succeed. The result is similar to `IsType()`, and if `IsType()` returns true, then
   /// this will return true. However, this method will also return true for other
   /// types such as an `ArrayHandleMultiplexer` that can contain the array.
   ///
@@ -651,24 +649,23 @@ private:
   }
 
 public:
-  ///@{
   /// Returns this array cast appropriately and stored in the given `ArrayHandle` type.
-  /// Throws an `ErrorBadType` if the stored array cannot be stored in the given array type.
-  /// Use the `CanConvert` method to determine if the array can be returned with the given type.
-  ///
+  /// Throws a `vtkm::cont::ErrorBadType` if the stored array cannot be stored in the given
+  /// array type. Use the `CanConvert()` method to determine if the array can be returned
+  /// with the given type.
   template <typename T, typename S>
   VTKM_CONT void AsArrayHandle(vtkm::cont::ArrayHandle<T, S>& array) const
   {
     this->BaseAsArrayHandle(array);
   }
-
+  /// @copydoc AsArrayHandle
   template <typename T>
   VTKM_CONT void AsArrayHandle(vtkm::cont::ArrayHandle<T>& array) const;
-
+  /// @copydoc AsArrayHandle
   template <typename T, typename... Ss>
   VTKM_CONT void AsArrayHandle(
     vtkm::cont::ArrayHandle<T, vtkm::cont::StorageTagMultiplexer<Ss...>>& array) const;
-
+  /// @copydoc AsArrayHandle
   template <typename TargetT, typename SourceT, typename SourceS>
   VTKM_CONT void AsArrayHandle(
     vtkm::cont::ArrayHandle<TargetT, vtkm::cont::StorageTagCast<SourceT, SourceS>>& array) const
@@ -677,7 +674,7 @@ public:
     array = vtkm::cont::ArrayHandleCast<TargetT, ContainedArrayType>(
       this->AsArrayHandle<ContainedArrayType>());
   }
-
+  /// @copydoc AsArrayHandle
   template <typename T>
   VTKM_CONT void AsArrayHandle(
     vtkm::cont::ArrayHandle<T, vtkm::cont::StorageTagRuntimeVec>& array) const
@@ -697,7 +694,7 @@ public:
       this->BaseAsArrayHandle(array);
     }
   }
-
+  /// @copydoc AsArrayHandle
   template <typename ArrayType>
   VTKM_CONT ArrayType AsArrayHandle() const
   {
@@ -706,12 +703,12 @@ public:
     this->AsArrayHandle(array);
     return array;
   }
-  ///@}
+
 #ifdef VTKM_MSVC
   VTKM_DEPRECATED_SUPPRESS_END
 #endif
 
-  /// \brief Deep copies data from another `UnknownArrayHandle`.
+  /// @brief Deep copies data from another `UnknownArrayHandle`.
   ///
   /// This method takes an `UnknownArrayHandle` and deep copies data from it.
   ///
@@ -720,60 +717,60 @@ public:
   ///
   void DeepCopyFrom(const vtkm::cont::UnknownArrayHandle& source);
 
-  /// \brief Deep copies data from another `UnknownArrayHandle`.
+  /// @brief Deep copies data from another `UnknownArrayHandle`.
   ///
   /// This method takes an `UnknownArrayHandle` and deep copies data from it.
   ///
   /// If this object does not point to an existing `ArrayHandle`, this const version
-  /// of `DeepCopyFrom` throws an exception.
+  /// of `DeepCopyFrom()` throws an exception.
   ///
   void DeepCopyFrom(const vtkm::cont::UnknownArrayHandle& source) const;
 
-  /// \brief Attempts a shallow copy of an array or a deep copy if that is not possible.
+  /// @brief Attempts a shallow copy of an array or a deep copy if that is not possible.
   ///
   /// This method takes an `UnknownArrayHandle` and attempts to perform a shallow copy.
   /// This shallow copy occurs if this object points to an `ArrayHandle` of the same type
   /// or does not point to any `ArrayHandle` at all. If this is not possible, then
   /// the array is deep copied.
   ///
-  /// This method is roughly equivalent to the `ArrayCopyShallowIfPossible` function
+  /// This method is roughly equivalent to the `vtkm::cont::ArrayCopyShallowIfPossible()` function
   /// (defined in `vtkm/cont/ArrayCopy.h`). However, this method can be used without
-  /// having to use a device compiler (whereas `ArrayCopyShallowIfPossible` does require
+  /// having to use a device compiler (whereas `vtkm::cont::ArrayCopyShallowIfPossible()` does require
   /// a device device compiler).
   ///
   void CopyShallowIfPossible(const vtkm::cont::UnknownArrayHandle& source);
 
-  /// \brief Attempts a shallow copy of an array or a deep copy if that is not possible.
+  /// @brief Attempts a shallow copy of an array or a deep copy if that is not possible.
   ///
   /// This method takes an `UnknownArrayHandle` and attempts to perform a shallow copy.
   /// This shallow copy occurs if this object points to an `ArrayHandle` of the same type.
   /// If the types are incompatible, then the array is deep copied.
   ///
   /// If this object does not point to an existing `ArrayHandle`, this const version
-  /// of `CopyShallowIfPossible` throws an exception.
+  /// of `CopyShallowIfPossible()` throws an exception.
   ///
-  /// This method is roughly equivalent to the `ArrayCopyShallowIfPossible` function
+  /// This method is roughly equivalent to the `vtkm::cont::ArrayCopyShallowIfPossible()` function
   /// (defined in `vtkm/cont/ArrayCopy.h`). However, this method can be used without
-  /// having to use a device compiler (whereas `ArrayCopyShallowIfPossible` does require
+  /// having to use a device compiler (whereas `vtkm::cont::ArrayCopyShallowIfPossible()` does require
   /// a device device compiler).
   ///
   void CopyShallowIfPossible(const vtkm::cont::UnknownArrayHandle& source) const;
 
-  /// \brief Extract a component of the array.
+  /// @brief Extract a component of the array.
   ///
   /// This method returns an array that holds the data for a given flat component of the data.
   /// The `BaseComponentType` has to be specified and must match the contained array (i.e.
-  /// the result of `IsBaseComponentType` must succeed for the given type).
+  /// the result of `IsBaseComponentType()` must succeed for the given type).
   ///
-  /// This method treats each value in the array as a flat `Vec` even if it is a `Vec` of
-  /// `Vec`s. For example, if the array actually holds values of type `Vec<Vec<T, 3>, 2>`,
-  /// it is treated as if it holds a `Vec<T, 6>`. See `vtkm::VecFlat` for details on how
-  /// vectors are flattened.
+  /// This method treats each value in the array as a flat `vtkm::Vec` even if it is a
+  /// `vtkm::Vec` of `Vec`s. For example, if the array actually holds values of type
+  /// `vtkm::Vec<vtkm::Vec<T, 3>, 2>`, it is treated as if it holds a `Vec<T, 6>`. See
+  /// `vtkm::VecFlat` for details on how vectors are flattened.
   ///
-  /// The point of using `ExtractComponent` over `AsArrayHandle` is that it drastically reduces
-  /// the amount of types you have to try. Most of the type the base component type is one of
+  /// The point of using `ExtractComponent()` over `AsArrayHandle()` is that it drastically reduces
+  /// the amount of types you have to try. Most of the time the base component type is one of
   /// the basic C types (i.e. `int`, `long`, `float`, etc.). You do not need to know what shape
-  /// the containing `Vec` is in, nor do you need to know the actual storage of the array.
+  /// the containing `vtkm::Vec` is in, nor do you need to know the actual storage of the array.
   ///
   /// Note that the type of the array returned is `ArrayHandleStride`. Using this type of
   /// array handle has a slight overhead over basic arrays like `ArrayHandleBasic` and
@@ -804,38 +801,39 @@ public:
     return ComponentArrayType(buffers);
   }
 
-  /// \brief Extract the array knowing only the component type of the array.
+  /// @brief Extract the array knowing only the component type of the array.
   ///
   /// This method returns an `ArrayHandle` that points to the data in the array. This method
-  /// differs from `AsArrayHandle` because you do not need to know the exact `ValueType` and
+  /// differs from `AsArrayHandle()` because you do not need to know the exact `ValueType` and
   /// `StorageTag` of the array. Instead, you only need to know the base component type.
   ///
-  /// `ExtractArrayFromComponents` works by calling the `ExtractComponent` method and then
+  /// `ExtractArrayFromComponents()` works by calling the `ExtractComponent()` method and then
   /// combining them together in a fancy `ArrayHandle`. This allows you to ignore the storage
   /// type of the underlying array as well as any `Vec` structure of the value type. However,
   /// it also places some limitations on how the data can be pulled from the data.
   ///
   /// First, you have to specify the base component type. This must match the data in the
-  /// underlying array (as reported by `IsBaseComponentType`).
+  /// underlying array (as reported by `IsBaseComponentType()`).
   ///
   /// Second, the array returned will have the `Vec`s flattened. For example, if the underlying
-  /// array has a `ValueType` of `Vec<Vec<T, 3>, 3>`, then this method will tread the data as
-  /// if it was `Vec<T, 9>`. There is no way to get an array with `Vec` of `Vec` values.
+  /// array has a `ValueType` of `vtkm::Vec<vtkm::Vec<T, 3>, 3>`, then this method will treat
+  /// the data as if it was `vtkm::Vec<T, 9>`. There is no way to get an array with `vtkm::Vec`
+  /// of `vtkm::Vec` values.
   ///
   /// Third, because the `Vec` length of the values in the returned `ArrayHandle` must be
-  /// determined at runtime, that can break many assumptions of using `Vec` objects. The
-  /// type is not going to be a `Vec<T,N>` type but rather an internal class that is intended
-  /// to behave like that. The type should behave mostly like a `Vec`, but will have some
+  /// determined at runtime, that can break many assumptions of using `vtkm::Vec` objects. The
+  /// type is not going to be a `vtkm::Vec<T,N>` type but rather an internal class that is intended
+  /// to behave like that. The type should behave mostly like a `vtkm::Vec`, but will have some
   /// differences that can lead to unexpected behavior. For example, this `Vec`-like object
   /// will not have a `NUM_COMPONENTS` constant static expression because it is not known
-  /// at compile time. (Use the `GetNumberOfComponents` method instead.) And for the same
+  /// at compile time. (Use the `GetNumberOfComponents()` method instead.) And for the same
   /// reason you will not be able to pass these objects to classes overloaded or templated
   /// on the `Vec` type. Also, these `Vec`-like objects cannot be created as new instances.
   /// Thus, you will likely have to iterate over all components rather than do operations on
   /// the whole `Vec`.
   ///
-  /// Fourth, because `ExtractArrayFromComponents` uses `ExtractComponent` to pull data from
-  /// the array (which in turn uses `ArrayExtractComponent`), there are some `ArrayHandle` types
+  /// Fourth, because `ExtractArrayFromComponents()` uses `ExtractComponent()` to pull data from
+  /// the array (which in turn uses `ArrayExtractComponent()`), there are some `ArrayHandle` types
   /// that will require copying data to a new array. This could be problematic in cases where
   /// you want to write to the array. To prevent data from being copied, set the optional
   ///  `allowCopy` to `vtkm::CopyFlag::Off`. This will cause an exception to be thrown if
@@ -858,9 +856,9 @@ public:
     return result;
   }
 
-  /// \brief Call a functor using the underlying array type.
+  /// @brief Call a functor using the underlying array type.
   ///
-  /// `CastAndCallForTypes` attempts to cast the held array to a specific value type,
+  /// `CastAndCallForTypes()` attempts to cast the held array to a specific value type,
   /// and then calls the given functor with the cast array. You must specify
   /// the `TypeList` and `StorageList` as template arguments.
   ///
@@ -870,9 +868,9 @@ public:
   template <typename TypeList, typename StorageList, typename Functor, typename... Args>
   VTKM_CONT void CastAndCallForTypes(Functor&& functor, Args&&... args) const;
 
-  /// \brief Call a functor using the underlying array type with a float cast fallback.
+  /// @brief Call a functor using the underlying array type with a float cast fallback.
   ///
-  /// `CastAndCallForTypesWithFloatFallback` attempts to cast the held array to a specific
+  /// `CastAndCallForTypesWithFloatFallback()` attempts to cast the held array to a specific
   /// value type, and then calls the given functor with the cast array. You must specify
   /// the `TypeList` and `StorageList` as template arguments.
   ///
@@ -880,28 +878,28 @@ public:
   /// passed to the functor after the converted `ArrayHandle`.
   ///
   /// If the underlying array does not match any of the requested array types, the
-  /// array is copied to a new `ArrayHandleBasic` with `FloatDefault` components
+  /// array is copied to a new `ArrayHandleBasic` with `vtkm::FloatDefault` components
   /// in its value and attempts to cast to those types.
   ///
   template <typename TypeList, typename StorageList, typename Functor, typename... Args>
   VTKM_CONT void CastAndCallForTypesWithFloatFallback(Functor&& functor, Args&&... args) const;
 
-  /// \brief Call a functor on an array extracted from the components.
+  /// @brief Call a functor on an array extracted from the components.
   ///
-  /// `CastAndCallWithExtractedArray` behaves similarly to `CastAndCallForTypes`.
+  /// `CastAndCallWithExtractedArray()` behaves similarly to `CastAndCallForTypes()`.
   /// It converts the contained data to an `ArrayHandle` and calls a functor with
   /// that `ArrayHandle` (and any number of optionally specified arguments).
   ///
-  /// The advantage of `CastAndCallWithExtractedArray` is that you do not need to
+  /// The advantage of `CastAndCallWithExtractedArray()` is that you do not need to
   /// specify any `TypeList` or `StorageList`. Instead, it internally uses
-  /// `ExtractArrayFromComponents` to work with most `ArrayHandle` types with only
-  /// about 10 instances of the functor. In contrast, calling `CastAndCallForTypes`
+  /// `ExtractArrayFromComponents()` to work with most `ArrayHandle` types with only
+  /// about 10 instances of the functor. In contrast, calling `CastAndCallForTypes()`
   /// with, for example, `VTKM_DEFAULT_TYPE_LIST` and `VTKM_DEFAULT_STORAGE_LIST`
   /// results in many more instances of the functor but handling many fewer types
   /// of `ArrayHandle`.
   ///
   /// There are, however, costs to using this method. Details of these costs are
-  /// documented for the `ExtractArrayFromComponents` method, but briefly they
+  /// documented for the `ExtractArrayFromComponents()` method, but briefly they
   /// are that `Vec` types get flattened, the resulting array has a strange `Vec`-like
   /// value type that has many limitations on its use, there is an overhead for
   /// retrieving each value from the array, and there is a potential that data
@@ -919,6 +917,7 @@ public:
   ///
   VTKM_CONT void ReleaseResources() const;
 
+  /// Prints a summary of the array's type, size, and contents.
   VTKM_CONT void PrintSummary(std::ostream& out, bool full = false) const;
 };
 

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