Add a form of CellInterpolate that operates on whole cell sets

The initial implementation of `CellInterpolate` takes arguments that are
expected from a topology map worklet. However, sometimes you want to
interplate cells that are queried from locators or otherwise come from a
`WholeCellSet` control signature argument.

A new form of `CellInterpolate` is added to handle this case.

docs/changelog/cell-interpolation.global.md

@@ -0,0 +1,8 @@
+## Add a form of CellInterpolate that operates on whole cell sets
+
+The initial implementation of `CellInterpolate` takes arguments that are
+expected from a topology map worklet. However, sometimes you want to
+interplate cells that are queried from locators or otherwise come from a
+`WholeCellSet` control signature argument.
+
+A new form of `CellInterpolate` is added to handle this case.

docs/users-guide/examples/GuideExampleCellOperations.cxx

@@ -11,9 +11,10 @@
 #include <vtkm/exec/CellInterpolate.h>
 #include <vtkm/exec/ParametricCoordinates.h>
 
-#include <vtkm/worklet/DispatcherMapTopology.h>
+#include <vtkm/worklet/WorkletMapField.h>
 #include <vtkm/worklet/WorkletMapTopology.h>
 
+#include <vtkm/cont/Invoker.h>
 #include <vtkm/cont/testing/MakeTestDataSet.h>
 #include <vtkm/cont/testing/Testing.h>
 
@@ -25,7 +26,7 @@ namespace
 ////
 struct CellCenters : vtkm::worklet::WorkletVisitCellsWithPoints
 {
-  using ControlSignature = void(CellSetIn,
+  using ControlSignature = void(CellSetIn inputCells,
                                 FieldInPoint inputField,
                                 FieldOutCell outputField);
   using ExecutionSignature = void(CellShape, PointCount, _2, _3);
@@ -52,6 +53,46 @@ struct CellCenters : vtkm::worklet::WorkletVisitCellsWithPoints
 //// END-EXAMPLE CellCenters
 ////
 
+////
+//// BEGIN-EXAMPLE CellLookupInterp
+////
+struct CellLookupInterp : vtkm::worklet::WorkletMapField
+{
+  using ControlSignature = void(WholeCellSetIn<> inputCells,
+                                WholeArrayIn inputField,
+                                FieldOut outputField);
+  using ExecutionSignature = void(InputIndex, _1, _2, _3);
+  using InputDomain = _3;
+
+  template<typename StructureType, typename FieldInPortalType, typename FieldOutType>
+  VTKM_EXEC void operator()(vtkm::Id index,
+                            const StructureType& structure,
+                            const FieldInPortalType& inputField,
+                            FieldOutType& outputField) const
+  {
+    // Normally you would use something like a locator to find the index to
+    // a cell that matches some query criteria. For demonstration purposes,
+    // we are just using a passed in index.
+    auto shape = structure.GetCellShape(index);
+    vtkm::IdComponent pointCount = structure.GetNumberOfIndices(index);
+
+    vtkm::Vec3f center;
+    vtkm::ErrorCode status =
+      vtkm::exec::ParametricCoordinatesCenter(pointCount, shape, center);
+    if (status != vtkm::ErrorCode::Success)
+    {
+      this->RaiseError(vtkm::ErrorString(status));
+      return;
+    }
+
+    auto pointIndices = structure.GetIndices(index);
+    vtkm::exec::CellInterpolate(pointIndices, inputField, center, shape, outputField);
+  }
+};
+////
+//// END-EXAMPLE CellLookupInterp
+////
+
 void TryCellCenters()
 {
   std::cout << "Trying CellCenters worklet." << std::endl;
@@ -62,14 +103,25 @@ void TryCellCenters()
   using ArrayType = vtkm::cont::ArrayHandle<vtkm::Float32>;
   ArrayType centers;
 
-  vtkm::worklet::DispatcherMapTopology<CellCenters> dispatcher;
-  dispatcher.Invoke(dataSet.GetCellSet(),
-                    dataSet.GetField("pointvar").GetData().AsArrayHandle<ArrayType>(),
-                    centers);
-
+  vtkm::cont::Invoker invoke;
+  invoke(CellCenters{},
+         dataSet.GetCellSet(),
+         dataSet.GetField("pointvar").GetData().AsArrayHandle<ArrayType>(),
+         centers);
   vtkm::cont::printSummary_ArrayHandle(centers, std::cout);
   std::cout << std::endl;
+  VTKM_TEST_ASSERT(centers.GetNumberOfValues() ==
+                     dataSet.GetCellSet().GetNumberOfCells(),
+                   "Bad number of cells.");
+  VTKM_TEST_ASSERT(test_equal(60.1875, centers.ReadPortal().Get(0)), "Bad first value.");
 
+  centers.Fill(0);
+  invoke(CellLookupInterp{},
+         dataSet.GetCellSet(),
+         dataSet.GetField("pointvar").GetData().AsArrayHandle<ArrayType>(),
+         centers);
+  vtkm::cont::printSummary_ArrayHandle(centers, std::cout);
+  std::cout << std::endl;
   VTKM_TEST_ASSERT(centers.GetNumberOfValues() ==
                      dataSet.GetCellSet().GetNumberOfCells(),
                    "Bad number of cells.");
@@ -122,11 +174,13 @@ void TryCellDerivatives()
   using ArrayType = vtkm::cont::ArrayHandle<vtkm::Float32>;
   vtkm::cont::ArrayHandle<vtkm::Vec3f_32> derivatives;
 
+  vtkm::cont::Invoker invoke;
   vtkm::worklet::DispatcherMapTopology<CellDerivatives> dispatcher;
-  dispatcher.Invoke(dataSet.GetCellSet(),
-                    dataSet.GetField("pointvar").GetData().AsArrayHandle<ArrayType>(),
-                    dataSet.GetCoordinateSystem().GetData(),
-                    derivatives);
+  invoke(CellDerivatives{},
+         dataSet.GetCellSet(),
+         dataSet.GetField("pointvar").GetData().AsArrayHandle<ArrayType>(),
+         dataSet.GetCoordinateSystem().GetData(),
+         derivatives);
 
   vtkm::cont::printSummary_ArrayHandle(derivatives, std::cout);
   std::cout << std::endl;

docs/users-guide/working-with-cells.rst

@@ -168,6 +168,16 @@ The :file:`vtkm/exec/CellInterpolate.h` header contains the function :func:`vtkm
    :file: GuideExampleCellOperations.cxx
    :caption: Interpolating field values to a cell's center.
 
+The previous form of :func:`vtkm::exec::CellInterpolate` is typically used within a :class:`vtkm::worklet::WorkletVisitCellsWithPoints` to interpolate within cells provided by the worklet.
+There is another form of :func:`vtkm::exec::CellInterpolate` that is used when provided with a cell set structure from a ``WholeCellSetIn``.
+Using ``WholeCellSetIn`` is described in more detail in :secref:`globals:Whole Cell Sets`, but the summary is that you can get the shape and incident points of any cell in the mesh.
+In this case, the alternate form of :func:`vtkm::exec::CellInterpolate` takes a ``Vec`` of the point indices for the cell and an array portal of all field values.
+
+.. doxygenfunction:: vtkm::exec::CellInterpolate(const IndicesVecType&, const FieldPortalType&, const vtkm::Vec<ParametricCoordType, 3>&, CellShapeTag, typename FieldPortalType::ValueType&)
+
+.. load-example:: CellLookupInterp
+   :file: GuideExampleCellOperations.cxx
+   :caption: Interpolating field values in a cell queried from a global structure.
 
 ------------------------------
 Derivatives

vtkm/exec/CellInterpolate.h

@@ -13,6 +13,8 @@
 #include <vtkm/CellShape.h>
 #include <vtkm/ErrorCode.h>
 #include <vtkm/VecAxisAlignedPointCoordinates.h>
+#include <vtkm/VecFromPortalPermute.h>
+#include <vtkm/cont/ArrayHandleUniformPointCoordinates.h>
 #include <vtkm/exec/FunctorBase.h>
 
 #include <lcl/lcl.h>
@@ -156,7 +158,7 @@ VTKM_EXEC vtkm::ErrorCode CellInterpolate(const vtkm::VecAxisAlignedPointCoordin
 /// \brief Interpolate a point field in a cell.
 ///
 /// Given the point field values for each node and the parametric coordinates
-/// of a point within the cell, interpolates the field to that point.
+/// of a location within the cell, interpolates the field to that location.
 ///
 /// @param[in]  pointFieldValues A list of field values for each point in the cell. This
 ///     usually comes from a `FieldInPoint` argument in a
@@ -184,6 +186,38 @@ VTKM_EXEC vtkm::ErrorCode CellInterpolate(const FieldVecType& pointFieldValues,
   return status;
 }
 
+//-----------------------------------------------------------------------------
+/// @brief Interpolate a point field in a cell.
+///
+/// Given the indices of the points for each node in a `Vec`, a portal to the point
+/// field values, and the parametric coordinates of a location within the cell, interpolates
+/// to that location.
+///
+/// @param[in]  pointIndices A list of point indices for each point in the cell. This
+///     usually comes from a `GetIndices()` call on the structure object provided by
+///     a `WholeCellSetIn` argument to a worklet.
+/// @param[in]  pointFieldPortal An array portal containing all the values in a point
+///     field array. This usually comes from a `WholeArrayIn` worklet argument.
+/// @param[in]  parametricCoords The parametric coordinates where you want to get the
+///     interpolaged field value for.
+/// @param[in]  shape A tag of type `CellShapeTag*` to identify the shape of the cell.
+/// @param[out] result Value to store the interpolated field.
+template <typename IndicesVecType,
+          typename FieldPortalType,
+          typename ParametricCoordType,
+          typename CellShapeTag>
+VTKM_EXEC vtkm::ErrorCode CellInterpolate(const IndicesVecType& pointIndices,
+                                          const FieldPortalType& pointFieldPortal,
+                                          const vtkm::Vec<ParametricCoordType, 3>& parametricCoords,
+                                          CellShapeTag shape,
+                                          typename FieldPortalType::ValueType& result)
+{
+  return CellInterpolate(vtkm::make_VecFromPortalPermute(&pointIndices, pointFieldPortal),
+                         parametricCoords,
+                         shape,
+                         result);
+}
+
 }
 } // namespace vtkm::exec
 

vtkm/exec/testing/UnitTestCellDerivative.cxx

@@ -17,7 +17,7 @@
 #include <vtkm/StaticAssert.h>
 #include <vtkm/VecVariable.h>
 
-#include <vtkm/testing/Testing.h>
+#include <vtkm/cont/testing/Testing.h>
 
 #include <ctime>
 #include <random>
@@ -330,5 +330,5 @@ void TestDerivative()
 
 int UnitTestCellDerivative(int argc, char* argv[])
 {
-  return vtkm::testing::Testing::Run(TestDerivative, argc, argv);
+  return vtkm::cont::testing::Testing::Run(TestDerivative, argc, argv);
 }

vtkm/exec/testing/UnitTestCellInterpolate.cxx

@@ -18,7 +18,7 @@
 #include <vtkm/VecAxisAlignedPointCoordinates.h>
 #include <vtkm/VecVariable.h>
 
-#include <vtkm/testing/Testing.h>
+#include <vtkm/cont/testing/Testing.h>
 
 #define CHECK_CALL(call) \
   VTKM_TEST_ASSERT((call) == vtkm::ErrorCode::Success, "Call resulted in error.")
@@ -91,6 +91,49 @@ struct TestInterpolateFunctor
                      "Interpolation at center not average value.");
   }
 
+  template <typename CellShapeTag, typename IndexVecType, typename FieldPortalType>
+  void DoTestWithIndices(CellShapeTag shape,
+                         const IndexVecType& pointIndices,
+                         const FieldPortalType& fieldValues) const
+  {
+    vtkm::IdComponent numPoints = pointIndices.GetNumberOfComponents();
+    if (numPoints < 1)
+    {
+      return;
+    }
+
+    FieldType averageValue = vtkm::TypeTraits<FieldType>::ZeroInitialization();
+    for (vtkm::IdComponent pointIndex = 0; pointIndex < numPoints; pointIndex++)
+    {
+      averageValue = averageValue + fieldValues.Get(pointIndices[pointIndex]);
+    }
+    averageValue = static_cast<ComponentType>(1.0 / numPoints) * averageValue;
+
+    for (vtkm::IdComponent pointIndex = 0; pointIndex < numPoints; pointIndex++)
+    {
+      vtkm::Vec3f pcoord;
+      CHECK_CALL(vtkm::exec::ParametricCoordinatesPoint(numPoints, pointIndex, shape, pcoord));
+      FieldType interpolatedValue;
+      CHECK_CALL(
+        vtkm::exec::CellInterpolate(pointIndices, fieldValues, pcoord, shape, interpolatedValue));
+
+      VTKM_TEST_ASSERT(test_equal(fieldValues.Get(pointIndices[pointIndex]), interpolatedValue),
+                       "Interpolation at point not point value.");
+    }
+
+    if (shape.Id != vtkm::CELL_SHAPE_POLY_LINE)
+    {
+      vtkm::Vec3f pcoord;
+      CHECK_CALL(vtkm::exec::ParametricCoordinatesCenter(numPoints, shape, pcoord));
+      FieldType interpolatedValue;
+      CHECK_CALL(
+        vtkm::exec::CellInterpolate(pointIndices, fieldValues, pcoord, shape, interpolatedValue));
+
+      VTKM_TEST_ASSERT(test_equal(averageValue, interpolatedValue),
+                       "Interpolation at center not average value.");
+    }
+  }
+
   template <typename CellShapeTag>
   void DoTest(CellShapeTag shape, vtkm::IdComponent numPoints) const
   {
@@ -102,6 +145,19 @@ struct TestInterpolateFunctor
     }
 
     this->DoTestWithField(shape, fieldValues);
+
+    vtkm::cont::ArrayHandle<FieldType> fieldArray;
+    fieldArray.Allocate(41);
+    SetPortal(fieldArray.WritePortal());
+
+    vtkm::VecVariable<vtkm::Id, MAX_POINTS> pointIndices;
+    for (vtkm::IdComponent pointIndex = 0; pointIndex < numPoints; pointIndex++)
+    {
+      vtkm::Id globalIndex = (7 + (13 * pointIndex)) % 41;
+      pointIndices.Append(globalIndex);
+    }
+
+    this->DoTestWithIndices(shape, pointIndices, fieldArray.ReadPortal());
   }
 
   template <typename CellShapeTag>
@@ -154,5 +210,5 @@ void TestInterpolate()
 
 int UnitTestCellInterpolate(int argc, char* argv[])
 {
-  return vtkm::testing::Testing::Run(TestInterpolate, argc, argv);
+  return vtkm::cont::testing::Testing::Run(TestInterpolate, argc, argv);
 }

vtkm/exec/testing/UnitTestParametricCoordinates.cxx

@@ -15,7 +15,7 @@
 #include <vtkm/StaticAssert.h>
 #include <vtkm/VecVariable.h>
 
-#include <vtkm/testing/Testing.h>
+#include <vtkm/cont/testing/Testing.h>
 
 #include <ctime>
 #include <random>
@@ -232,5 +232,5 @@ void TestAllPCoords()
 
 int UnitTestParametricCoordinates(int argc, char* argv[])
 {
-  return vtkm::testing::Testing::Run(TestAllPCoords, argc, argv);
+  return vtkm::cont::testing::Testing::Run(TestAllPCoords, argc, argv);
 }