Add CellLocatorChooser

`CellLocatorChooser` allows you to select a cell locator at compile
time. Unlike `CellLocatorGeneral`, you have to deduce what cell set
types you might have, but you don't have to have a large switch
statement in the middle of your worklet.

vtkm/cont/CMakeLists.txt

@@ -55,6 +55,7 @@ set(headers
   BoundsGlobalCompute.h
   CastAndCall.h
   CellLocatorBoundingIntervalHierarchy.h
+  CellLocatorChooser.h
   CellLocatorGeneral.h
   CellLocatorRectilinearGrid.h
   CellLocatorTwoLevel.h

vtkm/cont/CellLocator.h

@@ -29,9 +29,9 @@ namespace vtkm
 namespace cont
 {
 
-class VTKM_CONT_EXPORT VTKM_DEPRECATED(
-  1.6,
-  "CellLocator with virtual methods no longer supported. Use CellLocatorGeneral.") CellLocator
+class VTKM_CONT_EXPORT VTKM_DEPRECATED(1.6,
+                                       "CellLocator with virtual methods no longer supported. Use "
+                                       "CellLocatorGeneral or CellLocatorChooser.") CellLocator
   : public vtkm::cont::ExecutionObjectBase
 {
   VTKM_DEPRECATED_SUPPRESS_BEGIN

vtkm/cont/CellLocatorChooser.h

@@ -0,0 +1,65 @@
+//============================================================================
+//  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_cont_CellLocatorChooser_h
+#define vtk_m_cont_CellLocatorChooser_h
+
+#include <vtkm/cont/CellLocatorRectilinearGrid.h>
+#include <vtkm/cont/CellLocatorTwoLevel.h>
+#include <vtkm/cont/CellLocatorUniformGrid.h>
+#include <vtkm/cont/CellSetStructured.h>
+
+namespace vtkm
+{
+namespace cont
+{
+
+namespace detail
+{
+
+template <typename CellSetType, typename CoordinateSystemArrayType>
+struct CellLocatorChooserImpl
+{
+  using type = vtkm::cont::CellLocatorTwoLevel;
+};
+
+template <>
+struct CellLocatorChooserImpl<vtkm::cont::CellSetStructured<3>,
+                              vtkm::cont::ArrayHandleUniformPointCoordinates>
+{
+  using type = vtkm::cont::CellLocatorUniformGrid;
+};
+
+template <>
+struct CellLocatorChooserImpl<
+  vtkm::cont::CellSetStructured<3>,
+  vtkm::cont::ArrayHandleCartesianProduct<vtkm::cont::ArrayHandle<vtkm::FloatDefault>,
+                                          vtkm::cont::ArrayHandle<vtkm::FloatDefault>,
+                                          vtkm::cont::ArrayHandle<vtkm::FloatDefault>>>
+{
+  using type = vtkm::cont::CellLocatorRectilinearGrid;
+};
+
+} // namespace detail
+
+/// \brief A template to select an appropriate CellLocator based on CellSet type.
+///
+/// Given a concrete type for a `CellSet` subclass and a type of `ArrayHandle` for the
+/// coordinate system, `CellLocatorChooser` picks an appropriate `CellLocator` for that
+/// type of grid. It is a convenient class to use when you can resolve your templates
+/// to discover the type of data set being used for location.
+///
+template <typename CellSetType, typename CoordinateSystemArrayType>
+using CellLocatorChooser =
+  typename detail::CellLocatorChooserImpl<CellSetType, CoordinateSystemArrayType>::type;
+
+}
+} // namespace vtkm::cont
+
+#endif //vtk_m_cont_CellLocatorChooser_h

vtkm/cont/CellLocatorGeneral.h

@@ -26,6 +26,18 @@ namespace vtkm
 namespace cont
 {
 
+/// \brief A CellLocator that works generally well for any supported cell set.
+///
+/// `CellLocatorGeneral` creates a `CellLocator` that acts like a multiplexer to
+/// switch at runtime to any supported cell set. It is a convenient class to use
+/// when the type of `CellSet` cannot be determined at runtime.
+///
+/// Note that `CellLocatorGeneral` only supports a finite amount of `CellSet` types.
+/// Thus, it is possible to give it a cell set type that is not supported.
+///
+/// Also note that `CellLocatorGeneral` can add a significant amount of code inside
+/// of worklet that uses it, and this might cause some issues with some compilers.
+///
 class VTKM_CONT_EXPORT CellLocatorGeneral
   : public vtkm::cont::internal::CellLocatorBase<CellLocatorGeneral>
 {

vtkm/cont/testing/CMakeLists.txt

@@ -55,6 +55,7 @@ set(unit_tests
   UnitTestArrayHandleVirtual.cxx
   UnitTestArrayHandleXGCCoordinates.cxx
   UnitTestArrayPortalToIterators.cxx
+  UnitTestCellLocatorChooser.cxx
   UnitTestCellLocatorGeneral.cxx
   UnitTestCellSet.cxx
   UnitTestCellSetExplicit.cxx

vtkm/cont/testing/UnitTestCellLocatorChooser.cxx

@@ -0,0 +1,232 @@
+//============================================================================
+//  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.
+//============================================================================
+#include <vtkm/cont/CellLocatorChooser.h>
+
+#include <vtkm/cont/ArrayHandle.h>
+#include <vtkm/cont/DataSetBuilderRectilinear.h>
+#include <vtkm/cont/DataSetBuilderUniform.h>
+#include <vtkm/cont/testing/Testing.h>
+#include <vtkm/exec/CellInterpolate.h>
+#include <vtkm/worklet/DispatcherMapField.h>
+#include <vtkm/worklet/DispatcherMapTopology.h>
+#include <vtkm/worklet/ScatterPermutation.h>
+#include <vtkm/worklet/WorkletMapField.h>
+#include <vtkm/worklet/WorkletMapTopology.h>
+
+#include <ctime>
+#include <random>
+
+namespace
+{
+
+std::default_random_engine RandomGenerator;
+
+using PointType = vtkm::Vec3f;
+
+//-----------------------------------------------------------------------------
+vtkm::cont::DataSet MakeTestDataSetUniform()
+{
+  return vtkm::cont::DataSetBuilderUniform::Create(
+    vtkm::Id3{ 64 }, PointType{ -32.0f }, PointType{ 1.0f / 64.0f });
+}
+
+vtkm::cont::DataSet MakeTestDataSetRectilinear()
+{
+  std::uniform_real_distribution<vtkm::FloatDefault> coordGen(1.0f / 128.0f, 1.0f / 32.0f);
+
+  vtkm::cont::ArrayHandle<vtkm::FloatDefault> coords[3];
+  for (int i = 0; i < 3; ++i)
+  {
+    coords[i].Allocate(64);
+    auto portal = coords[i].WritePortal();
+
+    vtkm::FloatDefault cur = 0.0f;
+    for (vtkm::Id j = 0; j < portal.GetNumberOfValues(); ++j)
+    {
+      cur += coordGen(RandomGenerator);
+      portal.Set(j, cur);
+    }
+  }
+
+  return vtkm::cont::DataSetBuilderRectilinear::Create(coords[0], coords[1], coords[2]);
+}
+
+vtkm::cont::DataSet MakeTestDataSetCurvilinear()
+{
+  auto recti = MakeTestDataSetRectilinear();
+  auto coords = recti.GetCoordinateSystem().GetDataAsMultiplexer();
+
+  vtkm::cont::ArrayHandle<PointType> sheared;
+  sheared.Allocate(coords.GetNumberOfValues());
+
+  auto inPortal = coords.ReadPortal();
+  auto outPortal = sheared.WritePortal();
+  for (vtkm::Id i = 0; i < inPortal.GetNumberOfValues(); ++i)
+  {
+    auto val = inPortal.Get(i);
+    outPortal.Set(i, val + vtkm::make_Vec(val[1], val[2], val[0]));
+  }
+
+  vtkm::cont::DataSet curvi;
+  curvi.SetCellSet(recti.GetCellSet());
+  curvi.AddCoordinateSystem(vtkm::cont::CoordinateSystem("coords", sheared));
+
+  return curvi;
+}
+
+//-----------------------------------------------------------------------------
+class ParametricToWorldCoordinates : public vtkm::worklet::WorkletVisitCellsWithPoints
+{
+public:
+  using ControlSignature = void(CellSetIn cellset,
+                                FieldInPoint coords,
+                                FieldInOutCell pcs,
+                                FieldOutCell wcs);
+  using ExecutionSignature = void(CellShape, _2, _3, _4);
+
+  using ScatterType = vtkm::worklet::ScatterPermutation<>;
+
+  VTKM_CONT
+  static ScatterType MakeScatter(const vtkm::cont::ArrayHandle<vtkm::Id>& cellIds)
+  {
+    return ScatterType(cellIds);
+  }
+
+  template <typename CellShapeTagType, typename PointsVecType>
+  VTKM_EXEC void operator()(CellShapeTagType cellShape,
+                            PointsVecType points,
+                            const PointType& pc,
+                            PointType& wc) const
+  {
+    auto status = vtkm::exec::CellInterpolate(points, pc, cellShape, wc);
+    if (status != vtkm::ErrorCode::Success)
+    {
+      this->RaiseError(vtkm::ErrorString(status));
+    }
+  }
+};
+
+void GenerateRandomInput(const vtkm::cont::DataSet& ds,
+                         vtkm::Id count,
+                         vtkm::cont::ArrayHandle<vtkm::Id>& cellIds,
+                         vtkm::cont::ArrayHandle<PointType>& pcoords,
+                         vtkm::cont::ArrayHandle<PointType>& wcoords)
+{
+  vtkm::Id numberOfCells = ds.GetNumberOfCells();
+
+  std::uniform_int_distribution<vtkm::Id> cellIdGen(0, numberOfCells - 1);
+  std::uniform_real_distribution<vtkm::FloatDefault> pcoordGen(0.0f, 1.0f);
+
+  cellIds.Allocate(count);
+  pcoords.Allocate(count);
+  wcoords.Allocate(count);
+
+  for (vtkm::Id i = 0; i < count; ++i)
+  {
+    cellIds.WritePortal().Set(i, cellIdGen(RandomGenerator));
+
+    PointType pc{ pcoordGen(RandomGenerator),
+                  pcoordGen(RandomGenerator),
+                  pcoordGen(RandomGenerator) };
+    pcoords.WritePortal().Set(i, pc);
+  }
+
+  vtkm::worklet::DispatcherMapTopology<ParametricToWorldCoordinates> dispatcher(
+    ParametricToWorldCoordinates::MakeScatter(cellIds));
+  dispatcher.Invoke(
+    ds.GetCellSet(), ds.GetCoordinateSystem().GetDataAsMultiplexer(), pcoords, wcoords);
+}
+
+//-----------------------------------------------------------------------------
+class FindCellWorklet : public vtkm::worklet::WorkletMapField
+{
+public:
+  using ControlSignature = void(FieldIn points,
+                                ExecObject locator,
+                                FieldOut cellIds,
+                                FieldOut pcoords);
+  using ExecutionSignature = void(_1, _2, _3, _4);
+
+  template <typename LocatorType>
+  VTKM_EXEC void operator()(const vtkm::Vec3f& point,
+                            const LocatorType& locator,
+                            vtkm::Id& cellId,
+                            vtkm::Vec3f& pcoords) const
+  {
+    vtkm::ErrorCode status = locator.FindCell(point, cellId, pcoords);
+    if (status != vtkm::ErrorCode::Success)
+    {
+      this->RaiseError(vtkm::ErrorString(status));
+    }
+  }
+};
+
+template <typename CellSetType, typename CoordinateSystemArrayType>
+void TestWithDataSet(const vtkm::cont::DataSet& dataset)
+{
+  VTKM_TEST_ASSERT(dataset.GetCellSet().IsType<CellSetType>());
+  VTKM_TEST_ASSERT(dataset.GetCoordinateSystem().GetData().IsType<CoordinateSystemArrayType>());
+
+  vtkm::cont::CellLocatorChooser<CellSetType, CoordinateSystemArrayType> locator;
+  locator.SetCellSet(dataset.GetCellSet());
+  locator.SetCoordinates(dataset.GetCoordinateSystem());
+  locator.Update();
+
+  vtkm::cont::ArrayHandle<vtkm::Id> expCellIds;
+  vtkm::cont::ArrayHandle<PointType> expPCoords;
+  vtkm::cont::ArrayHandle<PointType> points;
+  GenerateRandomInput(dataset, 128, expCellIds, expPCoords, points);
+
+  vtkm::cont::ArrayHandle<vtkm::Id> cellIds;
+  vtkm::cont::ArrayHandle<PointType> pcoords;
+
+  vtkm::worklet::DispatcherMapField<FindCellWorklet> dispatcher;
+  dispatcher.Invoke(points, locator, cellIds, pcoords);
+
+  auto cellIdPortal = cellIds.ReadPortal();
+  auto expCellIdsPortal = expCellIds.ReadPortal();
+  auto pcoordsPortal = pcoords.ReadPortal();
+  auto expPCoordsPortal = expPCoords.ReadPortal();
+  for (vtkm::Id i = 0; i < 128; ++i)
+  {
+    VTKM_TEST_ASSERT(cellIdPortal.Get(i) == expCellIdsPortal.Get(i), "Incorrect cell ids");
+    VTKM_TEST_ASSERT(test_equal(pcoordsPortal.Get(i), expPCoordsPortal.Get(i), 1e-3),
+                     "Incorrect parameteric coordinates");
+  }
+
+  std::cout << "Passed.\n";
+}
+
+void TestCellLocatorChooser()
+{
+  std::cout << "Test UniformGrid dataset\n";
+  TestWithDataSet<vtkm::cont::CellSetStructured<3>, vtkm::cont::ArrayHandleUniformPointCoordinates>(
+    MakeTestDataSetUniform());
+
+  std::cout << "Test Rectilinear dataset\n";
+  TestWithDataSet<
+    vtkm::cont::CellSetStructured<3>,
+    vtkm::cont::ArrayHandleCartesianProduct<vtkm::cont::ArrayHandle<vtkm::FloatDefault>,
+                                            vtkm::cont::ArrayHandle<vtkm::FloatDefault>,
+                                            vtkm::cont::ArrayHandle<vtkm::FloatDefault>>>(
+    MakeTestDataSetRectilinear());
+
+  std::cout << "Test Curvilinear dataset\n";
+  TestWithDataSet<vtkm::cont::CellSetStructured<3>, vtkm::cont::ArrayHandle<PointType>>(
+    MakeTestDataSetCurvilinear());
+}
+
+} // anonymous namespace
+
+int UnitTestCellLocatorChooser(int argc, char* argv[])
+{
+  vtkm::cont::GetRuntimeDeviceTracker().ForceDevice(vtkm::cont::DeviceAdapterTagSerial{});
+  return vtkm::cont::testing::Testing::Run(TestCellLocatorChooser, argc, argv);
+}