Implement Flying Edges for all structured CellSets

In addition to  using uniform coordinates, the ContourFlyingEdges filter can now process any type of coordinate system, making the  filter use Flying Edges in more cases

docs/changelog/flying-edges-structured.md

@@ -0,0 +1,5 @@
+# Implement Flying Edges for structured cellsets with rectilinear and curvilinear coordinates
+
+When Flying Edges was introduced to compute contours of a 3D structured cellset, it could only process uniform coordinates. This limitation is now lifted : an alternative interpolation function can be used in the fourth pass of the algorithm in order to support rectilinear and curvilinear coordinate systems.
+
+Accordingly, the `Contour` filter now calls `ContourFlyingEdges` instead of `ContourMarchingCells` for these newly supported cases.

vtkm/filter/clean_grid/testing/UnitTestCleanGrid.cxx

@@ -13,7 +13,7 @@
 #include <vtkm/cont/ArrayCopy.h>
 #include <vtkm/cont/testing/MakeTestDataSet.h>
 #include <vtkm/cont/testing/Testing.h>
-#include <vtkm/filter/contour/Contour.h>
+#include <vtkm/filter/contour/ContourMarchingCells.h>
 
 namespace
 {
@@ -72,13 +72,7 @@ void TestPointMerging()
   vtkm::cont::testing::MakeTestDataSet makeDataSet;
   vtkm::cont::DataSet baseData = makeDataSet.Make3DUniformDataSet3(vtkm::Id3(4, 4, 4));
 
-  //Convert the baseData implicit points to explicit points, since the contour
+  vtkm::filter::contour::ContourMarchingCells marchingCubes;
-  //filter for uniform data always does point merging
-  vtkm::cont::ArrayHandle<vtkm::Vec3f> newcoords;
-  vtkm::cont::ArrayCopy(baseData.GetCoordinateSystem().GetData(), newcoords);
-  baseData.AddPointField(baseData.GetCoordinateSystemName(), newcoords);
-
-  vtkm::filter::contour::Contour marchingCubes;
   marchingCubes.SetIsoValue(0.05);
   marchingCubes.SetMergeDuplicatePoints(false);
   marchingCubes.SetActiveField("pointvar");

vtkm/filter/contour/Contour.cxx

@@ -34,11 +34,8 @@ vtkm::cont::DataSet Contour::DoExecute(const vtkm::cont::DataSet& inDataSet)
   auto inCoords = inDataSet.GetCoordinateSystem(this->GetActiveCoordinateSystemIndex()).GetData();
   std::unique_ptr<vtkm::filter::contour::AbstractContour> implementation;
 
-  // For now, Flying Edges is only used for 3D Structured CellSets,
+  // Flying Edges is only used for 3D Structured CellSets
-  // using Uniform coordinates.
+  if (inCellSet.template IsType<vtkm::cont::CellSetStructured<3>>())
-  if (inCellSet.template IsType<vtkm::cont::CellSetStructured<3>>() &&
-      inCoords.template IsType<
-        vtkm::cont::ArrayHandle<vtkm::Vec3f, vtkm::cont::StorageTagUniformPoints>>())
   {
     VTKM_LOG_S(vtkm::cont::LogLevel::Info, "Using flying edges");
     implementation.reset(new vtkm::filter::contour::ContourFlyingEdges);

vtkm/filter/contour/ContourFlyingEdges.cxx

@@ -44,13 +44,10 @@ vtkm::cont::DataSet ContourFlyingEdges::DoExecute(const vtkm::cont::DataSet& inD
   const vtkm::cont::CoordinateSystem& inCoords =
     inDataSet.GetCoordinateSystem(this->GetActiveCoordinateSystemIndex());
 
-  if (!inCellSet.template IsType<vtkm::cont::CellSetStructured<3>>() ||
+  if (!inCellSet.template IsType<vtkm::cont::CellSetStructured<3>>())
-      !inCoords.GetData()
-         .template IsType<
-           vtkm::cont::ArrayHandle<vtkm::Vec3f, vtkm::cont::StorageTagUniformPoints>>())
   {
     throw vtkm::cont::ErrorFilterExecution("This filter is only available for 3-Dimensional "
-                                           "Structured Cell Sets using uniform point coordinates.");
+                                           "Structured Cell Sets");
   }
 
   // Get the CellSet's known dynamic type
@@ -74,22 +71,11 @@ vtkm::cont::DataSet ContourFlyingEdges::DoExecute(const vtkm::cont::DataSet& inD
 
     if (this->GenerateNormals && !this->GetComputeFastNormals())
     {
-      outputCells = worklet.Run(
+      outputCells = worklet.Run(ivalues, inputCells, inCoords, concrete, vertices, normals);
-        ivalues,
-        inputCells,
-        inCoords.GetData().AsArrayHandle<vtkm::cont::ArrayHandleUniformPointCoordinates>(),
-        concrete,
-        vertices,
-        normals);
     }
     else
     {
-      outputCells = worklet.Run(
+      outputCells = worklet.Run(ivalues, inputCells, inCoords, concrete, vertices);
-        ivalues,
-        inputCells,
-        inCoords.GetData().AsArrayHandle<vtkm::cont::ArrayHandleUniformPointCoordinates>(),
-        concrete,
-        vertices);
     }
   };
 

vtkm/filter/contour/testing/UnitTestContourFilter.cxx

@@ -176,24 +176,12 @@ public:
   void TestUnsupportedFlyingEdges() const
   {
     vtkm::cont::testing::MakeTestDataSet maker;
-
-    vtkm::cont::DataSet rectilinearDataset = maker.Make3DRectilinearDataSet0();
     vtkm::cont::DataSet explicitDataSet = maker.Make3DExplicitDataSet0();
 
     vtkm::filter::contour::ContourFlyingEdges filter;
     filter.SetIsoValue(2.0);
     filter.SetActiveField("pointvar");
 
-    try
-    {
-      filter.Execute(rectilinearDataset);
-      VTKM_TEST_FAIL("Flying Edges filter should not run on datasets with rectilinear coordinates");
-    }
-    catch (vtkm::cont::ErrorFilterExecution&)
-    {
-      std::cout << "Execution successfully aborted" << std::endl;
-    }
-
     try
     {
       filter.Execute(explicitDataSet);
@@ -205,6 +193,58 @@ public:
     }
   }
 
+  template <typename ContourFilterType>
+  void TestNonUniformStructured() const
+  {
+    auto pathname =
+      vtkm::cont::testing::Testing::DataPath("rectilinear/simple_rectilinear1_ascii.vtk");
+    vtkm::io::VTKDataSetReader reader(pathname);
+    vtkm::cont::DataSet rectilinearDataset = reader.ReadDataSet();
+
+    // Single-cell contour
+    ContourFilterType filter;
+    filter.SetActiveField("var");
+    filter.SetIsoValue(2.0);
+    vtkm::cont::DataSet outputSingleCell = filter.Execute(rectilinearDataset);
+    auto coordinates = outputSingleCell.GetCoordinateSystem()
+                         .GetData()
+                         .AsArrayHandle<vtkm::cont::ArrayHandle<vtkm::Vec3f>>();
+
+    VTKM_TEST_ASSERT(outputSingleCell.GetNumberOfPoints() == 3,
+                     "Wrong number of points in rectilinear contour");
+    VTKM_TEST_ASSERT(outputSingleCell.GetNumberOfCells() == 1,
+                     "Wrong number of cells in rectilinear contour");
+    VTKM_TEST_ASSERT(outputSingleCell.GetCellSet().GetCellShape(0) == vtkm::CELL_SHAPE_TRIANGLE,
+                     "Wrong contour cell shape");
+
+    auto expectedCoordinates =
+      vtkm::cont::make_ArrayHandle<vtkm::Vec3f>({ vtkm::Vec3f{ 10.0f, -10.0f, 9.66341f },
+                                                  vtkm::Vec3f{ 9.30578f, -10.0f, 10.0f },
+                                                  vtkm::Vec3f{ 10.0f, -9.78842f, 10.0f } });
+    VTKM_TEST_ASSERT(test_equal_ArrayHandles(coordinates, expectedCoordinates),
+                     "Wrong contour coordinates");
+
+    // Generating normals triggers a different worklet for Flying Edges pass 4,
+    // But it should not change anything on the contour itself.
+    filter.SetGenerateNormals(true);
+    vtkm::cont::DataSet outputNormals = filter.Execute(rectilinearDataset);
+    coordinates = outputSingleCell.GetCoordinateSystem()
+                    .GetData()
+                    .AsArrayHandle<vtkm::cont::ArrayHandle<vtkm::Vec3f>>();
+    VTKM_TEST_ASSERT(test_equal_ArrayHandles(coordinates, expectedCoordinates),
+                     "Wrong contour coordinates");
+
+    // Full contour
+    filter.SetIsoValue(3.0);
+    filter.SetGenerateNormals(false);
+    vtkm::cont::DataSet output = filter.Execute(rectilinearDataset);
+
+    VTKM_TEST_ASSERT(output.GetNumberOfPoints() == 93,
+                     "Wrong number of points in rectilinear contour");
+    VTKM_TEST_ASSERT(output.GetNumberOfCells() == 144,
+                     "Wrong number of cells in rectilinear contour");
+  }
+
   void operator()() const
   {
     this->TestContourUniformGrid<vtkm::filter::contour::Contour>(72);
@@ -220,6 +260,10 @@ public:
     this->TestContourWedges<vtkm::filter::contour::Contour>();
     this->TestContourWedges<vtkm::filter::contour::ContourMarchingCells>();
 
+    this->TestNonUniformStructured<vtkm::filter::contour::Contour>();
+    this->TestNonUniformStructured<vtkm::filter::contour::ContourFlyingEdges>();
+    this->TestNonUniformStructured<vtkm::filter::contour::ContourMarchingCells>();
+
     this->TestUnsupportedFlyingEdges();
   }
 

vtkm/filter/contour/worklet/ContourFlyingEdges.h

@@ -67,13 +67,14 @@ public:
 
   // Filter called without normals generation
   template <typename ValueType,
+            typename CoordsType,
             typename StorageTagField,
             typename CoordinateType,
             typename StorageTagVertices>
   vtkm::cont::CellSetSingleType<> Run(
     const std::vector<ValueType>& isovalues,
     const vtkm::cont::CellSetStructured<3>& cells,
-    const vtkm::cont::ArrayHandleUniformPointCoordinates& coordinateSystem,
+    const CoordsType& coordinateSystem,
     const vtkm::cont::ArrayHandle<ValueType, StorageTagField>& input,
     vtkm::cont::ArrayHandle<vtkm::Vec<CoordinateType, 3>, StorageTagVertices>& vertices)
   {
@@ -87,6 +88,7 @@ public:
 
   // Filter called with normals generation
   template <typename ValueType,
+            typename CoordsType,
             typename StorageTagField,
             typename CoordinateType,
             typename StorageTagVertices,
@@ -94,7 +96,7 @@ public:
   vtkm::cont::CellSetSingleType<> Run(
     const std::vector<ValueType>& isovalues,
     const vtkm::cont::CellSetStructured<3>& cells,
-    const vtkm::cont::ArrayHandleUniformPointCoordinates& coordinateSystem,
+    const CoordsType& coordinateSystem,
     const vtkm::cont::ArrayHandle<ValueType, StorageTagField>& input,
     vtkm::cont::ArrayHandle<vtkm::Vec<CoordinateType, 3>, StorageTagVertices>& vertices,
     vtkm::cont::ArrayHandle<vtkm::Vec<CoordinateType, 3>, StorageTagNormals>& normals)

vtkm/filter/contour/worklet/contour/FlyingEdges.h

@@ -41,6 +41,7 @@ vtkm::Id extend_by(vtkm::cont::ArrayHandle<T, S>& handle, vtkm::Id size)
 
 //----------------------------------------------------------------------------
 template <typename ValueType,
+          typename CoordsType,
           typename StorageTagField,
           typename StorageTagVertices,
           typename StorageTagNormals,
@@ -48,7 +49,7 @@ template <typename ValueType,
           typename NormalType>
 vtkm::cont::CellSetSingleType<> execute(
   const vtkm::cont::CellSetStructured<3>& cells,
-  const vtkm::cont::ArrayHandleUniformPointCoordinates& coordinateSystem,
+  const CoordsType coordinateSystem,
   const std::vector<ValueType>& isovalues,
   const vtkm::cont::ArrayHandle<ValueType, StorageTagField>& inputField,
   vtkm::cont::ArrayHandle<vtkm::Vec<CoordinateType, 3>, StorageTagVertices>& points,
@@ -56,18 +57,10 @@ vtkm::cont::CellSetSingleType<> execute(
   vtkm::worklet::contour::CommonState& sharedState)
 {
   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());
+  edgeCases.Allocate(coordinateSystem.GetData().GetNumberOfValues());
 
   vtkm::cont::CellSetStructured<2> metaDataMesh2D;
   vtkm::cont::ArrayHandle<vtkm::Id> metaDataLinearSums; //per point of metaDataMesh
@@ -158,8 +151,7 @@ vtkm::cont::CellSetSingleType<> execute(
       {
         VTKM_LOG_SCOPE(vtkm::cont::LogLevel::Perf, "FlyingEdges Pass4");
 
-        launchComputePass4 pass4(
+        auto pass4 = launchComputePass4(pdims, multiContourCellOffset, multiContourPointOffset);
-          pdims, origin, spacing, multiContourCellOffset, multiContourPointOffset);
 
         detail::extend_by(points, newPointSize);
         if (sharedState.GenerateNormals)
@@ -171,6 +163,7 @@ vtkm::cont::CellSetSingleType<> execute(
                                        pass4,
                                        newPointSize,
                                        isoval,
+                                       coordinateSystem,
                                        inputField,
                                        edgeCases,
                                        metaDataMesh2D,

vtkm/filter/contour/worklet/contour/FlyingEdgesPass4.h

@@ -28,20 +28,14 @@ namespace flying_edges
 struct launchComputePass4
 {
   vtkm::Id3 PointDims;
-  vtkm::Vec3f Origin;
-  vtkm::Vec3f Spacing;
 
   vtkm::Id CellWriteOffset;
   vtkm::Id PointWriteOffset;
 
   launchComputePass4(const vtkm::Id3& pdims,
-                     const vtkm::Vec3f& origin,
-                     const vtkm::Vec3f& spacing,
                      vtkm::Id multiContourCellOffset,
                      vtkm::Id multiContourPointOffset)
     : PointDims(pdims)
-    , Origin(origin)
-    , Spacing(spacing)
     , CellWriteOffset(multiContourCellOffset)
     , PointWriteOffset(multiContourPointOffset)
   {
@@ -49,6 +43,7 @@ struct launchComputePass4
 
   template <typename DeviceAdapterTag,
             typename T,
+            typename CoordsType,
             typename StorageTagField,
             typename MeshSums,
             typename PointType,
@@ -56,6 +51,7 @@ struct launchComputePass4
   VTKM_CONT bool LaunchXAxis(DeviceAdapterTag device,
                              vtkm::Id vtkmNotUsed(newPointSize),
                              T isoval,
+                             CoordsType coordinateSystem,
                              const vtkm::cont::ArrayHandle<T, StorageTagField>& inputField,
                              vtkm::cont::ArrayHandle<vtkm::UInt8> edgeCases,
                              vtkm::cont::CellSetStructured<2>& metaDataMesh2D,
@@ -71,12 +67,8 @@ struct launchComputePass4
     vtkm::cont::Invoker invoke(device);
     if (sharedState.GenerateNormals)
     {
-      ComputePass4XWithNormals<T> worklet4(isoval,
+      ComputePass4XWithNormals<T> worklet4(
-                                           this->PointDims,
+        isoval, this->PointDims, this->CellWriteOffset, this->PointWriteOffset);
-                                           this->Origin,
-                                           this->Spacing,
-                                           this->CellWriteOffset,
-                                           this->PointWriteOffset);
       invoke(worklet4,
              metaDataMesh2D,
              metaDataSums,
@@ -84,6 +76,7 @@ struct launchComputePass4
              metaDataMax,
              metaDataNumTris,
              edgeCases,
+             coordinateSystem,
              inputField,
              triangle_topology,
              sharedState.InterpolationEdgeIds,
@@ -94,12 +87,8 @@ struct launchComputePass4
     }
     else
     {
-      ComputePass4X<T> worklet4(isoval,
+      ComputePass4X<T> worklet4(
-                                this->PointDims,
+        isoval, this->PointDims, this->CellWriteOffset, this->PointWriteOffset);
-                                this->Origin,
-                                this->Spacing,
-                                this->CellWriteOffset,
-                                this->PointWriteOffset);
       invoke(worklet4,
              metaDataMesh2D,
              metaDataSums,
@@ -107,6 +96,7 @@ struct launchComputePass4
              metaDataMax,
              metaDataNumTris,
              edgeCases,
+             coordinateSystem,
              inputField,
              triangle_topology,
              sharedState.InterpolationEdgeIds,
@@ -120,6 +110,7 @@ struct launchComputePass4
 
   template <typename DeviceAdapterTag,
             typename T,
+            typename CoordsType,
             typename StorageTagField,
             typename MeshSums,
             typename PointType,
@@ -127,6 +118,7 @@ struct launchComputePass4
   VTKM_CONT bool LaunchYAxis(DeviceAdapterTag device,
                              vtkm::Id newPointSize,
                              T isoval,
+                             CoordsType coordinateSystem,
                              const vtkm::cont::ArrayHandle<T, StorageTagField>& inputField,
                              vtkm::cont::ArrayHandle<vtkm::UInt8> edgeCases,
                              vtkm::cont::CellSetStructured<2>& metaDataMesh2D,
@@ -157,11 +149,8 @@ struct launchComputePass4
            sharedState.CellIdMap);
 
     //This needs to be done on array handle view ( start = this->PointWriteOffset, len = newPointSize)
-    ComputePass5Y<T> worklet5(this->PointDims,
+    ComputePass5Y<T> worklet5(this->PointDims, this->PointWriteOffset, sharedState.GenerateNormals);
-                              this->Origin,
+
-                              this->Spacing,
-                              this->PointWriteOffset,
-                              sharedState.GenerateNormals);
     invoke(worklet5,
            vtkm::cont::make_ArrayHandleView(
              sharedState.InterpolationEdgeIds, this->PointWriteOffset, newPointSize),
@@ -169,6 +158,7 @@ struct launchComputePass4
              sharedState.InterpolationWeights, this->PointWriteOffset, newPointSize),
            vtkm::cont::make_ArrayHandleView(points, this->PointWriteOffset, newPointSize),
            inputField,
+           coordinateSystem,
            normals);
 
     return true;

vtkm/filter/contour/worklet/contour/FlyingEdgesPass4X.h

@@ -32,9 +32,6 @@ struct ComputePass4X : public vtkm::worklet::WorkletVisitCellsWithPoints
 {
 
   vtkm::Id3 PointDims;
-  vtkm::Vec3f Origin;
-  vtkm::Vec3f Spacing;
-
   T IsoValue;
 
   vtkm::Id CellWriteOffset;
@@ -43,13 +40,9 @@ struct ComputePass4X : public vtkm::worklet::WorkletVisitCellsWithPoints
   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)
@@ -62,6 +55,7 @@ struct ComputePass4X : public vtkm::worklet::WorkletVisitCellsWithPoints
                                 FieldInPoint axis_maxs,
                                 WholeArrayIn cell_tri_count,
                                 WholeArrayIn edgeData,
+                                WholeArrayIn coords,
                                 WholeArrayIn data,
                                 WholeArrayOut connectivity,
                                 WholeArrayOut edgeIds,
@@ -69,13 +63,14 @@ struct ComputePass4X : public vtkm::worklet::WorkletVisitCellsWithPoints
                                 WholeArrayOut inputCellIds,
                                 WholeArrayOut points);
   using ExecutionSignature =
-    void(ThreadIndices, _2, _3, _4, _5, _6, _7, _8, _9, _10, _11, _12, WorkIndex);
+    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 WholeCoordsField,
             typename WholeDataField,
             typename WholeConnField,
             typename WholeEdgeIdField,
@@ -88,6 +83,7 @@ struct ComputePass4X : public vtkm::worklet::WorkletVisitCellsWithPoints
                             const FieldInPointId& axis_maxs,
                             const WholeTriField& cellTriCount,
                             const WholeEdgeField& edges,
+                            const WholeCoordsField& coords,
                             const WholeDataField& field,
                             const WholeConnField& conn,
                             const WholeEdgeIdField& interpolatedEdgeIds,
@@ -141,6 +137,7 @@ struct ComputePass4X : public vtkm::worklet::WorkletVisitCellsWithPoints
                          field,
                          interpolatedEdgeIds,
                          weights,
+                         coords,
                          points,
                          state.startPos,
                          increments,
@@ -158,12 +155,14 @@ struct ComputePass4X : public vtkm::worklet::WorkletVisitCellsWithPoints
   template <typename WholeDataField,
             typename WholeIEdgeField,
             typename WholeWeightField,
+            typename WholeCoordsField,
             typename WholePointField>
   VTKM_EXEC inline void Generate(const vtkm::Vec<vtkm::UInt8, 3>& boundaryStatus,
                                  const vtkm::Id3& ijk,
                                  const WholeDataField& field,
                                  const WholeIEdgeField& interpolatedEdgeIds,
                                  const WholeWeightField& weights,
+                                 const WholeCoordsField coords,
                                  const WholePointField& points,
                                  const vtkm::Id4& startPos,
                                  const vtkm::Id3& incs,
@@ -188,7 +187,7 @@ struct ComputePass4X : public vtkm::worklet::WorkletVisitCellsWithPoints
         interpolatedEdgeIds.Set(writeIndex, pos);
         weights.Set(writeIndex, static_cast<vtkm::FloatDefault>(t));
 
-        auto coord = this->InterpolateCoordinate(t, ijk, ijk + vtkm::Id3{ 1, 0, 0 });
+        auto coord = this->InterpolateCoordinate(coords, t, ijk, ijk + vtkm::Id3{ 1, 0, 0 });
         points.Set(writeIndex, coord);
       }
       if (edgeUses[4])
@@ -201,7 +200,7 @@ struct ComputePass4X : public vtkm::worklet::WorkletVisitCellsWithPoints
         interpolatedEdgeIds.Set(writeIndex, pos);
         weights.Set(writeIndex, static_cast<vtkm::FloatDefault>(t));
 
-        auto coord = this->InterpolateCoordinate(t, ijk, ijk + vtkm::Id3{ 0, 1, 0 });
+        auto coord = this->InterpolateCoordinate(coords, t, ijk, ijk + vtkm::Id3{ 0, 1, 0 });
         points.Set(writeIndex, coord);
       }
       if (edgeUses[8])
@@ -214,7 +213,7 @@ struct ComputePass4X : public vtkm::worklet::WorkletVisitCellsWithPoints
         interpolatedEdgeIds.Set(writeIndex, pos);
         weights.Set(writeIndex, static_cast<vtkm::FloatDefault>(t));
 
-        auto coord = this->InterpolateCoordinate(t, ijk, ijk + vtkm::Id3{ 0, 0, 1 });
+        auto coord = this->InterpolateCoordinate(coords, t, ijk, ijk + vtkm::Id3{ 0, 0, 1 });
         points.Set(writeIndex, coord);
       }
     }
@@ -231,30 +230,30 @@ struct ComputePass4X : public vtkm::worklet::WorkletVisitCellsWithPoints
     const bool onZ = boundaryStatus[AxisToSum::zindex] & FlyingEdges3D::MaxBoundary;
     if (onX) //+x boundary
     {
-      this->InterpolateEdge(ijk, pos[0], incs, 5, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
+      this->InterpolateEdge(ijk, pos[0], incs, 5, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, coords, points);
-      this->InterpolateEdge(ijk, pos[0], incs, 9, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
+      this->InterpolateEdge(ijk, pos[0], incs, 9, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, coords, points);
       if (onY) //+x +y
       {
-        this->InterpolateEdge(ijk, pos[0], incs, 11, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
+        this->InterpolateEdge(ijk, pos[0], incs, 11, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, coords, points);
       }
       if (onZ) //+x +z
       {
-        this->InterpolateEdge(ijk, pos[0], incs, 7, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
+        this->InterpolateEdge(ijk, pos[0], incs, 7, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, coords, points);
       }
     }
     if (onY) //+y boundary
     {
-      this->InterpolateEdge(ijk, pos[0], incs, 1, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
+      this->InterpolateEdge(ijk, pos[0], incs, 1, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, coords, points);
-      this->InterpolateEdge(ijk, pos[0], incs, 10, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
+      this->InterpolateEdge(ijk, pos[0], incs, 10, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, coords, points);
       if (onZ) //+y +z boundary
       {
-        this->InterpolateEdge(ijk, pos[0], incs, 3, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
+        this->InterpolateEdge(ijk, pos[0], incs, 3, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, coords, points);
       }
     }
     if (onZ) //+z boundary
     {
-      this->InterpolateEdge(ijk, pos[0], incs, 2, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
+      this->InterpolateEdge(ijk, pos[0], incs, 2, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, coords, points);
-      this->InterpolateEdge(ijk, pos[0], incs, 6, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
+      this->InterpolateEdge(ijk, pos[0], incs, 6, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, coords, points);
     }
     // clang-format on
   }
@@ -264,6 +263,7 @@ struct ComputePass4X : public vtkm::worklet::WorkletVisitCellsWithPoints
   template <typename WholeField,
             typename WholeIEdgeField,
             typename WholeWeightField,
+            typename WholeCoordsField,
             typename WholePointField>
   VTKM_EXEC inline void InterpolateEdge(const vtkm::Id3& ijk,
                                         vtkm::Id currentIdx,
@@ -274,6 +274,7 @@ struct ComputePass4X : public vtkm::worklet::WorkletVisitCellsWithPoints
                                         const WholeField& field,
                                         const WholeIEdgeField& interpolatedEdgeIds,
                                         const WholeWeightField& weights,
+                                        const WholeCoordsField& coords,
                                         const WholePointField& points) const
   {
     using AxisToSum = SumXAxis;
@@ -300,30 +301,49 @@ struct ComputePass4X : public vtkm::worklet::WorkletVisitCellsWithPoints
     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);
+    auto coord = this->InterpolateCoordinate(coords, t, ijk + offsets1, ijk + offsets2);
     points.Set(writeIndex, coord);
   }
 
+  // Fast interpolation method for uniform coordinates
   //----------------------------------------------------------------------------
-  inline VTKM_EXEC vtkm::Vec3f InterpolateCoordinate(T t,
+  inline VTKM_EXEC vtkm::Vec3f InterpolateCoordinate(
-                                                     const vtkm::Id3& ijk0,
+    vtkm::internal::ArrayPortalUniformPointCoordinates coords,
-                                                     const vtkm::Id3& ijk1) const
+    T t,
+    const vtkm::Id3& ijk0,
+    const vtkm::Id3& ijk1) const
   {
     return vtkm::Vec3f(
-      this->Origin[0] +
+      coords.GetOrigin()[0] +
-        this->Spacing[0] *
+        coords.GetSpacing()[0] *
           (static_cast<vtkm::FloatDefault>(ijk0[0]) +
            static_cast<vtkm::FloatDefault>(t) * static_cast<vtkm::FloatDefault>(ijk1[0] - ijk0[0])),
-      this->Origin[1] +
+      coords.GetOrigin()[1] +
-        this->Spacing[1] *
+        coords.GetSpacing()[1] *
           (static_cast<vtkm::FloatDefault>(ijk0[1]) +
            static_cast<vtkm::FloatDefault>(t) * static_cast<vtkm::FloatDefault>(ijk1[1] - ijk0[1])),
-      this->Origin[2] +
+      coords.GetOrigin()[2] +
-        this->Spacing[2] *
+        coords.GetSpacing()[2] *
           (static_cast<vtkm::FloatDefault>(ijk0[2]) +
            static_cast<vtkm::FloatDefault>(t) *
              static_cast<vtkm::FloatDefault>(ijk1[2] - ijk0[2])));
   }
+
+  // Interpolation for explicit coordinates
+  //----------------------------------------------------------------------------
+  template <typename CoordsPortal>
+  inline VTKM_EXEC vtkm::Vec3f InterpolateCoordinate(CoordsPortal coords,
+                                                     T t,
+                                                     const vtkm::Id3& ijk0,
+                                                     const vtkm::Id3& ijk1) const
+  {
+    return (1.0f - static_cast<vtkm::FloatDefault>(t)) *
+      coords.Get(ijk0[0] + this->PointDims[0] * ijk0[1] +
+                 this->PointDims[0] * this->PointDims[1] * ijk0[2]) +
+      static_cast<vtkm::FloatDefault>(t) *
+      coords.Get(ijk1[0] + this->PointDims[0] * ijk1[1] +
+                 this->PointDims[0] * this->PointDims[1] * ijk1[2]);
+  }
 };
 }
 }

vtkm/filter/contour/worklet/contour/FlyingEdgesPass4XWithNormals.h

@@ -31,9 +31,6 @@ struct ComputePass4XWithNormals : public vtkm::worklet::WorkletVisitCellsWithPoi
 {
 
   vtkm::Id3 PointDims;
-  vtkm::Vec3f Origin;
-  vtkm::Vec3f Spacing;
-
   T IsoValue;
 
   vtkm::Id CellWriteOffset;
@@ -42,13 +39,9 @@ struct ComputePass4XWithNormals : public vtkm::worklet::WorkletVisitCellsWithPoi
   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)
@@ -61,6 +54,7 @@ struct ComputePass4XWithNormals : public vtkm::worklet::WorkletVisitCellsWithPoi
                                 FieldInPoint axis_maxs,
                                 WholeArrayIn cell_tri_count,
                                 WholeArrayIn edgeData,
+                                WholeArrayIn coords,
                                 WholeArrayIn data,
                                 WholeArrayOut connectivity,
                                 WholeArrayOut edgeIds,
@@ -69,13 +63,14 @@ struct ComputePass4XWithNormals : public vtkm::worklet::WorkletVisitCellsWithPoi
                                 WholeArrayOut points,
                                 WholeArrayOut normals);
   using ExecutionSignature =
-    void(ThreadIndices, _2, _3, _4, _5, _6, _7, _8, _9, _10, _11, _12, _13, WorkIndex);
+    void(ThreadIndices, _2, _3, _4, _5, _6, _7, _8, _9, _10, _11, _12, _13, _14, WorkIndex);
 
   template <typename ThreadIndices,
             typename FieldInPointId3,
             typename FieldInPointId,
             typename WholeTriField,
             typename WholeEdgeField,
+            typename WholeCoordsField,
             typename WholeDataField,
             typename WholeConnField,
             typename WholeEdgeIdField,
@@ -89,6 +84,7 @@ struct ComputePass4XWithNormals : public vtkm::worklet::WorkletVisitCellsWithPoi
                             const FieldInPointId& axis_maxs,
                             const WholeTriField& cellTriCount,
                             const WholeEdgeField& edges,
+                            const WholeCoordsField& coords,
                             const WholeDataField& field,
                             const WholeConnField& conn,
                             const WholeEdgeIdField& interpolatedEdgeIds,
@@ -144,6 +140,7 @@ struct ComputePass4XWithNormals : public vtkm::worklet::WorkletVisitCellsWithPoi
                          field,
                          interpolatedEdgeIds,
                          weights,
+                         coords,
                          points,
                          normals,
                          state.startPos,
@@ -162,6 +159,7 @@ struct ComputePass4XWithNormals : public vtkm::worklet::WorkletVisitCellsWithPoi
   template <typename WholeDataField,
             typename WholeIEdgeField,
             typename WholeWeightField,
+            typename WholeCoordsField,
             typename WholePointField,
             typename WholeNormalField>
   VTKM_EXEC inline void Generate(const vtkm::Vec<vtkm::UInt8, 3>& boundaryStatus,
@@ -169,6 +167,7 @@ struct ComputePass4XWithNormals : public vtkm::worklet::WorkletVisitCellsWithPoi
                                  const WholeDataField& field,
                                  const WholeIEdgeField& interpolatedEdgeIds,
                                  const WholeWeightField& weights,
+                                 const WholeCoordsField coords,
                                  const WholePointField& points,
                                  const WholeNormalField& normals,
                                  const vtkm::Id4& startPos,
@@ -197,7 +196,7 @@ struct ComputePass4XWithNormals : public vtkm::worklet::WorkletVisitCellsWithPoi
         weights.Set(writeIndex, static_cast<vtkm::FloatDefault>(t));
 
         auto ijk1 = ijk + vtkm::Id3{ 1, 0, 0 };
-        auto coord = this->InterpolateCoordinate(t, ijk, ijk1);
+        auto coord = this->InterpolateCoordinate(coords, t, ijk, ijk1);
         points.Set(writeIndex, coord);
 
         //gradient generation
@@ -216,7 +215,7 @@ struct ComputePass4XWithNormals : public vtkm::worklet::WorkletVisitCellsWithPoi
         weights.Set(writeIndex, static_cast<vtkm::FloatDefault>(t));
 
         auto ijk1 = ijk + vtkm::Id3{ 0, 1, 0 };
-        auto coord = this->InterpolateCoordinate(t, ijk, ijk1);
+        auto coord = this->InterpolateCoordinate(coords, t, ijk, ijk1);
         points.Set(writeIndex, coord);
 
         //gradient generation
@@ -235,7 +234,7 @@ struct ComputePass4XWithNormals : public vtkm::worklet::WorkletVisitCellsWithPoi
         weights.Set(writeIndex, static_cast<vtkm::FloatDefault>(t));
 
         auto ijk1 = ijk + vtkm::Id3{ 0, 0, 1 };
-        auto coord = this->InterpolateCoordinate(t, ijk, ijk1);
+        auto coord = this->InterpolateCoordinate(coords, t, ijk, ijk1);
         points.Set(writeIndex, coord);
 
         //gradient generation
@@ -258,38 +257,38 @@ struct ComputePass4XWithNormals : public vtkm::worklet::WorkletVisitCellsWithPoi
     if (onX) //+x boundary
     {
       this->InterpolateEdge(
-          fullyInterior, ijk, pos[0], incs, 5, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points, normals);
+          fullyInterior, ijk, pos[0], incs, 5, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, coords, points, normals);
       this->InterpolateEdge(
-          fullyInterior, ijk, pos[0], incs, 9, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points, normals);
+          fullyInterior, ijk, pos[0], incs, 9, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, coords, points, normals);
       if (onY) //+x +y
       {
         this->InterpolateEdge(
-          fullyInterior, ijk, pos[0], incs, 11, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points, normals);
+          fullyInterior, ijk, pos[0], incs, 11, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, coords, points, normals);
       }
       if (onZ) //+x +z
       {
         this->InterpolateEdge(
-          fullyInterior, ijk, pos[0], incs, 7, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points, normals);
+          fullyInterior, ijk, pos[0], incs, 7, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, coords, points, normals);
       }
     }
     if (onY) //+y boundary
     {
       this->InterpolateEdge(
-        fullyInterior, ijk, pos[0], incs, 1, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points, normals);
+        fullyInterior, ijk, pos[0], incs, 1, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, coords, points, normals);
       this->InterpolateEdge(
-        fullyInterior, ijk, pos[0], incs, 10, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points, normals);
+        fullyInterior, ijk, pos[0], incs, 10, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, coords, points, normals);
       if (onZ) //+y +z boundary
       {
         this->InterpolateEdge(
-          fullyInterior, ijk, pos[0], incs, 3, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points, normals);
+          fullyInterior, ijk, pos[0], incs, 3, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, coords, points, normals);
       }
     }
     if (onZ) //+z boundary
     {
       this->InterpolateEdge(
-        fullyInterior, ijk, pos[0], incs, 2, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points, normals);
+        fullyInterior, ijk, pos[0], incs, 2, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, coords, points, normals);
       this->InterpolateEdge(
-        fullyInterior, ijk, pos[0], incs, 6, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points, normals);
+        fullyInterior, ijk, pos[0], incs, 6, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, coords, points, normals);
     }
     // clang-format on
   }
@@ -300,6 +299,7 @@ struct ComputePass4XWithNormals : public vtkm::worklet::WorkletVisitCellsWithPoi
             typename WholeIEdgeField,
             typename WholeWeightField,
             typename WholePointField,
+            typename WholeCoordsField,
             typename WholeNormalField>
   VTKM_EXEC inline void InterpolateEdge(bool fullyInterior,
                                         const vtkm::Id3& ijk,
@@ -311,6 +311,7 @@ struct ComputePass4XWithNormals : public vtkm::worklet::WorkletVisitCellsWithPoi
                                         const WholeField& field,
                                         const WholeIEdgeField& interpolatedEdgeIds,
                                         const WholeWeightField& weights,
+                                        const WholeCoordsField& coords,
                                         const WholePointField& points,
                                         const WholeNormalField& normals) const
   {
@@ -338,36 +339,55 @@ struct ComputePass4XWithNormals : public vtkm::worklet::WorkletVisitCellsWithPoi
     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);
+    auto coord = this->InterpolateCoordinate(coords, t, ijk + offsets1, ijk + offsets2);
     points.Set(writeIndex, coord);
 
-    auto g0 = this->ComputeGradient(fullyInterior, ijk + offsets1, incs, iEdge[0], field);
     auto g1 = this->ComputeGradient(fullyInterior, ijk + offsets2, incs, iEdge[1], field);
+    auto g0 = this->ComputeGradient(fullyInterior, ijk + offsets1, incs, iEdge[0], field);
     g1 = g0 + (t * (g1 - g0));
     normals.Set(writeIndex, vtkm::Normal(g1));
   }
 
+  // Fast interpolation method for uniform coordinates
   //----------------------------------------------------------------------------
-  inline VTKM_EXEC vtkm::Vec3f InterpolateCoordinate(T t,
+  inline VTKM_EXEC vtkm::Vec3f InterpolateCoordinate(
-                                                     const vtkm::Id3& ijk0,
+    const vtkm::internal::ArrayPortalUniformPointCoordinates& coords,
-                                                     const vtkm::Id3& ijk1) const
+    T t,
+    const vtkm::Id3& ijk0,
+    const vtkm::Id3& ijk1) const
   {
     return vtkm::Vec3f(
-      this->Origin[0] +
+      coords.GetOrigin()[0] +
-        this->Spacing[0] *
+        coords.GetSpacing()[0] *
           (static_cast<vtkm::FloatDefault>(ijk0[0]) +
            static_cast<vtkm::FloatDefault>(t) * static_cast<vtkm::FloatDefault>(ijk1[0] - ijk0[0])),
-      this->Origin[1] +
+      coords.GetOrigin()[1] +
-        this->Spacing[1] *
+        coords.GetSpacing()[1] *
           (static_cast<vtkm::FloatDefault>(ijk0[1]) +
            static_cast<vtkm::FloatDefault>(t) * static_cast<vtkm::FloatDefault>(ijk1[1] - ijk0[1])),
-      this->Origin[2] +
+      coords.GetOrigin()[2] +
-        this->Spacing[2] *
+        coords.GetSpacing()[2] *
           (static_cast<vtkm::FloatDefault>(ijk0[2]) +
            static_cast<vtkm::FloatDefault>(t) *
              static_cast<vtkm::FloatDefault>(ijk1[2] - ijk0[2])));
   }
 
+  // Interpolation for explicit coordinates
+  //----------------------------------------------------------------------------
+  template <typename CoordsPortal>
+  inline VTKM_EXEC vtkm::Vec3f InterpolateCoordinate(const CoordsPortal& coords,
+                                                     T t,
+                                                     const vtkm::Id3& ijk0,
+                                                     const vtkm::Id3& ijk1) const
+  {
+    return (1.0f - static_cast<vtkm::FloatDefault>(t)) *
+      coords.Get(ijk0[0] + this->PointDims[0] * ijk0[1] +
+                 this->PointDims[0] * this->PointDims[1] * ijk0[2]) +
+      static_cast<vtkm::FloatDefault>(t) *
+      coords.Get(ijk1[0] + this->PointDims[0] * ijk1[1] +
+                 this->PointDims[0] * this->PointDims[1] * ijk1[2]);
+  }
+
   //----------------------------------------------------------------------------
   template <typename WholeDataField>
   VTKM_EXEC vtkm::Vec3f ComputeGradient(bool fullyInterior,

vtkm/filter/contour/worklet/contour/FlyingEdgesPass4Y.h

@@ -276,17 +276,11 @@ struct ComputePass4Y : public vtkm::worklet::WorkletVisitCellsWithPoints
 template <typename T>
 struct ComputePass5Y : public vtkm::worklet::WorkletMapField
 {
-
+  vtkm::Id3 PointDims;
-  vtkm::internal::ArrayPortalUniformPointCoordinates Coordinates;
   vtkm::Id NormalWriteOffset;
 
-  ComputePass5Y() {}
+  ComputePass5Y(const vtkm::Id3& pdims, vtkm::Id normalWriteOffset, bool generateNormals)
-  ComputePass5Y(const vtkm::Id3& pdims,
+    : PointDims(pdims)
-                const vtkm::Vec3f& origin,
-                const vtkm::Vec3f& spacing,
-                vtkm::Id normalWriteOffset,
-                bool generateNormals)
-    : Coordinates(pdims, origin, spacing)
     , NormalWriteOffset(normalWriteOffset)
   {
     if (!generateNormals)
@@ -299,20 +293,25 @@ struct ComputePass5Y : public vtkm::worklet::WorkletMapField
                                 FieldIn interpWeight,
                                 FieldOut points,
                                 WholeArrayIn field,
+                                WholeArrayIn coords,
                                 WholeArrayOut normals);
-  using ExecutionSignature = void(_1, _2, _3, _4, _5, WorkIndex);
+  using ExecutionSignature = void(_1, _2, _3, _4, _5, _6, WorkIndex);
 
-  template <typename PT, typename WholeInputField, typename WholeNormalField>
+  template <typename PT,
+            typename WholeInputField,
+            typename WholeNormalField,
+            typename WholeCoordsField>
   VTKM_EXEC void operator()(const vtkm::Id2& interpEdgeIds,
                             vtkm::FloatDefault weight,
                             vtkm::Vec<PT, 3>& outPoint,
                             const WholeInputField& field,
+                            const WholeCoordsField& coords,
                             WholeNormalField& normals,
                             vtkm::Id oidx) const
   {
     {
-      vtkm::Vec3f point1 = this->Coordinates.Get(interpEdgeIds[0]);
+      vtkm::Vec3f point1 = coords.Get(interpEdgeIds[0]);
-      vtkm::Vec3f point2 = this->Coordinates.Get(interpEdgeIds[1]);
+      vtkm::Vec3f point2 = coords.Get(interpEdgeIds[1]);
       outPoint = vtkm::Lerp(point1, point2, weight);
     }
 
@@ -320,15 +319,13 @@ struct ComputePass5Y : public vtkm::worklet::WorkletMapField
     if (this->NormalWriteOffset >= 0)
     {
       vtkm::Vec<T, 3> g0, g1;
-      const vtkm::Id3& dims = this->Coordinates.GetDimensions();
+      vtkm::Id3 ijk{ interpEdgeIds[0] % this->PointDims[0],
-      vtkm::Id3 ijk{ interpEdgeIds[0] % dims[0],
+                     (interpEdgeIds[0] / this->PointDims[0]) % this->PointDims[1],
-                     (interpEdgeIds[0] / dims[0]) % dims[1],
+                     interpEdgeIds[0] / (this->PointDims[0] * this->PointDims[1]) };
-                     interpEdgeIds[0] / (dims[0] * dims[1]) };
 
       vtkm::worklet::gradient::StructuredPointGradient gradient;
-      vtkm::exec::BoundaryState boundary(ijk, dims);
+      vtkm::exec::BoundaryState boundary(ijk, this->PointDims);
-      vtkm::exec::FieldNeighborhood<vtkm::internal::ArrayPortalUniformPointCoordinates>
+      vtkm::exec::FieldNeighborhood<WholeCoordsField> coord_neighborhood(coords, boundary);
-        coord_neighborhood(this->Coordinates, boundary);
 
       vtkm::exec::FieldNeighborhood<WholeInputField> field_neighborhood(field, boundary);
 
@@ -337,9 +334,9 @@ struct ComputePass5Y : public vtkm::worklet::WorkletMapField
       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],
+      boundary.IJK = vtkm::Id3{ interpEdgeIds[1] % this->PointDims[0],
-                                (interpEdgeIds[1] / dims[0]) % dims[1],
+                                (interpEdgeIds[1] / this->PointDims[0]) % this->PointDims[1],
-                                interpEdgeIds[1] / (dims[0] * dims[1]) };
+                                interpEdgeIds[1] / (this->PointDims[0] * this->PointDims[1]) };
       gradient(boundary, coord_neighborhood, field_neighborhood, g1);
 
       vtkm::Vec3f n = vtkm::Lerp(g0, g1, weight);