Merge topic 'improve_flying_edge_perf'

251bd82b8 Significantly improve FlyingEdges performance across all devices
fa9373801 Rework FlyingEdges::Pass1 to handle NUMA and CUDA requirements.
769a10b47 FlyingEdge Normal and Point generation occurs in Pass4
93d87e06f Optimize StructuredPointGradient for non boundary points.

Acked-by: Kitware Robot <kwrobot@kitware.com>
Merge-request: !2080

vtkm/cont/ArrayHandleView.h

@@ -274,7 +274,7 @@ public:
     {
       throw vtkm::cont::ErrorBadValue(
         "The value array must be pre-allocated before it is used for the "
-        "output of ArrayHandlePermutation.");
+        "output of ArrayHandleView.");
     }
 
     return PortalExecution(

vtkm/exec/BoundaryState.h

@@ -176,6 +176,24 @@ struct BoundaryState
   }
   //@}
 
+  //@{
+  /// Takes a local neighborhood index (in the ranges of -neighborhood size to neighborhood size)
+  /// and returns the ijk of the equivalent point in the full data set. If the given value is out
+  /// of range, the returned value is undefined.
+  ///
+  VTKM_EXEC vtkm::Id3 NeighborIndexToFullIndex(const vtkm::IdComponent3& neighbor) const
+  {
+    return this->IJK + neighbor;
+  }
+
+  VTKM_EXEC vtkm::Id3 NeighborIndexToFullIndex(vtkm::IdComponent neighborI,
+                                               vtkm::IdComponent neighborJ,
+                                               vtkm::IdComponent neighborK) const
+  {
+    return this->NeighborIndexToFullIndex(vtkm::make_Vec(neighborI, neighborJ, neighborK));
+  }
+  //@}
+
   //@{
   /// Takes a local neighborhood index (in the ranges of -neighborhood size to
   /// neighborhood size), clamps it to the dataset bounds, and returns a new
@@ -221,6 +239,24 @@ struct BoundaryState
   }
   //@}
 
+  //@{
+  /// Takes a local neighborhood index (in the ranges of -neighborhood size to neighborhood size)
+  /// and returns the flat index of the equivalent point in the full data set. If the given value
+  /// is out of range, the result is undefined.
+  ///
+  VTKM_EXEC vtkm::Id NeighborIndexToFlatIndex(const vtkm::IdComponent3& neighbor) const
+  {
+    vtkm::Id3 full = this->IJK + neighbor;
+    return (full[2] * this->PointDimensions[1] + full[1]) * this->PointDimensions[0] + full[0];
+  }
+
+  VTKM_EXEC vtkm::Id NeighborIndexToFlatIndex(vtkm::IdComponent neighborI,
+                                              vtkm::IdComponent neighborJ,
+                                              vtkm::IdComponent neighborK) const
+  {
+    return this->NeighborIndexToFlatIndex(vtkm::make_Vec(neighborI, neighborJ, neighborK));
+  }
+  //@}
   vtkm::Id3 IJK;
   vtkm::Id3 PointDimensions;
 };

vtkm/exec/FieldNeighborhood.h

@@ -50,12 +50,24 @@ struct FieldNeighborhood
     return Portal.Get(this->Boundary->NeighborIndexToFlatIndexClamp(i, j, k));
   }
 
+  VTKM_EXEC
+  ValueType GetUnchecked(vtkm::IdComponent i, vtkm::IdComponent j, vtkm::IdComponent k) const
+  {
+    return Portal.Get(this->Boundary->NeighborIndexToFlatIndex(i, j, k));
+  }
+
   VTKM_EXEC
   ValueType Get(const vtkm::Id3& ijk) const
   {
     return Portal.Get(this->Boundary->NeighborIndexToFlatIndexClamp(ijk));
   }
 
+  VTKM_EXEC
+  ValueType GetUnchecked(const vtkm::Id3& ijk) const
+  {
+    return Portal.Get(this->Boundary->NeighborIndexToFlatIndex(ijk));
+  }
+
   vtkm::exec::BoundaryState const* const Boundary;
   FieldPortalType Portal;
 };
@@ -82,12 +94,24 @@ struct FieldNeighborhood<vtkm::internal::ArrayPortalUniformPointCoordinates>
     return Portal.Get(this->Boundary->NeighborIndexToFullIndexClamp(i, j, k));
   }
 
+  VTKM_EXEC
+  ValueType GetUnchecked(vtkm::IdComponent i, vtkm::IdComponent j, vtkm::IdComponent k) const
+  {
+    return Portal.Get(this->Boundary->NeighborIndexToFullIndex(i, j, k));
+  }
+
   VTKM_EXEC
   ValueType Get(const vtkm::IdComponent3& ijk) const
   {
     return Portal.Get(this->Boundary->NeighborIndexToFullIndexClamp(ijk));
   }
 
+  VTKM_EXEC
+  ValueType GetUnchecked(const vtkm::IdComponent3& ijk) const
+  {
+    return Portal.Get(this->Boundary->NeighborIndexToFullIndex(ijk));
+  }
+
   vtkm::exec::BoundaryState const* const Boundary;
   vtkm::internal::ArrayPortalUniformPointCoordinates Portal;
 };

vtkm/filter/testing/UnitTestContourFilterNormals.cxx

@@ -17,6 +17,8 @@
 
 #include <vtkm/filter/Contour.h>
 
+#include <vtkm/io/VTKDataSetWriter.h>
+
 namespace vtkm_ut_mc_normals
 {
 
@@ -70,6 +72,10 @@ void TestNormals(const vtkm::cont::DataSet& dataset, bool structured)
     { 0.770536f, -0.421248f, -0.478356f }, { -0.736036f, -0.445244f, -0.509910f },
     { 0.123446f, -0.887088f, -0.444788f }, { 0.133328f, -0.397444f, -0.907889f }
   };
+  //Calculated using FlyingEdges and Y axis iteration which causes
+  //the points to be in a different order
+  const vtkm::Id fe_y_alg_ordering[numVerts] = { 0, 1,  3,  5,  4, 6,  2,  7,
+                                                 9, 12, 10, 13, 8, 14, 11, 15 };
 
   //Calculated using normals of the output triangles
   const vtkm::Vec3f fast[numVerts] = {
@@ -83,6 +89,19 @@ void TestNormals(const vtkm::cont::DataSet& dataset, bool structured)
     { 0.2164f, -0.9401f, -0.2635f }, { -0.1589f, -0.1642f, -0.9735f }
   };
 
+  //When using the Y axis algorithm the cells are generated in a different
+  //order.
+  const vtkm::Vec3f fast_fe_y[numVerts] = {
+    { -0.243433f, -0.429741f, -0.869519f }, { 0.158904f, 0.164214f, -0.973542f },
+    { -0.895292f, -0.390217f, -0.214903f }, { -0.895057f, 0.134692f, -0.425125f },
+    { 0.829547f, -0.418793f, -0.36941f },   { 0.846705f, 0.425787f, -0.319054f },
+    { 0.253811f, -0.853394f, -0.4553f },    { -0.216381f, 0.940084f, -0.263478f },
+    { -0.848579f, -0.35602f, 0.391362f },   { -0.93948f, 0.252957f, 0.231065f },
+    { 0.831549f, -0.472663f, 0.291744f },   { 0.910494f, 0.0298277f, 0.412446f },
+    { -0.362862f, -0.815464f, 0.450944f },  { 0.107848f, 0.958544f, 0.263748f },
+    { 0.135131f, -0.437674f, 0.888921f },   { -0.286251f, 0.172078f, 0.942576f }
+  };
+
   vtkm::cont::ArrayHandle<vtkm::Vec3f> normals;
 
   vtkm::filter::Contour mc;
@@ -97,16 +116,23 @@ void TestNormals(const vtkm::cont::DataSet& dataset, bool structured)
   result.GetField("normals").GetData().CopyTo(normals);
   VTKM_TEST_ASSERT(normals.GetNumberOfValues() == numVerts,
                    "Wrong number of values in normals field");
-  auto normalsPortal = normals.ReadPortal();
+
-  for (vtkm::Id i = 0; i < numVerts; ++i)
+  //determine if we are using flying edge Y axis algorithm by checking the first normal value that differs
+  const bool using_fe_y_alg_ordering =
+    test_equal(normals.ReadPortal().Get(2), expected[fe_y_alg_ordering[2]], 0.001);
   {
-    VTKM_TEST_ASSERT(test_equal(normalsPortal.Get(i), expected[i], 0.001),
+    auto normalPotals = normals.ReadPortal();
-                     "Result (",
+    for (vtkm::Id i = 0; i < numVerts; ++i)
-                     normalsPortal.Get(i),
+    {
-                     ") does not match expected value (",
+      auto expected_v = !using_fe_y_alg_ordering ? expected[i] : expected[fe_y_alg_ordering[i]];
-                     expected[i],
+      VTKM_TEST_ASSERT(test_equal(normalPotals.Get(i), expected_v, 0.001),
-                     ") vert ",
+                       "Result (",
-                     i);
+                       normalPotals.Get(i),
+                       ") does not match expected value (",
+                       expected_v,
+                       ") vert ",
+                       i);
+    }
   }
 
   // Test the other normals generation method
@@ -114,6 +140,10 @@ void TestNormals(const vtkm::cont::DataSet& dataset, bool structured)
   {
     mc.SetComputeFastNormalsForStructured(true);
     expected = fast;
+    if (using_fe_y_alg_ordering)
+    {
+      expected = fast_fe_y;
+    }
   }
   else
   {
@@ -125,11 +155,20 @@ void TestNormals(const vtkm::cont::DataSet& dataset, bool structured)
   result.GetField("normals").GetData().CopyTo(normals);
   VTKM_TEST_ASSERT(normals.GetNumberOfValues() == numVerts,
                    "Wrong number of values in normals field");
-  normalsPortal = normals.ReadPortal();
+
-  for (vtkm::Id i = 0; i < numVerts; ++i)
   {
-    VTKM_TEST_ASSERT(test_equal(normalsPortal.Get(i), expected[i], 0.001),
+    auto normalPotals = normals.ReadPortal();
-                     "Result does not match expected values");
+    for (vtkm::Id i = 0; i < numVerts; ++i)
+    {
+      bool equal = test_equal(normalPotals.Get(i), expected[i], 0.001);
+      VTKM_TEST_ASSERT(equal,
+                       "Result (",
+                       normalPotals.Get(i),
+                       ") does not match expected value (",
+                       expected[i],
+                       ") vert ",
+                       i);
+    }
   }
 }
 

vtkm/worklet/contour/CMakeLists.txt

@@ -16,7 +16,10 @@ set(headers
   FlyingEdgesPass1.h
   FlyingEdgesPass2.h
   FlyingEdgesPass4.h
-  FlyingEdgesPass5.h
+  FlyingEdgesPass4Common.h
+  FlyingEdgesPass4X.h
+  FlyingEdgesPass4XWithNormals.h
+  FlyingEdgesPass4Y.h
   FlyingEdgesTables.h
   MarchingCellTables.h
   MarchingCells.h

vtkm/worklet/contour/FlyingEdges.h

@@ -16,7 +16,6 @@
 #include <vtkm/worklet/contour/FlyingEdgesPass1.h>
 #include <vtkm/worklet/contour/FlyingEdgesPass2.h>
 #include <vtkm/worklet/contour/FlyingEdgesPass4.h>
-#include <vtkm/worklet/contour/FlyingEdgesPass5.h>
 
 #include <vtkm/cont/ArrayHandleGroupVec.h>
 #include <vtkm/cont/Invoker.h>
@@ -30,32 +29,6 @@ namespace flying_edges
 
 namespace detail
 {
-inline vtkm::cont::CellSetStructured<3> make_metaDataMesh3D(SumXAxis, const vtkm::Id3& pdims)
-{
-  vtkm::cont::CellSetStructured<3> metaDataMesh;
-  metaDataMesh.SetPointDimensions(vtkm::Id3{ pdims[1], pdims[2], 1 });
-  return metaDataMesh;
-}
-inline vtkm::cont::CellSetStructured<2> make_metaDataMesh2D(SumXAxis, const vtkm::Id3& pdims)
-{
-  vtkm::cont::CellSetStructured<2> metaDataMesh;
-  metaDataMesh.SetPointDimensions(vtkm::Id2{ pdims[1], pdims[2] });
-  return metaDataMesh;
-}
-
-inline vtkm::cont::CellSetStructured<3> make_metaDataMesh3D(SumYAxis, const vtkm::Id3& pdims)
-{
-  vtkm::cont::CellSetStructured<3> metaDataMesh;
-  metaDataMesh.SetPointDimensions(vtkm::Id3{ pdims[0], pdims[2], 1 });
-  return metaDataMesh;
-}
-inline vtkm::cont::CellSetStructured<2> make_metaDataMesh2D(SumYAxis, const vtkm::Id3& pdims)
-{
-  vtkm::cont::CellSetStructured<2> metaDataMesh;
-  metaDataMesh.SetPointDimensions(vtkm::Id2{ pdims[0], pdims[2] });
-  return metaDataMesh;
-}
-
 template <typename T, typename S>
 vtkm::Id extend_by(vtkm::cont::ArrayHandle<T, S>& handle, vtkm::Id size)
 {
@@ -75,7 +48,6 @@ vtkm::Id extend_by(vtkm::cont::ArrayHandle<T, S>& handle, vtkm::Id size)
 }
 }
 
-
 //----------------------------------------------------------------------------
 template <typename ValueType,
           typename StorageTagField,
@@ -92,25 +64,21 @@ vtkm::cont::CellSetSingleType<> execute(
   vtkm::cont::ArrayHandle<vtkm::Vec<NormalType, 3>, StorageTagNormals>& normals,
   vtkm::worklet::contour::CommonState& sharedState)
 {
-  //Tasks:
-  //2. Refactor how we map fields.
-  //   We need the ability unload everything in SharedState after
-  //   we have mapped all fields
-  //3. Support switching AxisToSum by running this whole thing in a TryExecute
-  //   Passes 5 can ignore this
-
-  using AxisToSum = SumXAxis;
-
   vtkm::cont::Invoker invoke;
 
+  vtkm::Vec3f origin, spacing;
+  { //extract out the origin and spacing as these are needed for Pass4 to properly
+    //interpolate the new points
+    auto portal = coordinateSystem.ReadPortal();
+    origin = portal.GetOrigin();
+    spacing = portal.GetSpacing();
+  }
   auto pdims = cells.GetPointDimensions();
 
   vtkm::cont::ArrayHandle<vtkm::UInt8> edgeCases;
   edgeCases.Allocate(coordinateSystem.GetNumberOfValues());
 
-  vtkm::cont::CellSetStructured<3> metaDataMesh3D = detail::make_metaDataMesh3D(AxisToSum{}, pdims);
+  vtkm::cont::CellSetStructured<2> metaDataMesh2D;
-  vtkm::cont::CellSetStructured<2> metaDataMesh2D = detail::make_metaDataMesh2D(AxisToSum{}, pdims);
-
   vtkm::cont::ArrayHandle<vtkm::Id> metaDataLinearSums; //per point of metaDataMesh
   vtkm::cont::ArrayHandle<vtkm::Id> metaDataMin;        //per point of metaDataMesh
   vtkm::cont::ArrayHandle<vtkm::Id> metaDataMax;        //per point of metaDataMesh
@@ -138,19 +106,44 @@ vtkm::cont::CellSetSingleType<> execute(
     // figure out where intersections along the row begins and ends
     // (i.e., gather information for computational trimming).
     //
-    // We mark everything as below as it is faster than having the worklet to it
+    {
-    vtkm::cont::Algorithm::Fill(edgeCases, static_cast<vtkm::UInt8>(FlyingEdges3D::Below));
+      VTKM_LOG_SCOPE(vtkm::cont::LogLevel::Perf, "FlyingEdges Pass1");
-    ComputePass1<ValueType, AxisToSum> worklet1(isoval, pdims);
+
-    invoke(worklet1, metaDataMesh3D, metaDataSums, metaDataMin, metaDataMax, edgeCases, inputField);
+      // We have different logic for CUDA compared to Shared memory systems
+      // since this is the first touch of lots of the arrays, and will effect
+      // NUMA perf.
+      //
+      // Additionally CUDA does significantly better when you do an initial fill
+      // and write only non-below values
+      //
+      ComputePass1<ValueType> worklet1(isoval, pdims);
+      vtkm::cont::TryExecuteOnDevice(invoke.GetDevice(),
+                                     launchComputePass1{},
+                                     worklet1,
+                                     inputField,
+                                     edgeCases,
+                                     metaDataMesh2D,
+                                     metaDataSums,
+                                     metaDataMin,
+                                     metaDataMax);
+    }
 
     //----------------------------------------------------------------------------
     // PASS 2: Process a single row of voxels/cells. Count the number of other
     // axis intersections by topological reasoning from previous edge cases.
     // Determine the number of primitives (i.e., triangles) generated from this
     // row. Use computational trimming to reduce work.
-    ComputePass2<AxisToSum> worklet2(pdims);
+    {
-    invoke(
+      VTKM_LOG_SCOPE(vtkm::cont::LogLevel::Perf, "FlyingEdges Pass2");
-      worklet2, metaDataMesh2D, metaDataSums, metaDataMin, metaDataMax, metaDataNumTris, edgeCases);
+      ComputePass2 worklet2(pdims);
+      invoke(worklet2,
+             metaDataMesh2D,
+             metaDataSums,
+             metaDataMin,
+             metaDataMax,
+             metaDataNumTris,
+             edgeCases);
+    }
 
     //----------------------------------------------------------------------------
     // PASS 3: Compute the number of points and triangles that each edge
@@ -164,50 +157,42 @@ vtkm::cont::CellSetSingleType<> execute(
       detail::extend_by(sharedState.CellIdMap, sumTris);
 
 
-      auto newPointSize =
+      vtkm::Id newPointSize =
         vtkm::cont::Algorithm::ScanExclusive(metaDataLinearSums, metaDataLinearSums);
       detail::extend_by(sharedState.InterpolationEdgeIds, newPointSize);
       detail::extend_by(sharedState.InterpolationWeights, newPointSize);
 
       //----------------------------------------------------------------------------
       // PASS 4: Process voxel rows and generate topology, and interpolation state
-      ComputePass4<ValueType, AxisToSum> worklet4(
+      {
-        isoval, pdims, multiContourCellOffset, multiContourPointOffset);
+        VTKM_LOG_SCOPE(vtkm::cont::LogLevel::Perf, "FlyingEdges Pass4");
-      invoke(worklet4,
-             metaDataMesh2D,
-             metaDataSums,
-             metaDataMin,
-             metaDataMax,
-             metaDataNumTris,
-             edgeCases,
-             inputField,
-             triangle_topology,
-             sharedState.InterpolationEdgeIds,
-             sharedState.InterpolationWeights,
-             sharedState.CellIdMap);
-    }
 
-    //----------------------------------------------------------------------------
+        launchComputePass4 pass4(
-    // PASS 5: Convert the edge interpolation information to point and normals
+          pdims, origin, spacing, multiContourCellOffset, multiContourPointOffset);
-    vtkm::Vec3f origin, spacing;
-    { //extract out the origin and spacing as these are needed for Pass5 to properly
-      //interpolate the new points
-      auto portal = coordinateSystem.ReadPortal();
-      origin = portal.GetOrigin();
-      spacing = portal.GetSpacing();
-    }
-    if (sharedState.GenerateNormals)
-    {
-      normals.Allocate(sharedState.InterpolationEdgeIds.GetNumberOfValues());
-    }
 
-    ComputePass5<ValueType> worklet5(pdims, origin, spacing, sharedState.GenerateNormals);
+        detail::extend_by(points, newPointSize);
-    invoke(worklet5,
+        if (sharedState.GenerateNormals)
-           sharedState.InterpolationEdgeIds,
+        {
-           sharedState.InterpolationWeights,
+          detail::extend_by(normals, newPointSize);
-           points,
+        }
-           inputField,
+
-           normals);
+        vtkm::cont::TryExecuteOnDevice(invoke.GetDevice(),
+                                       pass4,
+                                       newPointSize,
+                                       isoval,
+                                       inputField,
+                                       edgeCases,
+                                       metaDataMesh2D,
+                                       metaDataSums,
+                                       metaDataMin,
+                                       metaDataMax,
+                                       metaDataNumTris,
+                                       sharedState,
+                                       triangle_topology,
+                                       points,
+                                       normals);
+      }
+    }
   }
 
   vtkm::cont::CellSetSingleType<> outputCells;

vtkm/worklet/contour/FlyingEdgesHelpers.h

@@ -55,15 +55,33 @@ struct SumYAxis
   static constexpr vtkm::Id zindex = 2;
 };
 
-VTKM_EXEC inline vtkm::Id compute_num_pts(SumXAxis, vtkm::Id nx, vtkm::Id vtkmNotUsed(ny))
+template <typename Device>
+struct select_AxisToSum
 {
-  return nx;
+  using type = SumXAxis;
-}
+};
-VTKM_EXEC inline vtkm::Id compute_num_pts(SumYAxis, vtkm::Id vtkmNotUsed(nx), vtkm::Id ny)
+
+template <>
+struct select_AxisToSum<vtkm::cont::DeviceAdapterTagCuda>
 {
-  return ny;
+  using type = SumYAxis;
+};
+
+inline vtkm::cont::CellSetStructured<2> make_metaDataMesh2D(SumXAxis, const vtkm::Id3& pdims)
+{
+  vtkm::cont::CellSetStructured<2> metaDataMesh;
+  metaDataMesh.SetPointDimensions(vtkm::Id2{ pdims[1], pdims[2] });
+  return metaDataMesh;
 }
 
+inline vtkm::cont::CellSetStructured<2> make_metaDataMesh2D(SumYAxis, const vtkm::Id3& pdims)
+{
+  vtkm::cont::CellSetStructured<2> metaDataMesh;
+  metaDataMesh.SetPointDimensions(vtkm::Id2{ pdims[0], pdims[2] });
+  return metaDataMesh;
+}
+
+
 VTKM_EXEC inline vtkm::Id3 compute_ijk(SumXAxis, const vtkm::Id3& executionSpaceIJK)
 {
   return vtkm::Id3{ 0, executionSpaceIJK[0], executionSpaceIJK[1] };

vtkm/worklet/contour/FlyingEdgesPass1.h

@@ -43,16 +43,37 @@ namespace flying_edges
 * that is formed by the current and next point.
 *
 */
-template <typename T, typename AxisToSum>
-struct ComputePass1 : public vtkm::worklet::WorkletPointNeighborhood
-{
 
-  vtkm::Id NumberOfPoints = 0;
+template <typename Device, typename WholeEdgeField>
+inline VTKM_EXEC void write_edge(Device,
+                                 vtkm::Id write_index,
+                                 WholeEdgeField& edges,
+                                 vtkm::UInt8 edgeCase)
+{
+  edges.Set(write_index, edgeCase);
+}
+
+template <typename WholeEdgeField>
+inline VTKM_EXEC void write_edge(vtkm::cont::DeviceAdapterTagCuda,
+                                 vtkm::Id write_index,
+                                 WholeEdgeField& edges,
+                                 vtkm::UInt8 edgeCase)
+{
+  if (edgeCase != FlyingEdges3D::Below)
+  {
+    edges.Set(write_index, edgeCase);
+  }
+}
+
+template <typename T>
+struct ComputePass1 : public vtkm::worklet::WorkletVisitPointsWithCells
+{
+  vtkm::Id3 PointDims;
   T IsoValue;
 
   ComputePass1() {}
   ComputePass1(T value, const vtkm::Id3& pdims)
-    : NumberOfPoints(compute_num_pts(AxisToSum{}, pdims[0], pdims[1]))
+    : PointDims(pdims)
     , IsoValue(value)
   {
   }
@@ -63,37 +84,40 @@ struct ComputePass1 : public vtkm::worklet::WorkletPointNeighborhood
                                 FieldOut axis_max,
                                 WholeArrayInOut edgeData,
                                 WholeArrayIn data);
-  using ExecutionSignature = void(Boundary, _2, _3, _4, _5, _6);
+  using ExecutionSignature = void(ThreadIndices, _2, _3, _4, _5, _6, Device);
   using InputDomain = _1;
 
-  template <typename WholeEdgeField, typename WholeDataField>
+  template <typename ThreadIndices,
-  VTKM_EXEC void operator()(const vtkm::exec::BoundaryState& boundary,
+            typename WholeEdgeField,
+            typename WholeDataField,
+            typename Device>
+  VTKM_EXEC void operator()(const ThreadIndices& threadIndices,
                             vtkm::Id3& axis_sum,
                             vtkm::Id& axis_min,
                             vtkm::Id& axis_max,
                             WholeEdgeField& edges,
-                            const WholeDataField& field) const
+                            const WholeDataField& field,
+                            Device device) const
   {
+    using AxisToSum = typename select_AxisToSum<Device>::type;
 
-    const vtkm::Id3 ijk = compute_ijk(AxisToSum{}, boundary.IJK);
+    const vtkm::Id3 ijk = compute_ijk(AxisToSum{}, threadIndices.GetInputIndex3D());
-    const vtkm::Id3 dims = compute_pdims(AxisToSum{}, boundary.PointDimensions, NumberOfPoints);
+    const vtkm::Id3 dims = this->PointDims;
     const vtkm::Id startPos = compute_start(AxisToSum{}, ijk, dims);
     const vtkm::Id offset = compute_inc(AxisToSum{}, dims);
 
     const T value = this->IsoValue;
-    axis_min = this->NumberOfPoints;
+    axis_min = this->PointDims[AxisToSum::xindex];
     axis_max = 0;
     T s1 = field.Get(startPos);
     T s0 = s1;
     axis_sum = { 0, 0, 0 };
-    for (vtkm::Id i = 0; i < NumberOfPoints - 1; ++i)
+    const vtkm::Id end = this->PointDims[AxisToSum::xindex] - 1;
+    for (vtkm::Id i = 0; i < end; ++i)
     {
       s0 = s1;
       s1 = field.Get(startPos + (offset * (i + 1)));
 
-      // We don't explicit write the Below case as that ruins performance.
-      // It is better to initially fill everything as Below and only
-      // write the exceptions
       vtkm::UInt8 edgeCase = FlyingEdges3D::Below;
       if (s0 >= value)
       {
@@ -103,16 +127,14 @@ struct ComputePass1 : public vtkm::worklet::WorkletPointNeighborhood
       {
         edgeCase |= FlyingEdges3D::RightAbove;
       }
-      if (edgeCase != FlyingEdges3D::Below)
+
-      {
+      write_edge(device, startPos + (offset * i), edges, edgeCase);
-        edges.Set(startPos + (offset * i), edgeCase);
-      }
 
       if (edgeCase == FlyingEdges3D::LeftAbove || edgeCase == FlyingEdges3D::RightAbove)
       {
         axis_sum[AxisToSum::xindex] += 1; // increment number of intersections along axis
         axis_max = i + 1;
-        if (axis_min == this->NumberOfPoints)
+        if (axis_min == (end + 1))
         {
           axis_min = i;
         }
@@ -120,6 +142,40 @@ struct ComputePass1 : public vtkm::worklet::WorkletPointNeighborhood
     }
   }
 };
+
+struct launchComputePass1
+{
+  template <typename DeviceAdapterTag, typename T, typename StorageTagField, typename... Args>
+  VTKM_CONT bool operator()(DeviceAdapterTag device,
+                            const ComputePass1<T>& worklet,
+                            const vtkm::cont::ArrayHandle<T, StorageTagField>& inputField,
+                            vtkm::cont::ArrayHandle<vtkm::UInt8> edgeCases,
+                            vtkm::cont::CellSetStructured<2>& metaDataMesh2D,
+                            Args&&... args) const
+  {
+    vtkm::cont::Invoker invoke(device);
+    metaDataMesh2D = make_metaDataMesh2D(SumXAxis{}, worklet.PointDims);
+
+    invoke(worklet, metaDataMesh2D, std::forward<Args>(args)..., edgeCases, inputField);
+    return true;
+  }
+
+  template <typename T, typename StorageTagField, typename... Args>
+  VTKM_CONT bool operator()(vtkm::cont::DeviceAdapterTagCuda device,
+                            const ComputePass1<T>& worklet,
+                            const vtkm::cont::ArrayHandle<T, StorageTagField>& inputField,
+                            vtkm::cont::ArrayHandle<vtkm::UInt8> edgeCases,
+                            vtkm::cont::CellSetStructured<2>& metaDataMesh2D,
+                            Args&&... args) const
+  {
+    vtkm::cont::Invoker invoke(device);
+    metaDataMesh2D = make_metaDataMesh2D(SumYAxis{}, worklet.PointDims);
+
+    vtkm::cont::Algorithm::Fill(edgeCases, static_cast<vtkm::UInt8>(FlyingEdges3D::Below));
+    invoke(worklet, metaDataMesh2D, std::forward<Args>(args)..., edgeCases, inputField);
+    return true;
+  }
+};
 }
 }
 }

vtkm/worklet/contour/FlyingEdgesPass2.h

@@ -23,7 +23,6 @@ namespace worklet
 namespace flying_edges
 {
 
-template <typename AxisToSum>
 struct ComputePass2 : public vtkm::worklet::WorkletVisitCellsWithPoints
 {
   vtkm::Id3 PointDims;
@@ -40,20 +39,24 @@ struct ComputePass2 : public vtkm::worklet::WorkletVisitCellsWithPoints
                                 FieldInPoint axis_maxs,
                                 FieldOutCell cell_tri_count,
                                 WholeArrayIn edgeData);
-  using ExecutionSignature = void(ThreadIndices, _2, _3, _4, _5, _6);
+  using ExecutionSignature = void(ThreadIndices, _2, _3, _4, _5, _6, Device);
   using InputDomain = _1;
 
   template <typename ThreadIndices,
             typename WholeSumField,
             typename FieldInPointId,
-            typename WholeEdgeField>
+            typename WholeEdgeField,
+            typename Device>
   VTKM_EXEC void operator()(const ThreadIndices& threadIndices,
                             const WholeSumField& axis_sums,
                             const FieldInPointId& axis_mins,
                             const FieldInPointId& axis_maxs,
                             vtkm::Int32& cell_tri_count,
-                            const WholeEdgeField& edges) const
+                            const WholeEdgeField& edges,
+                            Device) const
   {
+    using AxisToSum = typename select_AxisToSum<Device>::type;
+
     // Pass 2. Traverse all cells in the meta data plane. This allows us to
     // easily grab the four edge cases bounding this voxel-row
 
@@ -134,7 +137,8 @@ struct ComputePass2 : public vtkm::worklet::WorkletVisitCellsWithPoints
         sums[AxisToSum::zindex] += edgeUses[8];
 
         // handle boundary
-        this->CountBoundaryEdgeUses(onBoundary, edgeUses, sums, adj_row_sum, adj_col_sum);
+        this->CountBoundaryEdgeUses(
+          AxisToSum{}, onBoundary, edgeUses, sums, adj_row_sum, adj_col_sum);
       }
     }
 
@@ -157,7 +161,9 @@ struct ComputePass2 : public vtkm::worklet::WorkletVisitCellsWithPoints
   //
   // Only on these boundaries do we write to the metaData of our neighbor
   // as it is safe as those
-  VTKM_EXEC inline void CountBoundaryEdgeUses(vtkm::Vec<bool, 3> onBoundary,
+  template <typename AxisToSum>
+  VTKM_EXEC inline void CountBoundaryEdgeUses(AxisToSum,
+                                              vtkm::Vec<bool, 3> onBoundary,
                                               vtkm::UInt8 const* const edgeUses,
                                               vtkm::Id3& sums,
                                               vtkm::Id3& adj_row_sum,

vtkm/worklet/contour/FlyingEdgesPass4.h

@@ -13,9 +13,10 @@
 #ifndef vtk_m_worklet_contour_flyingedges_pass4_h
 #define vtk_m_worklet_contour_flyingedges_pass4_h
 
-
+#include <vtkm/worklet/contour/FlyingEdgesPass4Common.h>
-#include <vtkm/worklet/contour/FlyingEdgesHelpers.h>
+#include <vtkm/worklet/contour/FlyingEdgesPass4X.h>
-#include <vtkm/worklet/contour/FlyingEdgesTables.h>
+#include <vtkm/worklet/contour/FlyingEdgesPass4XWithNormals.h>
+#include <vtkm/worklet/contour/FlyingEdgesPass4Y.h>
 
 namespace vtkm
 {
@@ -24,428 +25,152 @@ namespace worklet
 namespace flying_edges
 {
 
-VTKM_EXEC inline vtkm::Id3 compute_incs3d(const vtkm::Id3& dims)
+struct launchComputePass4
 {
-  return vtkm::Id3{ 1, dims[0], (dims[0] * dims[1]) };
-}
-
-template <typename T, typename WholeField, typename EdgeField, typename WeightField>
-VTKM_EXEC inline void interpolate_weight(T value,
-                                         const vtkm::Id2& iEdge,
-                                         vtkm::Id writeIndex,
-                                         const WholeField& field,
-                                         EdgeField& interpolatedEdgeIds,
-                                         WeightField& weights)
-{
-  interpolatedEdgeIds.Set(writeIndex, iEdge);
-  T s0 = field.Get(iEdge[0]);
-  T s1 = field.Get(iEdge[1]);
-
-  auto t = (value - s0) / (s1 - s0);
-  weights.Set(writeIndex, static_cast<vtkm::FloatDefault>(t));
-}
-
-VTKM_EXEC inline bool case_includes_axes(vtkm::UInt8 const* const edgeUses)
-{
-  return (edgeUses[0] != 0 || edgeUses[4] != 0 || edgeUses[8] != 0);
-}
-
-template <typename WholeConnField, typename WholeCellIdField>
-VTKM_EXEC inline void generate_tris(vtkm::Id inputCellId,
-                                    vtkm::UInt8 edgeCase,
-                                    vtkm::UInt8 numTris,
-                                    vtkm::Id* edgeIds,
-                                    vtkm::Id& triId,
-                                    const WholeConnField& conn,
-                                    const WholeCellIdField& cellIds)
-{
-  auto* edges = data::GetTriEdgeCases(edgeCase);
-  vtkm::Id edgeIndex = 1;
-  vtkm::Id index = static_cast<vtkm::Id>(triId) * 3;
-  for (vtkm::UInt8 i = 0; i < numTris; ++i)
-  {
-    cellIds.Set(triId + i, inputCellId);
-
-    //We use edgeIndex, edgeIndex+2, edgeIndex+1 to keep
-    //the same winding for the triangles that marching celss
-    //produced. By keeping the winding the same we make sure
-    //that 'fast' normals are consistent with the marching
-    //cells version
-    conn.Set(index, edgeIds[edges[edgeIndex]]);
-    conn.Set(index + 1, edgeIds[edges[edgeIndex + 2]]);
-    conn.Set(index + 2, edgeIds[edges[edgeIndex + 1]]);
-    index += 3;
-    edgeIndex += 3;
-  }
-  triId += numTris;
-}
-
-// Helper function to set up the point ids on voxel edges.
-//----------------------------------------------------------------------------
-template <typename AxisToSum, typename FieldInPointId3>
-VTKM_EXEC inline void init_voxelIds(AxisToSum,
-                                    vtkm::Id writeOffset,
-                                    vtkm::UInt8 edgeCase,
-                                    const FieldInPointId3& axis_sums,
-                                    vtkm::Id* edgeIds)
-{
-  auto* edgeUses = data::GetEdgeUses(edgeCase);
-  edgeIds[0] = writeOffset + axis_sums[0][AxisToSum::xindex]; // x-edges
-  edgeIds[1] = writeOffset + axis_sums[1][AxisToSum::xindex];
-  edgeIds[2] = writeOffset + axis_sums[3][AxisToSum::xindex];
-  edgeIds[3] = writeOffset + axis_sums[2][AxisToSum::xindex];
-  edgeIds[4] = writeOffset + axis_sums[0][AxisToSum::yindex]; // y-edges
-  edgeIds[5] = edgeIds[4] + edgeUses[4];
-  edgeIds[6] = writeOffset + axis_sums[3][AxisToSum::yindex];
-  edgeIds[7] = edgeIds[6] + edgeUses[6];
-  edgeIds[8] = writeOffset + axis_sums[0][AxisToSum::zindex]; // z-edges
-  edgeIds[9] = edgeIds[8] + edgeUses[8];
-  edgeIds[10] = writeOffset + axis_sums[1][AxisToSum::zindex];
-  edgeIds[11] = edgeIds[10] + edgeUses[10];
-}
-
-// Helper function to advance the point ids along voxel rows.
-//----------------------------------------------------------------------------
-VTKM_EXEC inline void advance_voxelIds(vtkm::UInt8 const* const edgeUses, vtkm::Id* edgeIds)
-{
-  edgeIds[0] += edgeUses[0]; // x-edges
-  edgeIds[1] += edgeUses[1];
-  edgeIds[2] += edgeUses[2];
-  edgeIds[3] += edgeUses[3];
-  edgeIds[4] += edgeUses[4]; // y-edges
-  edgeIds[5] = edgeIds[4] + edgeUses[5];
-  edgeIds[6] += edgeUses[6];
-  edgeIds[7] = edgeIds[6] + edgeUses[7];
-  edgeIds[8] += edgeUses[8]; // z-edges
-  edgeIds[9] = edgeIds[8] + edgeUses[9];
-  edgeIds[10] += edgeUses[10];
-  edgeIds[11] = edgeIds[10] + edgeUses[11];
-}
-
-template <typename T, typename AxisToSum>
-struct ComputePass4 : public vtkm::worklet::WorkletVisitCellsWithPoints
-{
-
   vtkm::Id3 PointDims;
-  T IsoValue;
+  vtkm::Vec3f Origin;
+  vtkm::Vec3f Spacing;
+
   vtkm::Id CellWriteOffset;
   vtkm::Id PointWriteOffset;
 
-  ComputePass4() {}
+  launchComputePass4(const vtkm::Id3& pdims,
-  ComputePass4(T value,
+                     const vtkm::Vec3f& origin,
-               const vtkm::Id3& pdims,
+                     const vtkm::Vec3f& spacing,
-               vtkm::Id multiContourCellOffset,
+                     vtkm::Id multiContourCellOffset,
-               vtkm::Id multiContourPointOffset)
+                     vtkm::Id multiContourPointOffset)
     : PointDims(pdims)
-    , IsoValue(value)
+    , Origin(origin)
+    , Spacing(spacing)
     , CellWriteOffset(multiContourCellOffset)
     , PointWriteOffset(multiContourPointOffset)
   {
   }
 
-  using ControlSignature = void(CellSetIn,
+  template <typename DeviceAdapterTag,
-                                FieldInPoint axis_sums,
+            typename T,
-                                FieldInPoint axis_mins,
+            typename StorageTagField,
-                                FieldInPoint axis_maxs,
+            typename MeshSums,
-                                WholeArrayIn cell_tri_count,
+            typename PointType,
-                                WholeArrayIn edgeData,
+            typename NormalType>
-                                WholeArrayIn data,
+  VTKM_CONT bool operator()(DeviceAdapterTag device,
-                                WholeArrayOut connectivity,
+                            vtkm::Id vtkmNotUsed(newPointSize),
-                                WholeArrayOut edgeIds,
+                            T isoval,
-                                WholeArrayOut weights,
+                            const vtkm::cont::ArrayHandle<T, StorageTagField>& inputField,
-                                WholeArrayOut inputCellIds);
+                            vtkm::cont::ArrayHandle<vtkm::UInt8> edgeCases,
-  using ExecutionSignature =
+                            vtkm::cont::CellSetStructured<2>& metaDataMesh2D,
-    void(ThreadIndices, _2, _3, _4, _5, _6, _7, _8, _9, _10, _11, WorkIndex);
+                            const MeshSums& metaDataSums,
-
+                            const vtkm::cont::ArrayHandle<vtkm::Id>& metaDataMin,
-  template <typename ThreadIndices,
+                            const vtkm::cont::ArrayHandle<vtkm::Id>& metaDataMax,
-            typename FieldInPointId3,
+                            const vtkm::cont::ArrayHandle<vtkm::Int32>& metaDataNumTris,
-            typename FieldInPointId,
+                            vtkm::worklet::contour::CommonState& sharedState,
-            typename WholeTriField,
+                            vtkm::cont::ArrayHandle<vtkm::Id>& triangle_topology,
-            typename WholeEdgeField,
+                            PointType& points,
-            typename WholeDataField,
+                            NormalType& normals) const
-            typename WholeConnField,
-            typename WholeEdgeIdField,
-            typename WholeWeightField,
-            typename WholeCellIdField>
-  VTKM_EXEC void operator()(const ThreadIndices& threadIndices,
-                            const FieldInPointId3& axis_sums,
-                            const FieldInPointId& axis_mins,
-                            const FieldInPointId& axis_maxs,
-                            const WholeTriField& cellTriCount,
-                            const WholeEdgeField& edges,
-                            const WholeDataField& field,
-                            const WholeConnField& conn,
-                            const WholeEdgeIdField& interpolatedEdgeIds,
-                            const WholeWeightField& weights,
-                            const WholeCellIdField& inputCellIds,
-                            vtkm::Id oidx) const
   {
-    //This works as cellTriCount was computed with ScanExtended
+    vtkm::cont::Invoker invoke(device);
-    //and therefore has one more entry than the number of cells
+    if (sharedState.GenerateNormals)
-    vtkm::Id cell_tri_offset = cellTriCount.Get(oidx);
-    vtkm::Id next_tri_offset = cellTriCount.Get(oidx + 1);
-    if (cell_tri_offset == next_tri_offset)
-    { //we produce nothing
-      return;
-    }
-    cell_tri_offset += this->CellWriteOffset;
-
-    // find adjusted trim values.
-    vtkm::Id left = vtkm::Min(axis_mins[0], axis_mins[1]);
-    left = vtkm::Min(left, axis_mins[2]);
-    left = vtkm::Min(left, axis_mins[3]);
-
-    vtkm::Id right = vtkm::Max(axis_maxs[0], axis_maxs[1]);
-    right = vtkm::Max(right, axis_maxs[2]);
-    right = vtkm::Max(right, axis_maxs[3]);
-
-    vtkm::Id3 ijk = compute_ijk(AxisToSum{}, threadIndices.GetInputIndex3D());
-    const vtkm::Id3 pdims = this->PointDims;
-    const vtkm::Id4 startPos = compute_neighbor_starts(AxisToSum{}, ijk, pdims);
-    const vtkm::Id axis_inc = compute_inc(AxisToSum{}, pdims);
-
-    if (left == pdims[AxisToSum::xindex] && right == 0)
     {
-      //verify that we have nothing to generate and early terminate.
+      ComputePass4XWithNormals<T> worklet4(isoval,
-      bool mins_same = (axis_mins[0] == axis_mins[1] && axis_mins[0] == axis_mins[2] &&
+                                           this->PointDims,
-                        axis_mins[0] == axis_mins[3]);
+                                           this->Origin,
-      bool maxs_same = (axis_maxs[0] == axis_maxs[1] && axis_maxs[0] == axis_maxs[2] &&
+                                           this->Spacing,
-                        axis_maxs[0] == axis_maxs[3]);
+                                           this->CellWriteOffset,
-      if (mins_same && maxs_same)
+                                           this->PointWriteOffset);
-      {
+      invoke(worklet4,
-        return;
+             metaDataMesh2D,
-      }
+             metaDataSums,
-      else
+             metaDataMin,
-      {
+             metaDataMax,
-        left = 0;
+             metaDataNumTris,
-        right = pdims[AxisToSum::xindex] - 1;
+             edgeCases,
-      }
+             inputField,
+             triangle_topology,
+             sharedState.InterpolationEdgeIds,
+             sharedState.InterpolationWeights,
+             sharedState.CellIdMap,
+             points,
+             normals);
     }
-
+    else
-    // The trim edges may need adjustment if the contour travels between rows
-    // of edges (without intersecting these edges). This means checking
-    // whether the trim faces at (left,right) made up of the edges intersect
-    // the contour.
-    adjustTrimBounds(pdims[AxisToSum::xindex] - 1, edges, startPos, axis_inc, left, right);
-    if (left == right)
     {
-      return;
+      ComputePass4X<T> worklet4(isoval,
+                                this->PointDims,
+                                this->Origin,
+                                this->Spacing,
+                                this->CellWriteOffset,
+                                this->PointWriteOffset);
+      invoke(worklet4,
+             metaDataMesh2D,
+             metaDataSums,
+             metaDataMin,
+             metaDataMax,
+             metaDataNumTris,
+             edgeCases,
+             inputField,
+             triangle_topology,
+             sharedState.InterpolationEdgeIds,
+             sharedState.InterpolationWeights,
+             sharedState.CellIdMap,
+             points);
     }
 
-    const vtkm::UInt8 yLoc =
+    return true;
-      (ijk[AxisToSum::yindex] < 1
-         ? FlyingEdges3D::MinBoundary
-         : (ijk[AxisToSum::yindex] >= (pdims[AxisToSum::yindex] - 2) ? FlyingEdges3D::MaxBoundary
-                                                                     : FlyingEdges3D::Interior));
-    const vtkm::UInt8 zLoc =
-      (ijk[AxisToSum::zindex] < 1
-         ? FlyingEdges3D::MinBoundary
-         : (ijk[AxisToSum::zindex] >= (pdims[AxisToSum::zindex] - 2) ? FlyingEdges3D::MaxBoundary
-                                                                     : FlyingEdges3D::Interior));
-    const vtkm::UInt8 yzLoc = static_cast<vtkm::UInt8>((yLoc << 2) | (zLoc << 4));
-
-
-
-    const vtkm::Id3 increments = compute_incs3d(pdims);
-    vtkm::Id edgeIds[12];
-    auto edgeCase = getEdgeCase(edges, startPos, (axis_inc * left));
-    init_voxelIds(AxisToSum{}, this->PointWriteOffset, edgeCase, axis_sums, edgeIds);
-    for (vtkm::Id i = left; i < right; ++i) // run along the trimmed voxels
-    {
-      ijk[AxisToSum::xindex] = i;
-      edgeCase = getEdgeCase(edges, startPos, (axis_inc * i));
-      vtkm::UInt8 numTris = data::GetNumberOfPrimitives(edgeCase);
-      if (numTris > 0)
-      {
-        //compute what the current cellId is
-        vtkm::Id cellId = compute_start(AxisToSum{}, ijk, pdims - vtkm::Id3{ 1, 1, 1 });
-
-        // Start by generating triangles for this case
-        generate_tris(cellId, edgeCase, numTris, edgeIds, cell_tri_offset, conn, inputCellIds);
-
-        // Now generate edgeIds and weights along voxel axes if needed. Remember to take
-        // boundary into account.
-        vtkm::UInt8 loc = static_cast<vtkm::UInt8>(
-          yzLoc | (i < 1 ? FlyingEdges3D::MinBoundary
-                         : (i >= (pdims[AxisToSum::xindex] - 2) ? FlyingEdges3D::MaxBoundary
-                                                                : FlyingEdges3D::Interior)));
-        auto* edgeUses = data::GetEdgeUses(edgeCase);
-        if (loc != FlyingEdges3D::Interior || case_includes_axes(edgeUses))
-        {
-          this->GenerateWeights(loc,
-                                field,
-                                interpolatedEdgeIds,
-                                weights,
-                                startPos,
-                                increments,
-                                (axis_inc * i),
-                                edgeUses,
-                                edgeIds);
-        }
-        advance_voxelIds(edgeUses, edgeIds);
-      }
-    }
   }
 
-  //----------------------------------------------------------------------------
+  template <typename T,
-  template <typename WholeDataField, typename WholeIEdgeField, typename WholeWeightField>
+            typename StorageTagField,
-  VTKM_EXEC inline void GenerateWeights(vtkm::UInt8 loc,
+            typename MeshSums,
-                                        const WholeDataField& field,
+            typename PointType,
-                                        const WholeIEdgeField& interpolatedEdgeIds,
+            typename NormalType>
-                                        const WholeWeightField& weights,
+  VTKM_CONT bool operator()(vtkm::cont::DeviceAdapterTagCuda device,
-                                        const vtkm::Id4& startPos,
+                            vtkm::Id newPointSize,
-                                        const vtkm::Id3& incs,
+                            T isoval,
-                                        vtkm::Id offset,
+                            const vtkm::cont::ArrayHandle<T, StorageTagField>& inputField,
-                                        vtkm::UInt8 const* const edgeUses,
+                            vtkm::cont::ArrayHandle<vtkm::UInt8> edgeCases,
-                                        vtkm::Id* edgeIds) const
+                            vtkm::cont::CellSetStructured<2>& metaDataMesh2D,
+                            const MeshSums& metaDataSums,
+                            const vtkm::cont::ArrayHandle<vtkm::Id>& metaDataMin,
+                            const vtkm::cont::ArrayHandle<vtkm::Id>& metaDataMax,
+                            const vtkm::cont::ArrayHandle<vtkm::Int32>& metaDataNumTris,
+                            vtkm::worklet::contour::CommonState& sharedState,
+                            vtkm::cont::ArrayHandle<vtkm::Id>& triangle_topology,
+                            PointType& points,
+                            NormalType& normals) const
   {
-    vtkm::Id2 pos(startPos[0] + offset, 0);
+    vtkm::cont::Invoker invoke(device);
-    //EdgesUses 0,4,8 work for Y axis
-    if (edgeUses[0])
-    { // edgesUses[0] == i axes edge
-      pos[1] = startPos[0] + offset + incs[AxisToSum::xindex];
-      interpolate_weight(this->IsoValue, pos, edgeIds[0], field, interpolatedEdgeIds, weights);
-    }
-    if (edgeUses[4])
-    { // edgesUses[4] == j axes edge
-      pos[1] = startPos[1] + offset;
-      interpolate_weight(this->IsoValue, pos, edgeIds[4], field, interpolatedEdgeIds, weights);
-    }
-    if (edgeUses[8])
-    { // edgesUses[8] == k axes edge
-      pos[1] = startPos[2] + offset;
-      interpolate_weight(this->IsoValue, pos, edgeIds[8], field, interpolatedEdgeIds, weights);
-    }
-    // On the boundary cells special work has to be done to cover the partial
-    // cell axes. These are boundary situations where the voxel axes is not
-    // fully formed. These situations occur on the +x,+y,+z volume
-    // boundaries. (The other cases fall through the default: case which is
-    // expected.)
-    //
-    // Note that loc is one of 27 regions in the volume, with (0,1,2)
-    // indicating (interior, min, max) along coordinate axes.
-    switch (loc)
-    {
-      case 2:
-      case 6:
-      case 18:
-      case 22: //+x
-        this->InterpolateEdge(
-          pos[0], incs, 5, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
-        this->InterpolateEdge(
-          pos[0], incs, 9, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
-        break;
-      case 8:
-      case 9:
-      case 24:
-      case 25: //+y
-        this->InterpolateEdge(
-          pos[0], incs, 1, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
-        this->InterpolateEdge(
-          pos[0], incs, 10, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
-        break;
-      case 32:
-      case 33:
-      case 36:
-      case 37: //+z
-        this->InterpolateEdge(
-          pos[0], incs, 2, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
-        this->InterpolateEdge(
-          pos[0], incs, 6, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
-        break;
-      case 10:
-      case 26: //+x +y
-        this->InterpolateEdge(
-          pos[0], incs, 1, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
-        this->InterpolateEdge(
-          pos[0], incs, 5, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
-        this->InterpolateEdge(
-          pos[0], incs, 9, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
-        this->InterpolateEdge(
-          pos[0], incs, 10, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
-        this->InterpolateEdge(
-          pos[0], incs, 11, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
-        break;
-      case 34:
-      case 38: //+x +z
-        this->InterpolateEdge(
-          pos[0], incs, 2, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
-        this->InterpolateEdge(
-          pos[0], incs, 5, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
-        this->InterpolateEdge(
-          pos[0], incs, 9, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
-        this->InterpolateEdge(
-          pos[0], incs, 6, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
-        this->InterpolateEdge(
-          pos[0], incs, 7, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
-        break;
-      case 40:
-      case 41: //+y +z
-        this->InterpolateEdge(
-          pos[0], incs, 1, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
-        this->InterpolateEdge(
-          pos[0], incs, 2, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
-        this->InterpolateEdge(
-          pos[0], incs, 3, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
-        this->InterpolateEdge(
-          pos[0], incs, 6, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
-        this->InterpolateEdge(
-          pos[0], incs, 10, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
-        break;
-      case 42: //+x +y +z happens no more than once per volume
-        this->InterpolateEdge(
-          pos[0], incs, 1, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
-        this->InterpolateEdge(
-          pos[0], incs, 2, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
-        this->InterpolateEdge(
-          pos[0], incs, 3, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
-        this->InterpolateEdge(
-          pos[0], incs, 5, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
-        this->InterpolateEdge(
-          pos[0], incs, 9, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
-        this->InterpolateEdge(
-          pos[0], incs, 10, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
-        this->InterpolateEdge(
-          pos[0], incs, 11, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
-        this->InterpolateEdge(
-          pos[0], incs, 6, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
-        this->InterpolateEdge(
-          pos[0], incs, 7, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
-        break;
-      default: // interior, or -x,-y,-z boundaries
-        return;
-    }
-  }
 
-  // Indicate whether voxel axes need processing for this case.
+    ComputePass4Y<T> worklet4(
-  //----------------------------------------------------------------------------
+      isoval, this->PointDims, this->CellWriteOffset, this->PointWriteOffset);
-  template <typename WholeField, typename WholeIEdgeField, typename WholeWeightField>
+    invoke(worklet4,
-  VTKM_EXEC inline void InterpolateEdge(vtkm::Id currentIdx,
+           metaDataMesh2D,
-                                        const vtkm::Id3& incs,
+           metaDataSums,
-                                        vtkm::Id edgeNum,
+           metaDataMin,
-                                        vtkm::UInt8 const* const edgeUses,
+           metaDataMax,
-                                        vtkm::Id* edgeIds,
+           metaDataNumTris,
-                                        const WholeField& field,
+           edgeCases,
-                                        const WholeIEdgeField& interpolatedEdgeIds,
+           inputField,
-                                        const WholeWeightField& weights) const
+           triangle_topology,
-  {
+           sharedState.InterpolationEdgeIds,
-    // if this edge is not used then get out
+           sharedState.InterpolationWeights,
-    if (!edgeUses[edgeNum])
+           sharedState.CellIdMap);
-    {
-      return;
-    }
-    const vtkm::Id writeIndex = edgeIds[edgeNum];
-    vtkm::Id2 iEdge;
 
-    // build the edge information
+    //This needs to be done on array handle view ( start = this->PointWriteOffset, len = newPointSize)
-    vtkm::Vec<vtkm::UInt8, 2> verts = data::GetVertMap(edgeNum);
+    ComputePass5Y<T> worklet5(this->PointDims,
+                              this->Origin,
+                              this->Spacing,
+                              this->PointWriteOffset,
+                              sharedState.GenerateNormals);
+    invoke(worklet5,
+           vtkm::cont::make_ArrayHandleView(
+             sharedState.InterpolationEdgeIds, this->PointWriteOffset, newPointSize),
+           vtkm::cont::make_ArrayHandleView(
+             sharedState.InterpolationWeights, this->PointWriteOffset, newPointSize),
+           vtkm::cont::make_ArrayHandleView(points, this->PointWriteOffset, newPointSize),
+           inputField,
+           normals);
 
-    vtkm::Id3 offsets = data::GetVertOffsets(AxisToSum{}, verts[0]);
+    return true;
-    iEdge[0] = currentIdx + vtkm::Dot(offsets, incs);
-    offsets = data::GetVertOffsets(AxisToSum{}, verts[1]);
-    iEdge[1] = currentIdx + vtkm::Dot(offsets, incs);
-
-    interpolate_weight(this->IsoValue, iEdge, writeIndex, field, interpolatedEdgeIds, weights);
   }
 };
 }

vtkm/worklet/contour/FlyingEdgesPass4Common.h

@@ -0,0 +1,190 @@
+
+//============================================================================
+//  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.
+//============================================================================
+
+
+#ifndef vtk_m_worklet_contour_flyingedges_pass4_common_h
+#define vtk_m_worklet_contour_flyingedges_pass4_common_h
+
+#include <vtkm/worklet/contour/FlyingEdgesHelpers.h>
+#include <vtkm/worklet/contour/FlyingEdgesTables.h>
+
+namespace vtkm
+{
+namespace worklet
+{
+namespace flying_edges
+{
+
+VTKM_EXEC inline vtkm::Id3 compute_incs3d(const vtkm::Id3& dims)
+{
+  return vtkm::Id3{ 1, dims[0], (dims[0] * dims[1]) };
+}
+
+
+VTKM_EXEC inline bool case_includes_axes(vtkm::UInt8 const* const edgeUses)
+{
+  return (edgeUses[0] != 0 || edgeUses[4] != 0 || edgeUses[8] != 0);
+}
+
+template <typename WholeConnField, typename WholeCellIdField>
+VTKM_EXEC inline void generate_tris(vtkm::Id inputCellId,
+                                    vtkm::UInt8 edgeCase,
+                                    vtkm::UInt8 numTris,
+                                    vtkm::Id* edgeIds,
+                                    vtkm::Id& triId,
+                                    const WholeConnField& conn,
+                                    const WholeCellIdField& cellIds)
+{
+  auto* edges = data::GetTriEdgeCases(edgeCase);
+  vtkm::Id edgeIndex = 1;
+  vtkm::Id index = static_cast<vtkm::Id>(triId) * 3;
+  for (vtkm::UInt8 i = 0; i < numTris; ++i)
+  {
+    cellIds.Set(triId + i, inputCellId);
+
+    //This keeps the same winding for the triangles that marching cells
+    //produced. By keeping the winding the same we make sure
+    //that 'fast' normals are consistent with the marching
+    //cells version
+    conn.Set(index, edgeIds[edges[edgeIndex]]);
+    conn.Set(index + 1, edgeIds[edges[edgeIndex + 2]]);
+    conn.Set(index + 2, edgeIds[edges[edgeIndex + 1]]);
+    index += 3;
+    edgeIndex += 3;
+  }
+  triId += numTris;
+}
+
+
+// Helper function to set up the point ids on voxel edges.
+//----------------------------------------------------------------------------
+template <typename AxisToSum, typename FieldInPointId3>
+VTKM_EXEC inline void init_voxelIds(AxisToSum,
+                                    vtkm::Id writeOffset,
+                                    vtkm::UInt8 edgeCase,
+                                    const FieldInPointId3& axis_sums,
+                                    vtkm::Id* edgeIds)
+{
+  auto* edgeUses = data::GetEdgeUses(edgeCase);
+  edgeIds[0] = writeOffset + axis_sums[0][AxisToSum::xindex]; // x-edges
+  edgeIds[1] = writeOffset + axis_sums[1][AxisToSum::xindex];
+  edgeIds[2] = writeOffset + axis_sums[3][AxisToSum::xindex];
+  edgeIds[3] = writeOffset + axis_sums[2][AxisToSum::xindex];
+  edgeIds[4] = writeOffset + axis_sums[0][AxisToSum::yindex]; // y-edges
+  edgeIds[5] = edgeIds[4] + edgeUses[4];
+  edgeIds[6] = writeOffset + axis_sums[3][AxisToSum::yindex];
+  edgeIds[7] = edgeIds[6] + edgeUses[6];
+  edgeIds[8] = writeOffset + axis_sums[0][AxisToSum::zindex]; // z-edges
+  edgeIds[9] = edgeIds[8] + edgeUses[8];
+  edgeIds[10] = writeOffset + axis_sums[1][AxisToSum::zindex];
+  edgeIds[11] = edgeIds[10] + edgeUses[10];
+}
+
+// Helper function to advance the point ids along voxel rows.
+//----------------------------------------------------------------------------
+VTKM_EXEC inline void advance_voxelIds(vtkm::UInt8 const* const edgeUses, vtkm::Id* edgeIds)
+{
+  edgeIds[0] += edgeUses[0]; // x-edges
+  edgeIds[1] += edgeUses[1];
+  edgeIds[2] += edgeUses[2];
+  edgeIds[3] += edgeUses[3];
+  edgeIds[4] += edgeUses[4]; // y-edges
+  edgeIds[5] = edgeIds[4] + edgeUses[5];
+  edgeIds[6] += edgeUses[6];
+  edgeIds[7] = edgeIds[6] + edgeUses[7];
+  edgeIds[8] += edgeUses[8]; // z-edges
+  edgeIds[9] = edgeIds[8] + edgeUses[9];
+  edgeIds[10] += edgeUses[10];
+  edgeIds[11] = edgeIds[10] + edgeUses[11];
+}
+
+//----------------------------------------------------------------------------
+struct Pass4TrimState
+{
+  vtkm::Id left, right;
+  vtkm::Id3 ijk;
+  vtkm::Id4 startPos;
+  vtkm::Id axis_inc;
+  vtkm::UInt8 yzLoc;
+  bool valid = true;
+
+  template <typename AxisToSum,
+            typename ThreadIndices,
+            typename FieldInPointId,
+            typename WholeEdgeField>
+  VTKM_EXEC Pass4TrimState(AxisToSum,
+                           const vtkm::Id3& pdims,
+                           const ThreadIndices& threadIndices,
+                           const FieldInPointId& axis_mins,
+                           const FieldInPointId& axis_maxs,
+                           const WholeEdgeField& edges)
+  {
+    // find adjusted trim values.
+    left = vtkm::Min(axis_mins[0], axis_mins[1]);
+    left = vtkm::Min(left, axis_mins[2]);
+    left = vtkm::Min(left, axis_mins[3]);
+
+    right = vtkm::Max(axis_maxs[0], axis_maxs[1]);
+    right = vtkm::Max(right, axis_maxs[2]);
+    right = vtkm::Max(right, axis_maxs[3]);
+
+    ijk = compute_ijk(AxisToSum{}, threadIndices.GetInputIndex3D());
+
+    startPos = compute_neighbor_starts(AxisToSum{}, ijk, pdims);
+    axis_inc = compute_inc(AxisToSum{}, pdims);
+
+    if (left == pdims[AxisToSum::xindex] && right == 0)
+    {
+      //verify that we have nothing to generate and early terminate.
+      bool mins_same = (axis_mins[0] == axis_mins[1] && axis_mins[0] == axis_mins[2] &&
+                        axis_mins[0] == axis_mins[3]);
+      bool maxs_same = (axis_maxs[0] == axis_maxs[1] && axis_maxs[0] == axis_maxs[2] &&
+                        axis_maxs[0] == axis_maxs[3]);
+      if (mins_same && maxs_same)
+      {
+        valid = false;
+        return;
+      }
+      else
+      {
+        left = 0;
+        right = pdims[AxisToSum::xindex] - 1;
+      }
+    }
+
+    // The trim edges may need adjustment if the contour travels between rows
+    // of edges (without intersecting these edges). This means checking
+    // whether the trim faces at (left,right) made up of the edges intersect
+    // the contour.
+    adjustTrimBounds(pdims[AxisToSum::xindex] - 1, edges, startPos, axis_inc, left, right);
+    if (left == right)
+    {
+      valid = false;
+      return;
+    }
+
+    const vtkm::UInt8 yLoc =
+      (ijk[AxisToSum::yindex] < 1
+         ? FlyingEdges3D::MinBoundary
+         : (ijk[AxisToSum::yindex] >= (pdims[AxisToSum::yindex] - 2) ? FlyingEdges3D::MaxBoundary
+                                                                     : FlyingEdges3D::Interior));
+    const vtkm::UInt8 zLoc =
+      (ijk[AxisToSum::zindex] < 1
+         ? FlyingEdges3D::MinBoundary
+         : (ijk[AxisToSum::zindex] >= (pdims[AxisToSum::zindex] - 2) ? FlyingEdges3D::MaxBoundary
+                                                                     : FlyingEdges3D::Interior));
+    yzLoc = static_cast<vtkm::UInt8>((yLoc << 2) | (zLoc << 4));
+  }
+};
+}
+}
+}
+#endif

vtkm/worklet/contour/FlyingEdgesPass4X.h

@@ -0,0 +1,388 @@
+
+//============================================================================
+//  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.
+//============================================================================
+
+
+#ifndef vtk_m_worklet_contour_flyingedges_pass4x_h
+#define vtk_m_worklet_contour_flyingedges_pass4x_h
+
+
+#include <vtkm/worklet/contour/FlyingEdgesHelpers.h>
+#include <vtkm/worklet/contour/FlyingEdgesTables.h>
+
+namespace vtkm
+{
+namespace worklet
+{
+namespace flying_edges
+{
+
+template <typename T>
+struct ComputePass4X : public vtkm::worklet::WorkletVisitCellsWithPoints
+{
+
+  vtkm::Id3 PointDims;
+  vtkm::Vec3f Origin;
+  vtkm::Vec3f Spacing;
+
+  T IsoValue;
+
+  vtkm::Id CellWriteOffset;
+  vtkm::Id PointWriteOffset;
+
+  ComputePass4X() {}
+  ComputePass4X(T value,
+                const vtkm::Id3& pdims,
+                const vtkm::Vec3f& origin,
+                const vtkm::Vec3f& spacing,
+                vtkm::Id multiContourCellOffset,
+                vtkm::Id multiContourPointOffset)
+    : PointDims(pdims)
+    , Origin(origin)
+    , Spacing(spacing)
+    , IsoValue(value)
+    , CellWriteOffset(multiContourCellOffset)
+    , PointWriteOffset(multiContourPointOffset)
+  {
+  }
+
+  using ControlSignature = void(CellSetIn,
+                                FieldInPoint axis_sums,
+                                FieldInPoint axis_mins,
+                                FieldInPoint axis_maxs,
+                                WholeArrayIn cell_tri_count,
+                                WholeArrayIn edgeData,
+                                WholeArrayIn data,
+                                WholeArrayOut connectivity,
+                                WholeArrayOut edgeIds,
+                                WholeArrayOut weights,
+                                WholeArrayOut inputCellIds,
+                                WholeArrayOut points);
+  using ExecutionSignature =
+    void(ThreadIndices, _2, _3, _4, _5, _6, _7, _8, _9, _10, _11, _12, WorkIndex);
+
+  template <typename ThreadIndices,
+            typename FieldInPointId3,
+            typename FieldInPointId,
+            typename WholeTriField,
+            typename WholeEdgeField,
+            typename WholeDataField,
+            typename WholeConnField,
+            typename WholeEdgeIdField,
+            typename WholeWeightField,
+            typename WholeCellIdField,
+            typename WholePointField>
+  VTKM_EXEC void operator()(const ThreadIndices& threadIndices,
+                            const FieldInPointId3& axis_sums,
+                            const FieldInPointId& axis_mins,
+                            const FieldInPointId& axis_maxs,
+                            const WholeTriField& cellTriCount,
+                            const WholeEdgeField& edges,
+                            const WholeDataField& field,
+                            const WholeConnField& conn,
+                            const WholeEdgeIdField& interpolatedEdgeIds,
+                            const WholeWeightField& weights,
+                            const WholeCellIdField& inputCellIds,
+                            const WholePointField& points,
+                            vtkm::Id oidx) const
+  {
+    using AxisToSum = SumXAxis;
+
+    //This works as cellTriCount was computed with ScanExtended
+    //and therefore has one more entry than the number of cells
+    vtkm::Id cell_tri_offset = cellTriCount.Get(oidx);
+    vtkm::Id next_tri_offset = cellTriCount.Get(oidx + 1);
+    if (cell_tri_offset == next_tri_offset)
+    { //we produce nothing
+      return;
+    }
+    cell_tri_offset += this->CellWriteOffset;
+
+    const Pass4TrimState state(
+      AxisToSum{}, this->PointDims, threadIndices, axis_mins, axis_maxs, edges);
+    if (!state.valid)
+    {
+      return;
+    }
+
+    const vtkm::Id3 pdims = this->PointDims;
+    const vtkm::Id3 increments = compute_incs3d(pdims);
+    vtkm::Id edgeIds[12];
+
+    auto edgeCase = getEdgeCase(edges, state.startPos, (state.axis_inc * state.left));
+    init_voxelIds(AxisToSum{}, this->PointWriteOffset, edgeCase, axis_sums, edgeIds);
+    for (vtkm::Id i = state.left; i < state.right; ++i) // run along the trimmed voxels
+    {
+      auto ijk = state.ijk;
+      ijk[AxisToSum::xindex] = i;
+      edgeCase = getEdgeCase(edges, state.startPos, (state.axis_inc * i));
+      vtkm::UInt8 numTris = data::GetNumberOfPrimitives(edgeCase);
+      if (numTris > 0)
+      {
+        //compute what the current cellId is
+        vtkm::Id cellId = compute_start(AxisToSum{}, ijk, pdims - vtkm::Id3{ 1, 1, 1 });
+
+        // Start by generating triangles for this case
+        generate_tris(cellId, edgeCase, numTris, edgeIds, cell_tri_offset, conn, inputCellIds);
+
+        // Now generate edgeIds and weights along voxel axes if needed. Remember to take
+        // boundary into account.
+        vtkm::UInt8 loc = static_cast<vtkm::UInt8>(
+          state.yzLoc | (i < 1 ? FlyingEdges3D::MinBoundary
+                               : (i >= (pdims[AxisToSum::xindex] - 2) ? FlyingEdges3D::MaxBoundary
+                                                                      : FlyingEdges3D::Interior)));
+        auto* edgeUses = data::GetEdgeUses(edgeCase);
+        if (loc != FlyingEdges3D::Interior || case_includes_axes(edgeUses))
+        {
+          this->Generate(loc,
+                         ijk,
+                         field,
+                         interpolatedEdgeIds,
+                         weights,
+                         points,
+                         state.startPos,
+                         increments,
+                         (state.axis_inc * i),
+                         edgeUses,
+                         edgeIds);
+        }
+        advance_voxelIds(edgeUses, edgeIds);
+      }
+    }
+  }
+
+  //----------------------------------------------------------------------------
+  template <typename WholeDataField,
+            typename WholeIEdgeField,
+            typename WholeWeightField,
+            typename WholePointField>
+  VTKM_EXEC inline void Generate(vtkm::UInt8 loc,
+                                 const vtkm::Id3& ijk,
+                                 const WholeDataField& field,
+                                 const WholeIEdgeField& interpolatedEdgeIds,
+                                 const WholeWeightField& weights,
+                                 const WholePointField& points,
+                                 const vtkm::Id4& startPos,
+                                 const vtkm::Id3& incs,
+                                 vtkm::Id offset,
+                                 vtkm::UInt8 const* const edgeUses,
+                                 vtkm::Id* edgeIds) const
+  {
+    using AxisToSum = SumXAxis;
+
+    vtkm::Id2 pos(startPos[0] + offset, 0);
+    {
+      auto s0 = field.Get(pos[0]);
+
+      //EdgesUses 0,4,8 work for Y axis
+      if (edgeUses[0])
+      { // edgesUses[0] == i axes edge
+        auto writeIndex = edgeIds[0];
+        pos[1] = startPos[0] + offset + incs[AxisToSum::xindex];
+        auto s1 = field.Get(pos[1]);
+        T t = static_cast<T>((this->IsoValue - s0) / (s1 - s0));
+
+        interpolatedEdgeIds.Set(writeIndex, pos);
+        weights.Set(writeIndex, static_cast<vtkm::FloatDefault>(t));
+
+        auto coord = this->InterpolateCoordinate(t, ijk, ijk + vtkm::Id3{ 1, 0, 0 });
+        points.Set(writeIndex, coord);
+      }
+      if (edgeUses[4])
+      { // edgesUses[4] == j axes edge
+        auto writeIndex = edgeIds[4];
+        pos[1] = startPos[1] + offset;
+        auto s1 = field.Get(pos[1]);
+        T t = static_cast<T>((this->IsoValue - s0) / (s1 - s0));
+
+        interpolatedEdgeIds.Set(writeIndex, pos);
+        weights.Set(writeIndex, static_cast<vtkm::FloatDefault>(t));
+
+        auto coord = this->InterpolateCoordinate(t, ijk, ijk + vtkm::Id3{ 0, 1, 0 });
+        points.Set(writeIndex, coord);
+      }
+      if (edgeUses[8])
+      { // edgesUses[8] == k axes edge
+        auto writeIndex = edgeIds[8];
+        pos[1] = startPos[2] + offset;
+        auto s1 = field.Get(pos[1]);
+        T t = static_cast<T>((this->IsoValue - s0) / (s1 - s0));
+
+        interpolatedEdgeIds.Set(writeIndex, pos);
+        weights.Set(writeIndex, static_cast<vtkm::FloatDefault>(t));
+
+        auto coord = this->InterpolateCoordinate(t, ijk, ijk + vtkm::Id3{ 0, 0, 1 });
+        points.Set(writeIndex, coord);
+      }
+    }
+    // On the boundary cells special work has to be done to cover the partial
+    // cell axes. These are boundary situations where the voxel axes is not
+    // fully formed. These situations occur on the +x,+y,+z volume
+    // boundaries. (The other cases fall through the default: case which is
+    // expected.)
+    //
+    // Note that loc is one of 27 regions in the volume, with (0,1,2)
+    // indicating (interior, min, max) along coordinate axes.
+    switch (loc)
+    {
+      case 2:
+      case 6:
+      case 18:
+      case 22: //+x
+        this->InterpolateEdge(
+          ijk, pos[0], incs, 5, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
+        this->InterpolateEdge(
+          ijk, pos[0], incs, 9, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
+        break;
+      case 8:
+      case 9:
+      case 24:
+      case 25: //+y
+        this->InterpolateEdge(
+          ijk, pos[0], incs, 1, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
+        this->InterpolateEdge(
+          ijk, pos[0], incs, 10, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
+        break;
+      case 32:
+      case 33:
+      case 36:
+      case 37: //+z
+        this->InterpolateEdge(
+          ijk, pos[0], incs, 2, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
+        this->InterpolateEdge(
+          ijk, pos[0], incs, 6, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
+        break;
+      case 10:
+      case 26: //+x +y
+        this->InterpolateEdge(
+          ijk, pos[0], incs, 1, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
+        this->InterpolateEdge(
+          ijk, pos[0], incs, 5, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
+        this->InterpolateEdge(
+          ijk, pos[0], incs, 9, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
+        this->InterpolateEdge(
+          ijk, pos[0], incs, 10, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
+        this->InterpolateEdge(
+          ijk, pos[0], incs, 11, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
+        break;
+      case 34:
+      case 38: //+x +z
+        this->InterpolateEdge(
+          ijk, pos[0], incs, 2, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
+        this->InterpolateEdge(
+          ijk, pos[0], incs, 5, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
+        this->InterpolateEdge(
+          ijk, pos[0], incs, 9, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
+        this->InterpolateEdge(
+          ijk, pos[0], incs, 6, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
+        this->InterpolateEdge(
+          ijk, pos[0], incs, 7, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
+        break;
+      case 40:
+      case 41: //+y +z
+        this->InterpolateEdge(
+          ijk, pos[0], incs, 1, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
+        this->InterpolateEdge(
+          ijk, pos[0], incs, 2, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
+        this->InterpolateEdge(
+          ijk, pos[0], incs, 3, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
+        this->InterpolateEdge(
+          ijk, pos[0], incs, 6, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
+        this->InterpolateEdge(
+          ijk, pos[0], incs, 10, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
+        break;
+      case 42: //+x +y +z happens no more than once per volume
+        this->InterpolateEdge(
+          ijk, pos[0], incs, 1, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
+        this->InterpolateEdge(
+          ijk, pos[0], incs, 2, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
+        this->InterpolateEdge(
+          ijk, pos[0], incs, 3, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
+        this->InterpolateEdge(
+          ijk, pos[0], incs, 5, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
+        this->InterpolateEdge(
+          ijk, pos[0], incs, 9, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
+        this->InterpolateEdge(
+          ijk, pos[0], incs, 10, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
+        this->InterpolateEdge(
+          ijk, pos[0], incs, 11, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
+        this->InterpolateEdge(
+          ijk, pos[0], incs, 6, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
+        this->InterpolateEdge(
+          ijk, pos[0], incs, 7, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
+        break;
+      default: // interior, or -x,-y,-z boundaries
+        return;
+    }
+  }
+
+  // Indicate whether voxel axes need processing for this case.
+  //----------------------------------------------------------------------------
+  template <typename WholeField,
+            typename WholeIEdgeField,
+            typename WholeWeightField,
+            typename WholePointField>
+  VTKM_EXEC inline void InterpolateEdge(const vtkm::Id3& ijk,
+                                        vtkm::Id currentIdx,
+                                        const vtkm::Id3& incs,
+                                        vtkm::Id edgeNum,
+                                        vtkm::UInt8 const* const edgeUses,
+                                        vtkm::Id* edgeIds,
+                                        const WholeField& field,
+                                        const WholeIEdgeField& interpolatedEdgeIds,
+                                        const WholeWeightField& weights,
+                                        const WholePointField& points) const
+  {
+    using AxisToSum = SumXAxis;
+
+    // if this edge is not used then get out
+    if (!edgeUses[edgeNum])
+    {
+      return;
+    }
+    const vtkm::Id writeIndex = edgeIds[edgeNum];
+
+    // build the edge information
+    vtkm::Vec<vtkm::UInt8, 2> verts = data::GetVertMap(edgeNum);
+
+    vtkm::Id3 offsets1 = data::GetVertOffsets(AxisToSum{}, verts[0]);
+    vtkm::Id3 offsets2 = data::GetVertOffsets(AxisToSum{}, verts[1]);
+
+    vtkm::Id2 iEdge(currentIdx + vtkm::Dot(offsets1, incs), currentIdx + vtkm::Dot(offsets2, incs));
+
+    interpolatedEdgeIds.Set(writeIndex, iEdge);
+
+    auto s0 = field.Get(iEdge[0]);
+    auto s1 = field.Get(iEdge[1]);
+    T t = static_cast<T>((this->IsoValue - s0) / (s1 - s0));
+    weights.Set(writeIndex, static_cast<vtkm::FloatDefault>(t));
+
+    auto coord = this->InterpolateCoordinate(t, ijk + offsets1, ijk + offsets2);
+    points.Set(writeIndex, coord);
+  }
+
+  //----------------------------------------------------------------------------
+  inline VTKM_EXEC vtkm::Vec3f InterpolateCoordinate(T t,
+                                                     const vtkm::Id3& ijk0,
+                                                     const vtkm::Id3& ijk1) const
+  {
+    return vtkm::Vec3f(
+      this->Origin[0] +
+        this->Spacing[0] * static_cast<vtkm::FloatDefault>(ijk0[0] + t * (ijk1[0] - ijk0[0])),
+      this->Origin[1] +
+        this->Spacing[1] * static_cast<vtkm::FloatDefault>(ijk0[1] + t * (ijk1[1] - ijk0[1])),
+      this->Origin[2] +
+        this->Spacing[2] * static_cast<vtkm::FloatDefault>(ijk0[2] + t * (ijk1[2] - ijk0[2])));
+  }
+};
+}
+}
+}
+#endif

vtkm/worklet/contour/FlyingEdgesPass4XWithNormals.h

@@ -0,0 +1,468 @@
+
+//============================================================================
+//  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.
+//============================================================================
+
+
+#ifndef vtk_m_worklet_contour_flyingedges_pass4x_with_norms_h
+#define vtk_m_worklet_contour_flyingedges_pass4x_with_norms_h
+
+
+#include <vtkm/worklet/contour/FlyingEdgesHelpers.h>
+#include <vtkm/worklet/contour/FlyingEdgesTables.h>
+
+#include <vtkm/worklet/contour/FlyingEdgesPass4.h>
+
+namespace vtkm
+{
+namespace worklet
+{
+namespace flying_edges
+{
+
+template <typename T>
+struct ComputePass4XWithNormals : public vtkm::worklet::WorkletVisitCellsWithPoints
+{
+
+  vtkm::Id3 PointDims;
+  vtkm::Vec3f Origin;
+  vtkm::Vec3f Spacing;
+
+  T IsoValue;
+
+  vtkm::Id CellWriteOffset;
+  vtkm::Id PointWriteOffset;
+
+  ComputePass4XWithNormals() {}
+  ComputePass4XWithNormals(T value,
+                           const vtkm::Id3& pdims,
+                           const vtkm::Vec3f& origin,
+                           const vtkm::Vec3f& spacing,
+                           vtkm::Id multiContourCellOffset,
+                           vtkm::Id multiContourPointOffset)
+    : PointDims(pdims)
+    , Origin(origin)
+    , Spacing(spacing)
+    , IsoValue(value)
+    , CellWriteOffset(multiContourCellOffset)
+    , PointWriteOffset(multiContourPointOffset)
+  {
+  }
+
+  using ControlSignature = void(CellSetIn,
+                                FieldInPoint axis_sums,
+                                FieldInPoint axis_mins,
+                                FieldInPoint axis_maxs,
+                                WholeArrayIn cell_tri_count,
+                                WholeArrayIn edgeData,
+                                WholeArrayIn data,
+                                WholeArrayOut connectivity,
+                                WholeArrayOut edgeIds,
+                                WholeArrayOut weights,
+                                WholeArrayOut inputCellIds,
+                                WholeArrayOut points,
+                                WholeArrayOut normals);
+  using ExecutionSignature =
+    void(ThreadIndices, _2, _3, _4, _5, _6, _7, _8, _9, _10, _11, _12, _13, WorkIndex);
+
+  template <typename ThreadIndices,
+            typename FieldInPointId3,
+            typename FieldInPointId,
+            typename WholeTriField,
+            typename WholeEdgeField,
+            typename WholeDataField,
+            typename WholeConnField,
+            typename WholeEdgeIdField,
+            typename WholeWeightField,
+            typename WholeCellIdField,
+            typename WholePointField,
+            typename WholeNormalsField>
+  VTKM_EXEC void operator()(const ThreadIndices& threadIndices,
+                            const FieldInPointId3& axis_sums,
+                            const FieldInPointId& axis_mins,
+                            const FieldInPointId& axis_maxs,
+                            const WholeTriField& cellTriCount,
+                            const WholeEdgeField& edges,
+                            const WholeDataField& field,
+                            const WholeConnField& conn,
+                            const WholeEdgeIdField& interpolatedEdgeIds,
+                            const WholeWeightField& weights,
+                            const WholeCellIdField& inputCellIds,
+                            const WholePointField& points,
+                            const WholeNormalsField& normals,
+                            vtkm::Id oidx) const
+  {
+    using AxisToSum = SumXAxis;
+
+    //This works as cellTriCount was computed with ScanExtended
+    //and therefore has one more entry than the number of cells
+    vtkm::Id cell_tri_offset = cellTriCount.Get(oidx);
+    vtkm::Id next_tri_offset = cellTriCount.Get(oidx + 1);
+    if (cell_tri_offset == next_tri_offset)
+    { //we produce nothing
+      return;
+    }
+    cell_tri_offset += this->CellWriteOffset;
+
+    const Pass4TrimState state(
+      AxisToSum{}, this->PointDims, threadIndices, axis_mins, axis_maxs, edges);
+    if (!state.valid)
+    {
+      return;
+    }
+
+    const vtkm::Id3 pdims = this->PointDims;
+    const vtkm::Id3 increments = compute_incs3d(pdims);
+    vtkm::Id edgeIds[12];
+
+    auto edgeCase = getEdgeCase(edges, state.startPos, (state.axis_inc * state.left));
+    init_voxelIds(AxisToSum{}, this->PointWriteOffset, edgeCase, axis_sums, edgeIds);
+    for (vtkm::Id i = state.left; i < state.right; ++i) // run along the trimmed voxels
+    {
+      auto ijk = state.ijk;
+      ijk[AxisToSum::xindex] = i;
+      edgeCase = getEdgeCase(edges, state.startPos, (state.axis_inc * i));
+      vtkm::UInt8 numTris = data::GetNumberOfPrimitives(edgeCase);
+      if (numTris > 0)
+      {
+        //compute what the current cellId is
+        vtkm::Id cellId = compute_start(AxisToSum{}, ijk, pdims - vtkm::Id3{ 1, 1, 1 });
+
+        // Start by generating triangles for this case
+        generate_tris(cellId, edgeCase, numTris, edgeIds, cell_tri_offset, conn, inputCellIds);
+
+        // Now generate edgeIds and weights along voxel axes if needed. Remember to take
+        // boundary into account.
+        vtkm::UInt8 loc = static_cast<vtkm::UInt8>(
+          state.yzLoc | (i < 1 ? FlyingEdges3D::MinBoundary
+                               : (i >= (pdims[AxisToSum::xindex] - 2) ? FlyingEdges3D::MaxBoundary
+                                                                      : FlyingEdges3D::Interior)));
+        auto* edgeUses = data::GetEdgeUses(edgeCase);
+        if (loc != FlyingEdges3D::Interior || case_includes_axes(edgeUses))
+        {
+          this->Generate(loc,
+                         ijk,
+                         field,
+                         interpolatedEdgeIds,
+                         weights,
+                         points,
+                         normals,
+                         state.startPos,
+                         increments,
+                         (state.axis_inc * i),
+                         edgeUses,
+                         edgeIds);
+        }
+        advance_voxelIds(edgeUses, edgeIds);
+      }
+    }
+  }
+
+  //----------------------------------------------------------------------------
+  template <typename WholeDataField,
+            typename WholeIEdgeField,
+            typename WholeWeightField,
+            typename WholePointField,
+            typename WholeNormalField>
+  VTKM_EXEC inline void Generate(vtkm::UInt8 loc,
+                                 const vtkm::Id3& ijk,
+                                 const WholeDataField& field,
+                                 const WholeIEdgeField& interpolatedEdgeIds,
+                                 const WholeWeightField& weights,
+                                 const WholePointField& points,
+                                 const WholeNormalField& normals,
+                                 const vtkm::Id4& startPos,
+                                 const vtkm::Id3& incs,
+                                 vtkm::Id offset,
+                                 vtkm::UInt8 const* const edgeUses,
+                                 vtkm::Id* edgeIds) const
+  {
+    using AxisToSum = SumXAxis;
+
+    vtkm::Id2 pos(startPos[0] + offset, 0);
+    {
+      auto s0 = field.Get(pos[0]);
+      auto g0 = this->ComputeGradient(loc, ijk, incs, pos[0], field);
+
+      //EdgesUses 0,4,8 work for Y axis
+      if (edgeUses[0])
+      { // edgesUses[0] == i axes edge
+        auto writeIndex = edgeIds[0];
+        pos[1] = startPos[0] + offset + incs[AxisToSum::xindex];
+        auto s1 = field.Get(pos[1]);
+        T t = static_cast<T>((this->IsoValue - s0) / (s1 - s0));
+
+        interpolatedEdgeIds.Set(writeIndex, pos);
+        weights.Set(writeIndex, static_cast<vtkm::FloatDefault>(t));
+
+        auto ijk1 = ijk + vtkm::Id3{ 1, 0, 0 };
+        auto coord = this->InterpolateCoordinate(t, ijk, ijk1);
+        points.Set(writeIndex, coord);
+
+        //gradient generation
+        auto g1 = this->ComputeGradient(loc, ijk1, incs, pos[1], field);
+        g1 = g0 + (t * (g1 - g0));
+        normals.Set(writeIndex, vtkm::Normal(g1));
+      }
+      if (edgeUses[4])
+      { // edgesUses[4] == j axes edge
+        auto writeIndex = edgeIds[4];
+        pos[1] = startPos[1] + offset;
+        auto s1 = field.Get(pos[1]);
+        T t = static_cast<T>((this->IsoValue - s0) / (s1 - s0));
+
+        interpolatedEdgeIds.Set(writeIndex, pos);
+        weights.Set(writeIndex, static_cast<vtkm::FloatDefault>(t));
+
+        auto ijk1 = ijk + vtkm::Id3{ 0, 1, 0 };
+        auto coord = this->InterpolateCoordinate(t, ijk, ijk1);
+        points.Set(writeIndex, coord);
+
+        //gradient generation
+        auto g1 = this->ComputeGradient(loc, ijk1, incs, pos[1], field);
+        g1 = g0 + (t * (g1 - g0));
+        normals.Set(writeIndex, vtkm::Normal(g1));
+      }
+      if (edgeUses[8])
+      { // edgesUses[8] == k axes edge
+        auto writeIndex = edgeIds[8];
+        pos[1] = startPos[2] + offset;
+        auto s1 = field.Get(pos[1]);
+        T t = static_cast<T>((this->IsoValue - s0) / (s1 - s0));
+
+        interpolatedEdgeIds.Set(writeIndex, pos);
+        weights.Set(writeIndex, static_cast<vtkm::FloatDefault>(t));
+
+        auto ijk1 = ijk + vtkm::Id3{ 0, 0, 1 };
+        auto coord = this->InterpolateCoordinate(t, ijk, ijk1);
+        points.Set(writeIndex, coord);
+
+        //gradient generation
+        auto g1 = this->ComputeGradient(loc, ijk1, incs, pos[1], field);
+        g1 = g0 + (t * (g1 - g0));
+        normals.Set(writeIndex, vtkm::Normal(g1));
+      }
+    }
+    // On the boundary cells special work has to be done to cover the partial
+    // cell axes. These are boundary situations where the voxel axes is not
+    // fully formed. These situations occur on the +x,+y,+z volume
+    // boundaries. (The other cases fall through the default: case which is
+    // expected.)
+    //
+    // Note that loc is one of 27 regions in the volume, with (0,1,2)
+    // indicating (interior, min, max) along coordinate axes.
+
+    // clang-format off
+    switch (loc)
+    {
+      case 2:
+      case 6:
+      case 18:
+      case 22: //+x
+        this->InterpolateEdge(
+          loc, ijk, pos[0], incs, 5, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points, normals);
+        this->InterpolateEdge(
+          loc, ijk, pos[0], incs, 9, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points, normals);
+        break;
+      case 8:
+      case 9:
+      case 24:
+      case 25: //+y
+        this->InterpolateEdge(
+          loc, ijk, pos[0], incs, 1, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points, normals);
+        this->InterpolateEdge(
+          loc, ijk, pos[0], incs, 10, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points, normals);
+        break;
+      case 32:
+      case 33:
+      case 36:
+      case 37: //+z
+        this->InterpolateEdge(
+          loc, ijk, pos[0], incs, 2, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points, normals);
+        this->InterpolateEdge(
+          loc, ijk, pos[0], incs, 6, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points, normals);
+        break;
+      case 10:
+      case 26: //+x +y
+        this->InterpolateEdge(
+          loc, ijk, pos[0], incs, 1, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points, normals);
+        this->InterpolateEdge(
+          loc, ijk, pos[0], incs, 5, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points, normals);
+        this->InterpolateEdge(
+          loc, ijk, pos[0], incs, 9, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points, normals);
+        this->InterpolateEdge(
+          loc, ijk, pos[0], incs, 10, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points, normals);
+        this->InterpolateEdge(
+          loc, ijk, pos[0], incs, 11, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points, normals);
+        break;
+      case 34:
+      case 38: //+x +z
+        this->InterpolateEdge(
+          loc, ijk, pos[0], incs, 2, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points, normals);
+        this->InterpolateEdge(
+          loc, ijk, pos[0], incs, 5, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points, normals);
+        this->InterpolateEdge(
+          loc, ijk, pos[0], incs, 9, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points, normals);
+        this->InterpolateEdge(
+          loc, ijk, pos[0], incs, 6, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points, normals);
+        this->InterpolateEdge(
+          loc, ijk, pos[0], incs, 7, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points, normals);
+        break;
+      case 40:
+      case 41: //+y +z
+        this->InterpolateEdge(
+          loc, ijk, pos[0], incs, 1, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points, normals);
+        this->InterpolateEdge(
+          loc, ijk, pos[0], incs, 2, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points, normals);
+        this->InterpolateEdge(
+          loc, ijk, pos[0], incs, 3, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points, normals);
+        this->InterpolateEdge(
+          loc, ijk, pos[0], incs, 6, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points, normals);
+        this->InterpolateEdge(
+          loc, ijk, pos[0], incs, 10, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points, normals);
+        break;
+      case 42: //+x +y +z happens no more than once per volume
+        this->InterpolateEdge(
+          loc, ijk, pos[0], incs, 1, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points, normals);
+        this->InterpolateEdge(
+          loc, ijk, pos[0], incs, 2, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points, normals);
+        this->InterpolateEdge(
+          loc, ijk, pos[0], incs, 3, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points, normals);
+        this->InterpolateEdge(
+          loc, ijk, pos[0], incs, 5, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points, normals);
+        this->InterpolateEdge(
+          loc, ijk, pos[0], incs, 9, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points, normals);
+        this->InterpolateEdge(
+          loc, ijk, pos[0], incs, 10, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points, normals);
+        this->InterpolateEdge(
+          loc, ijk, pos[0], incs, 11, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points, normals);
+        this->InterpolateEdge(
+          loc, ijk, pos[0], incs, 6, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points, normals);
+        this->InterpolateEdge(
+          loc, ijk, pos[0], incs, 7, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points, normals);
+        break;
+      default: // interior, or -x,-y,-z boundaries
+        return;
+    }
+    // clang-format on
+  }
+
+  // Indicate whether voxel axes need processing for this case.
+  //----------------------------------------------------------------------------
+  template <typename WholeField,
+            typename WholeIEdgeField,
+            typename WholeWeightField,
+            typename WholePointField,
+            typename WholeNormalField>
+  VTKM_EXEC inline void InterpolateEdge(vtkm::UInt8 loc,
+                                        const vtkm::Id3& ijk,
+                                        vtkm::Id currentIdx,
+                                        const vtkm::Id3& incs,
+                                        vtkm::Id edgeNum,
+                                        vtkm::UInt8 const* const edgeUses,
+                                        vtkm::Id* edgeIds,
+                                        const WholeField& field,
+                                        const WholeIEdgeField& interpolatedEdgeIds,
+                                        const WholeWeightField& weights,
+                                        const WholePointField& points,
+                                        const WholeNormalField& normals) const
+  {
+    using AxisToSum = SumXAxis;
+
+    // if this edge is not used then get out
+    if (!edgeUses[edgeNum])
+    {
+      return;
+    }
+    const vtkm::Id writeIndex = edgeIds[edgeNum];
+
+    // build the edge information
+    vtkm::Vec<vtkm::UInt8, 2> verts = data::GetVertMap(edgeNum);
+
+    vtkm::Id3 offsets1 = data::GetVertOffsets(AxisToSum{}, verts[0]);
+    vtkm::Id3 offsets2 = data::GetVertOffsets(AxisToSum{}, verts[1]);
+
+    vtkm::Id2 iEdge(currentIdx + vtkm::Dot(offsets1, incs), currentIdx + vtkm::Dot(offsets2, incs));
+
+    interpolatedEdgeIds.Set(writeIndex, iEdge);
+
+    auto s0 = field.Get(iEdge[0]);
+    auto s1 = field.Get(iEdge[1]);
+    T t = static_cast<T>((this->IsoValue - s0) / (s1 - s0));
+    weights.Set(writeIndex, static_cast<vtkm::FloatDefault>(t));
+
+    auto coord = this->InterpolateCoordinate(t, ijk + offsets1, ijk + offsets2);
+    points.Set(writeIndex, coord);
+
+    auto g0 = this->ComputeGradient(loc, ijk + offsets1, incs, iEdge[0], field);
+    auto g1 = this->ComputeGradient(loc, ijk + offsets2, incs, iEdge[1], field);
+    g1 = g0 + (t * (g1 - g0));
+    normals.Set(writeIndex, vtkm::Normal(g1));
+  }
+
+  //----------------------------------------------------------------------------
+  inline VTKM_EXEC vtkm::Vec3f InterpolateCoordinate(T t,
+                                                     const vtkm::Id3& ijk0,
+                                                     const vtkm::Id3& ijk1) const
+  {
+    return vtkm::Vec3f(
+      this->Origin[0] +
+        this->Spacing[0] * static_cast<vtkm::FloatDefault>(ijk0[0] + t * (ijk1[0] - ijk0[0])),
+      this->Origin[1] +
+        this->Spacing[1] * static_cast<vtkm::FloatDefault>(ijk0[1] + t * (ijk1[1] - ijk0[1])),
+      this->Origin[2] +
+        this->Spacing[2] * static_cast<vtkm::FloatDefault>(ijk0[2] + t * (ijk1[2] - ijk0[2])));
+  }
+
+  //----------------------------------------------------------------------------
+  template <typename WholeDataField>
+  VTKM_EXEC vtkm::Vec3f ComputeGradient(vtkm::UInt8 loc,
+                                        const vtkm::Id3& ijk,
+                                        const vtkm::Id3& incs,
+                                        vtkm::Id pos,
+                                        const WholeDataField& field) const
+  {
+    if (loc == FlyingEdges3D::Interior)
+    {
+      vtkm::Vec3f g = {
+        static_cast<vtkm::FloatDefault>(field.Get(pos + incs[0]) - field.Get(pos - incs[0])) * 0.5f,
+        static_cast<vtkm::FloatDefault>(field.Get(pos + incs[1]) - field.Get(pos - incs[1])) * 0.5f,
+        static_cast<vtkm::FloatDefault>(field.Get(pos + incs[2]) - field.Get(pos - incs[2])) * 0.5f
+      };
+      return g;
+    }
+
+    //We are on some boundary edge
+    auto s = field.Get(pos);
+    vtkm::Vec3f g;
+    for (int i = 0; i < 3; ++i)
+    {
+      if (ijk[i] == 0)
+      {
+        g[i] = static_cast<vtkm::FloatDefault>(field.Get(pos + incs[i]) - s);
+      }
+      else if (ijk[i] >= (this->PointDims[i] - 1))
+      {
+        g[i] = static_cast<vtkm::FloatDefault>(s - field.Get(pos - incs[i]));
+      }
+      else
+      {
+        g[i] =
+          static_cast<vtkm::FloatDefault>(field.Get(pos + incs[i]) - field.Get(pos - incs[i])) *
+          0.5f;
+      }
+    }
+
+    return g;
+  }
+};
+}
+}
+}
+#endif

vtkm/worklet/contour/FlyingEdgesPass4Y.h

@@ -0,0 +1,425 @@
+
+//============================================================================
+//  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.
+//============================================================================
+
+
+#ifndef vtk_m_worklet_contour_flyingedges_pass4y_h
+#define vtk_m_worklet_contour_flyingedges_pass4y_h
+
+
+#include <vtkm/worklet/contour/FlyingEdgesHelpers.h>
+#include <vtkm/worklet/contour/FlyingEdgesTables.h>
+
+#include <vtkm/VectorAnalysis.h>
+#include <vtkm/worklet/gradient/StructuredPointGradient.h>
+
+namespace vtkm
+{
+namespace worklet
+{
+namespace flying_edges
+{
+
+
+template <typename T>
+struct ComputePass4Y : public vtkm::worklet::WorkletVisitCellsWithPoints
+{
+
+  vtkm::Id3 PointDims;
+
+  T IsoValue;
+
+  vtkm::Id CellWriteOffset;
+  vtkm::Id PointWriteOffset;
+
+  ComputePass4Y() {}
+  ComputePass4Y(T value,
+                const vtkm::Id3& pdims,
+                vtkm::Id multiContourCellOffset,
+                vtkm::Id multiContourPointOffset)
+    : PointDims(pdims)
+    , IsoValue(value)
+    , CellWriteOffset(multiContourCellOffset)
+    , PointWriteOffset(multiContourPointOffset)
+  {
+  }
+
+  using ControlSignature = void(CellSetIn,
+                                FieldInPoint axis_sums,
+                                FieldInPoint axis_mins,
+                                FieldInPoint axis_maxs,
+                                WholeArrayIn cell_tri_count,
+                                WholeArrayIn edgeData,
+                                WholeArrayIn data,
+                                WholeArrayOut connectivity,
+                                WholeArrayOut edgeIds,
+                                WholeArrayOut weights,
+                                WholeArrayOut inputCellIds);
+  using ExecutionSignature =
+    void(ThreadIndices, _2, _3, _4, _5, _6, _7, _8, _9, _10, _11, WorkIndex);
+
+  template <typename ThreadIndices,
+            typename FieldInPointId3,
+            typename FieldInPointId,
+            typename WholeTriField,
+            typename WholeEdgeField,
+            typename WholeDataField,
+            typename WholeConnField,
+            typename WholeEdgeIdField,
+            typename WholeWeightField,
+            typename WholeCellIdField>
+  VTKM_EXEC void operator()(const ThreadIndices& threadIndices,
+                            const FieldInPointId3& axis_sums,
+                            const FieldInPointId& axis_mins,
+                            const FieldInPointId& axis_maxs,
+                            const WholeTriField& cellTriCount,
+                            const WholeEdgeField& edges,
+                            const WholeDataField& field,
+                            const WholeConnField& conn,
+                            const WholeEdgeIdField& interpolatedEdgeIds,
+                            const WholeWeightField& weights,
+                            const WholeCellIdField& inputCellIds,
+                            vtkm::Id oidx) const
+  {
+    using AxisToSum = SumYAxis;
+
+    //This works as cellTriCount was computed with ScanExtended
+    //and therefore has one more entry than the number of cells
+    vtkm::Id cell_tri_offset = cellTriCount.Get(oidx);
+    vtkm::Id next_tri_offset = cellTriCount.Get(oidx + 1);
+    if (cell_tri_offset == next_tri_offset)
+    { //we produce nothing
+      return;
+    }
+    cell_tri_offset += this->CellWriteOffset;
+
+    const Pass4TrimState state(
+      AxisToSum{}, this->PointDims, threadIndices, axis_mins, axis_maxs, edges);
+    if (!state.valid)
+    {
+      return;
+    }
+
+    const vtkm::Id3 pdims = this->PointDims;
+    const vtkm::Id3 increments = compute_incs3d(pdims);
+    vtkm::Id edgeIds[12];
+
+    auto edgeCase = getEdgeCase(edges, state.startPos, (state.axis_inc * state.left));
+    init_voxelIds(AxisToSum{}, this->PointWriteOffset, edgeCase, axis_sums, edgeIds);
+    for (vtkm::Id i = state.left; i < state.right; ++i) // run along the trimmed voxels
+    {
+      auto ijk = state.ijk;
+      ijk[AxisToSum::xindex] = i;
+      edgeCase = getEdgeCase(edges, state.startPos, (state.axis_inc * i));
+      vtkm::UInt8 numTris = data::GetNumberOfPrimitives(edgeCase);
+      if (numTris > 0)
+      {
+        //compute what the current cellId is
+        vtkm::Id cellId = compute_start(AxisToSum{}, ijk, pdims - vtkm::Id3{ 1, 1, 1 });
+
+        // Start by generating triangles for this case
+        generate_tris(cellId, edgeCase, numTris, edgeIds, cell_tri_offset, conn, inputCellIds);
+
+        // Now generate edgeIds and weights along voxel axes if needed. Remember to take
+        // boundary into account.
+        vtkm::UInt8 loc = static_cast<vtkm::UInt8>(
+          state.yzLoc | (i < 1 ? FlyingEdges3D::MinBoundary
+                               : (i >= (pdims[AxisToSum::xindex] - 2) ? FlyingEdges3D::MaxBoundary
+                                                                      : FlyingEdges3D::Interior)));
+        auto* edgeUses = data::GetEdgeUses(edgeCase);
+        if (loc != FlyingEdges3D::Interior || case_includes_axes(edgeUses))
+        {
+          this->Generate(loc,
+                         field,
+                         interpolatedEdgeIds,
+                         weights,
+                         state.startPos,
+                         increments,
+                         (state.axis_inc * i),
+                         edgeUses,
+                         edgeIds);
+        }
+        advance_voxelIds(edgeUses, edgeIds);
+      }
+    }
+  }
+
+  //----------------------------------------------------------------------------
+  template <typename WholeDataField, typename WholeIEdgeField, typename WholeWeightField>
+  VTKM_EXEC inline void Generate(vtkm::UInt8 loc,
+                                 const WholeDataField& field,
+                                 const WholeIEdgeField& interpolatedEdgeIds,
+                                 const WholeWeightField& weights,
+                                 const vtkm::Id4& startPos,
+                                 const vtkm::Id3& incs,
+                                 vtkm::Id offset,
+                                 vtkm::UInt8 const* const edgeUses,
+                                 vtkm::Id* edgeIds) const
+  {
+    using AxisToSum = SumYAxis;
+
+    vtkm::Id2 pos(startPos[0] + offset, 0);
+    {
+      auto s0 = field.Get(pos[0]);
+
+      //EdgesUses 0,4,8 work for Y axis
+      if (edgeUses[0])
+      { // edgesUses[0] == i axes edge
+        auto writeIndex = edgeIds[0];
+        pos[1] = startPos[0] + offset + incs[AxisToSum::xindex];
+        auto s1 = field.Get(pos[1]);
+        T t = static_cast<T>((this->IsoValue - s0) / (s1 - s0));
+
+        interpolatedEdgeIds.Set(writeIndex, pos);
+        weights.Set(writeIndex, static_cast<vtkm::FloatDefault>(t));
+      }
+      if (edgeUses[4])
+      { // edgesUses[4] == j axes edge
+        auto writeIndex = edgeIds[4];
+        pos[1] = startPos[1] + offset;
+        auto s1 = field.Get(pos[1]);
+        T t = static_cast<T>((this->IsoValue - s0) / (s1 - s0));
+
+        interpolatedEdgeIds.Set(writeIndex, pos);
+        weights.Set(writeIndex, static_cast<vtkm::FloatDefault>(t));
+      }
+      if (edgeUses[8])
+      { // edgesUses[8] == k axes edge
+        auto writeIndex = edgeIds[8];
+        pos[1] = startPos[2] + offset;
+        auto s1 = field.Get(pos[1]);
+        T t = static_cast<T>((this->IsoValue - s0) / (s1 - s0));
+
+        interpolatedEdgeIds.Set(writeIndex, pos);
+        weights.Set(writeIndex, static_cast<vtkm::FloatDefault>(t));
+      }
+    }
+    // On the boundary cells special work has to be done to cover the partial
+    // cell axes. These are boundary situations where the voxel axes is not
+    // fully formed. These situations occur on the +x,+y,+z volume
+    // boundaries. (The other cases fall through the default: case which is
+    // expected.)
+    //
+    // Note that loc is one of 27 regions in the volume, with (0,1,2)
+    // indicating (interior, min, max) along coordinate axes.
+    switch (loc)
+    {
+      case 2:
+      case 6:
+      case 18:
+      case 22: //+x
+        this->InterpolateEdge(
+          pos[0], incs, 5, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
+        this->InterpolateEdge(
+          pos[0], incs, 9, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
+        break;
+      case 8:
+      case 9:
+      case 24:
+      case 25: //+y
+        this->InterpolateEdge(
+          pos[0], incs, 1, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
+        this->InterpolateEdge(
+          pos[0], incs, 10, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
+        break;
+      case 32:
+      case 33:
+      case 36:
+      case 37: //+z
+        this->InterpolateEdge(
+          pos[0], incs, 2, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
+        this->InterpolateEdge(
+          pos[0], incs, 6, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
+        break;
+      case 10:
+      case 26: //+x +y
+        this->InterpolateEdge(
+          pos[0], incs, 1, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
+        this->InterpolateEdge(
+          pos[0], incs, 5, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
+        this->InterpolateEdge(
+          pos[0], incs, 9, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
+        this->InterpolateEdge(
+          pos[0], incs, 10, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
+        this->InterpolateEdge(
+          pos[0], incs, 11, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
+        break;
+      case 34:
+      case 38: //+x +z
+        this->InterpolateEdge(
+          pos[0], incs, 2, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
+        this->InterpolateEdge(
+          pos[0], incs, 5, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
+        this->InterpolateEdge(
+          pos[0], incs, 9, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
+        this->InterpolateEdge(
+          pos[0], incs, 6, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
+        this->InterpolateEdge(
+          pos[0], incs, 7, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
+        break;
+      case 40:
+      case 41: //+y +z
+        this->InterpolateEdge(
+          pos[0], incs, 1, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
+        this->InterpolateEdge(
+          pos[0], incs, 2, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
+        this->InterpolateEdge(
+          pos[0], incs, 3, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
+        this->InterpolateEdge(
+          pos[0], incs, 6, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
+        this->InterpolateEdge(
+          pos[0], incs, 10, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
+        break;
+      case 42: //+x +y +z happens no more than once per volume
+        this->InterpolateEdge(
+          pos[0], incs, 1, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
+        this->InterpolateEdge(
+          pos[0], incs, 2, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
+        this->InterpolateEdge(
+          pos[0], incs, 3, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
+        this->InterpolateEdge(
+          pos[0], incs, 5, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
+        this->InterpolateEdge(
+          pos[0], incs, 9, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
+        this->InterpolateEdge(
+          pos[0], incs, 10, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
+        this->InterpolateEdge(
+          pos[0], incs, 11, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
+        this->InterpolateEdge(
+          pos[0], incs, 6, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
+        this->InterpolateEdge(
+          pos[0], incs, 7, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
+        break;
+      default: // interior, or -x,-y,-z boundaries
+        return;
+    }
+  }
+
+  // Indicate whether voxel axes need processing for this case.
+  //----------------------------------------------------------------------------
+  template <typename WholeField, typename WholeIEdgeField, typename WholeWeightField>
+  VTKM_EXEC inline void InterpolateEdge(vtkm::Id currentIdx,
+                                        const vtkm::Id3& incs,
+                                        vtkm::Id edgeNum,
+                                        vtkm::UInt8 const* const edgeUses,
+                                        vtkm::Id* edgeIds,
+                                        const WholeField& field,
+                                        const WholeIEdgeField& interpolatedEdgeIds,
+                                        const WholeWeightField& weights) const
+  {
+    using AxisToSum = SumYAxis;
+
+    // if this edge is not used then get out
+    if (!edgeUses[edgeNum])
+    {
+      return;
+    }
+    const vtkm::Id writeIndex = edgeIds[edgeNum];
+
+    // build the edge information
+    vtkm::Vec<vtkm::UInt8, 2> verts = data::GetVertMap(edgeNum);
+
+    vtkm::Id3 offsets1 = data::GetVertOffsets(AxisToSum{}, verts[0]);
+    vtkm::Id3 offsets2 = data::GetVertOffsets(AxisToSum{}, verts[1]);
+
+    vtkm::Id2 iEdge(currentIdx + vtkm::Dot(offsets1, incs), currentIdx + vtkm::Dot(offsets2, incs));
+
+    interpolatedEdgeIds.Set(writeIndex, iEdge);
+
+    auto s0 = field.Get(iEdge[0]);
+    auto s1 = field.Get(iEdge[1]);
+    T t = static_cast<T>((this->IsoValue - s0) / (s1 - s0));
+    weights.Set(writeIndex, static_cast<vtkm::FloatDefault>(t));
+  }
+};
+
+template <typename T>
+struct ComputePass5Y : public vtkm::worklet::WorkletMapField
+{
+
+  vtkm::internal::ArrayPortalUniformPointCoordinates Coordinates;
+  vtkm::Id NormalWriteOffset;
+
+  ComputePass5Y() {}
+  ComputePass5Y(const vtkm::Id3& pdims,
+                const vtkm::Vec3f& origin,
+                const vtkm::Vec3f& spacing,
+                vtkm::Id normalWriteOffset,
+                bool generateNormals)
+    : Coordinates(pdims, origin, spacing)
+    , NormalWriteOffset(normalWriteOffset)
+  {
+    if (!generateNormals)
+    {
+      this->NormalWriteOffset = -1;
+    }
+  }
+
+  using ControlSignature = void(FieldIn interpEdgeIds,
+                                FieldIn interpWeight,
+                                FieldOut points,
+                                WholeArrayIn field,
+                                WholeArrayOut normals);
+  using ExecutionSignature = void(_1, _2, _3, _4, _5, WorkIndex);
+
+  template <typename PT, typename WholeInputField, typename WholeNormalField>
+  VTKM_EXEC void operator()(const vtkm::Id2& interpEdgeIds,
+                            vtkm::FloatDefault weight,
+                            vtkm::Vec<PT, 3>& outPoint,
+                            const WholeInputField& field,
+                            WholeNormalField& normals,
+                            vtkm::Id oidx) const
+  {
+    {
+      vtkm::Vec3f point1 = this->Coordinates.Get(interpEdgeIds[0]);
+      vtkm::Vec3f point2 = this->Coordinates.Get(interpEdgeIds[1]);
+      outPoint = vtkm::Lerp(point1, point2, weight);
+    }
+
+    //NormalWriteOffset of -1 means no normals
+    if (this->NormalWriteOffset >= 0)
+    {
+      vtkm::Vec<T, 3> g0, g1;
+      const vtkm::Id3& dims = this->Coordinates.GetDimensions();
+      vtkm::Id3 ijk{ interpEdgeIds[0] % dims[0],
+                     (interpEdgeIds[0] / dims[0]) % dims[1],
+                     interpEdgeIds[0] / (dims[0] * dims[1]) };
+
+      vtkm::worklet::gradient::StructuredPointGradient gradient;
+      vtkm::exec::BoundaryState boundary(ijk, dims);
+      vtkm::exec::FieldNeighborhood<vtkm::internal::ArrayPortalUniformPointCoordinates>
+        coord_neighborhood(this->Coordinates, boundary);
+
+      vtkm::exec::FieldNeighborhood<WholeInputField> field_neighborhood(field, boundary);
+
+
+      //compute the gradient at point 1
+      gradient(boundary, coord_neighborhood, field_neighborhood, g0);
+
+      //compute the gradient at point 2. This optimization can be optimized
+      boundary.IJK = vtkm::Id3{ interpEdgeIds[1] % dims[0],
+                                (interpEdgeIds[1] / dims[0]) % dims[1],
+                                interpEdgeIds[1] / (dims[0] * dims[1]) };
+      gradient(boundary, coord_neighborhood, field_neighborhood, g1);
+
+      vtkm::Vec3f n = vtkm::Lerp(g0, g1, weight);
+      const auto mag2 = vtkm::MagnitudeSquared(n);
+      if (mag2 > 0.)
+      {
+        n = n * vtkm::RSqrt(mag2);
+      }
+      normals.Set(this->NormalWriteOffset + oidx, n);
+    }
+  }
+};
+}
+}
+}
+#endif

vtkm/worklet/contour/FlyingEdgesPass5.h

@@ -1,106 +0,0 @@
-
-//============================================================================
-//  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.
-//============================================================================
-
-
-#ifndef vtk_m_worklet_contour_flyingedges_pass5_h
-#define vtk_m_worklet_contour_flyingedges_pass5_h
-
-
-#include <vtkm/worklet/contour/FlyingEdgesHelpers.h>
-#include <vtkm/worklet/contour/FlyingEdgesTables.h>
-
-#include <vtkm/VectorAnalysis.h>
-#include <vtkm/worklet/gradient/StructuredPointGradient.h>
-
-namespace vtkm
-{
-namespace worklet
-{
-namespace flying_edges
-{
-
-template <typename T>
-struct ComputePass5 : public vtkm::worklet::WorkletMapField
-{
-
-  vtkm::internal::ArrayPortalUniformPointCoordinates Coordinates;
-  bool GenerateNormals;
-
-  ComputePass5() {}
-  ComputePass5(const vtkm::Id3& pdims,
-               const vtkm::Vec3f& origin,
-               const vtkm::Vec3f& spacing,
-               bool generateNormals)
-    : Coordinates(pdims, origin, spacing)
-    , GenerateNormals(generateNormals)
-  {
-  }
-
-  using ControlSignature = void(FieldIn interpEdgeIds,
-                                FieldIn interpWeight,
-                                FieldOut points,
-                                WholeArrayIn field,
-                                WholeArrayOut normals);
-  using ExecutionSignature = void(_1, _2, _3, _4, _5, WorkIndex);
-
-  template <typename PT, typename WholeInputField, typename WholeNormalField>
-  VTKM_EXEC void operator()(const vtkm::Id2& interpEdgeIds,
-                            vtkm::FloatDefault weight,
-                            vtkm::Vec<PT, 3>& outPoint,
-                            const WholeInputField& field,
-                            WholeNormalField& normals,
-                            vtkm::Id oidx) const
-  {
-    {
-      vtkm::Vec3f point1 = this->Coordinates.Get(interpEdgeIds[0]);
-      vtkm::Vec3f point2 = this->Coordinates.Get(interpEdgeIds[1]);
-      outPoint = vtkm::Lerp(point1, point2, weight);
-    }
-
-    if (this->GenerateNormals)
-    {
-      vtkm::Vec<T, 3> g0, g1;
-      const vtkm::Id3& dims = this->Coordinates.GetDimensions();
-      vtkm::Id3 ijk{ interpEdgeIds[0] % dims[0],
-                     (interpEdgeIds[0] / dims[0]) % dims[1],
-                     interpEdgeIds[0] / (dims[0] * dims[1]) };
-
-      vtkm::worklet::gradient::StructuredPointGradient gradient;
-      vtkm::exec::BoundaryState boundary(ijk, dims);
-      vtkm::exec::FieldNeighborhood<vtkm::internal::ArrayPortalUniformPointCoordinates>
-        coord_neighborhood(this->Coordinates, boundary);
-
-      vtkm::exec::FieldNeighborhood<WholeInputField> field_neighborhood(field, boundary);
-
-
-      //compute the gradient at point 1
-      gradient(boundary, coord_neighborhood, field_neighborhood, g0);
-
-      //compute the gradient at point 2. This optimization can be optimized
-      boundary.IJK = vtkm::Id3{ interpEdgeIds[1] % dims[0],
-                                (interpEdgeIds[1] / dims[0]) % dims[1],
-                                interpEdgeIds[1] / (dims[0] * dims[1]) };
-      gradient(boundary, coord_neighborhood, field_neighborhood, g1);
-
-      vtkm::Vec3f n = vtkm::Lerp(g0, g1, weight);
-      const auto mag2 = vtkm::MagnitudeSquared(n);
-      if (mag2 > 0.)
-      {
-        n = n * vtkm::RSqrt(mag2);
-      }
-      normals.Set(oidx, n);
-    }
-  }
-};
-}
-}
-}
-#endif

vtkm/worklet/gradient/StructuredPointGradient.h

@@ -45,15 +45,19 @@ struct StructuredPointGradient : public vtkm::worklet::WorkletPointNeighborhood
     using OT = typename GradientOutType::ComponentType;
 
     vtkm::Vec<CT, 3> xi, eta, zeta;
-    this->Jacobian(inputPoints, boundary, xi, eta, zeta); //store the metrics in xi,eta,zeta
+    vtkm::Vec<bool, 3> onBoundary{ !boundary.IsRadiusInXBoundary(1),
+                                   !boundary.IsRadiusInYBoundary(1),
+                                   !boundary.IsRadiusInZBoundary(1) };
+
+    this->Jacobian(inputPoints, onBoundary, xi, eta, zeta); //store the metrics in xi,eta,zeta
 
     auto dxi = inputField.Get(1, 0, 0) - inputField.Get(-1, 0, 0);
     auto deta = inputField.Get(0, 1, 0) - inputField.Get(0, -1, 0);
     auto dzeta = inputField.Get(0, 0, 1) - inputField.Get(0, 0, -1);
 
-    dxi = (boundary.IsRadiusInXBoundary(1) ? dxi * 0.5f : dxi);
+    dxi = (onBoundary[0] ? dxi : dxi * 0.5f);
-    deta = (boundary.IsRadiusInYBoundary(1) ? deta * 0.5f : deta);
+    deta = (onBoundary[1] ? deta : deta * 0.5f);
-    dzeta = (boundary.IsRadiusInZBoundary(1) ? dzeta * 0.5f : dzeta);
+    dzeta = (onBoundary[2] ? dzeta : dzeta * 0.5f);
 
     outputGradient[0] = static_cast<OT>(xi[0] * dxi + eta[0] * deta + zeta[0] * dzeta);
     outputGradient[1] = static_cast<OT>(xi[1] * dxi + eta[1] * deta + zeta[1] * dzeta);
@@ -75,24 +79,44 @@ struct StructuredPointGradient : public vtkm::worklet::WorkletPointNeighborhood
     using CoordType = typename PointsIn::ValueType;
     using OT = typename GradientOutType::ComponentType;
 
-
     CoordType r = inputPoints.Portal.GetSpacing();
 
-    r[0] = (boundary.IsRadiusInXBoundary(1) ? r[0] * 0.5f : r[0]);
-    r[1] = (boundary.IsRadiusInYBoundary(1) ? r[1] * 0.5f : r[1]);
-    r[2] = (boundary.IsRadiusInZBoundary(1) ? r[2] * 0.5f : r[2]);
-
-    auto dx = inputField.Get(1, 0, 0) - inputField.Get(-1, 0, 0);
-    auto dy = inputField.Get(0, 1, 0) - inputField.Get(0, -1, 0);
-    auto dz = inputField.Get(0, 0, 1) - inputField.Get(0, 0, -1);
-
 #if (defined(VTKM_CUDA) && defined(VTKM_GCC))
 #pragma GCC diagnostic push
 #pragma GCC diagnostic ignored "-Wconversion"
 #endif
-    outputGradient[0] = static_cast<OT>(dx * r[0]);
+    if (boundary.IsRadiusInXBoundary(1))
-    outputGradient[1] = static_cast<OT>(dy * r[1]);
+    {
-    outputGradient[2] = static_cast<OT>(dz * r[2]);
+      auto dx = inputField.GetUnchecked(1, 0, 0) - inputField.GetUnchecked(-1, 0, 0);
+      outputGradient[0] = static_cast<OT>(dx * (r[0] * 0.5f));
+    }
+    else
+    {
+      auto dx = inputField.Get(1, 0, 0) - inputField.Get(-1, 0, 0);
+      outputGradient[0] = static_cast<OT>(dx * r[0]);
+    }
+
+    if (boundary.IsRadiusInYBoundary(1))
+    {
+      auto dy = inputField.GetUnchecked(0, 1, 0) - inputField.GetUnchecked(0, -1, 0);
+      outputGradient[1] = static_cast<OT>(dy * r[1] * 0.5f);
+    }
+    else
+    {
+      auto dy = inputField.Get(0, 1, 0) - inputField.Get(0, -1, 0);
+      outputGradient[1] = static_cast<OT>(dy * (r[1]));
+    }
+
+    if (boundary.IsRadiusInZBoundary(1))
+    {
+      auto dz = inputField.GetUnchecked(0, 0, 1) - inputField.GetUnchecked(0, 0, -1);
+      outputGradient[2] = static_cast<OT>(dz * r[2] * 0.5f);
+    }
+    else
+    {
+      auto dz = inputField.Get(0, 0, 1) - inputField.Get(0, 0, -1);
+      outputGradient[2] = static_cast<OT>(dz * (r[2]));
+    }
 #if (defined(VTKM_CUDA) && defined(VTKM_GCC))
 #pragma GCC diagnostic pop
 #endif
@@ -103,20 +127,41 @@ struct StructuredPointGradient : public vtkm::worklet::WorkletPointNeighborhood
   //will be float,3 even when T is a 3 component field
   template <typename PointsIn, typename CT>
   VTKM_EXEC void Jacobian(const PointsIn& inputPoints,
-                          const vtkm::exec::BoundaryState& boundary,
+                          const vtkm::Vec<bool, 3>& onBoundary,
                           vtkm::Vec<CT, 3>& m_xi,
                           vtkm::Vec<CT, 3>& m_eta,
                           vtkm::Vec<CT, 3>& m_zeta) const
   {
     using CoordType = typename PointsIn::ValueType;
+    CoordType xi, eta, zeta;
 
-    CoordType xi = inputPoints.Get(1, 0, 0) - inputPoints.Get(-1, 0, 0);
-    CoordType eta = inputPoints.Get(0, 1, 0) - inputPoints.Get(0, -1, 0);
-    CoordType zeta = inputPoints.Get(0, 0, 1) - inputPoints.Get(0, 0, -1);
 
-    xi = (boundary.IsRadiusInXBoundary(1) ? xi * 0.5f : xi);
+    if (onBoundary[0])
-    eta = (boundary.IsRadiusInYBoundary(1) ? eta * 0.5f : eta);
+    {
-    zeta = (boundary.IsRadiusInZBoundary(1) ? zeta * 0.5f : zeta);
+      xi = (inputPoints.Get(1, 0, 0) - inputPoints.Get(-1, 0, 0)) * 0.5f;
+    }
+    else
+    {
+      xi = inputPoints.GetUnchecked(1, 0, 0) - inputPoints.GetUnchecked(-1, 0, 0);
+    }
+
+    if (onBoundary[1])
+    {
+      eta = (inputPoints.Get(0, 1, 0) - inputPoints.Get(0, -1, 0)) * 0.5f;
+    }
+    else
+    {
+      eta = inputPoints.GetUnchecked(0, 1, 0) - inputPoints.GetUnchecked(0, -1, 0);
+    }
+
+    if (onBoundary[2])
+    {
+      zeta = (inputPoints.Get(0, 0, 1) - inputPoints.Get(0, 0, -1)) * 0.5f;
+    }
+    else
+    {
+      zeta = inputPoints.GetUnchecked(0, 0, 1) - inputPoints.GetUnchecked(0, 0, -1);
+    }
 
     CT aj = xi[0] * eta[1] * zeta[2] + xi[1] * eta[2] * zeta[0] + xi[2] * eta[0] * zeta[1] -
       xi[2] * eta[1] * zeta[0] - xi[1] * eta[0] * zeta[2] - xi[0] * eta[2] * zeta[1];

vtkm/worklet/testing/UnitTestOrientNormals.cxx

@@ -33,6 +33,7 @@
 #include <vtkm/cont/CoordinateSystem.h>
 #include <vtkm/cont/DataSet.h>
 
+#include <vtkm/filter/CleanGrid.h>
 #include <vtkm/filter/Contour.h>
 #include <vtkm/filter/PolicyBase.h>
 #include <vtkm/filter/SurfaceNormals.h>
@@ -56,14 +57,16 @@ vtkm::cont::DataSet CreateDataSet(bool pointNormals, bool cellNormals)
   wavelet.SetMagnitude({ 5 });
   auto dataSet = wavelet.Execute();
 
-  // Cut a contour
+  vtkm::filter::CleanGrid toGrid;
+
+  // unstructured grid contour
   vtkm::filter::Contour contour;
   contour.SetActiveField("scalars", vtkm::cont::Field::Association::POINTS);
   contour.SetNumberOfIsoValues(1);
   contour.SetIsoValue(192);
   contour.SetMergeDuplicatePoints(true);
   contour.SetGenerateNormals(false);
-  dataSet = contour.Execute(dataSet);
+  dataSet = contour.Execute(toGrid.Execute(dataSet));
 
   vtkm::filter::SurfaceNormals normals;
   normals.SetGeneratePointNormals(pointNormals);