Adding updated CellLocator / PointLocator APIs

- Adding API files
- Adding back Manish's BoundingIntervalHierarchy search structure
- Updating CMakeLists.txt to accomodate these changes
- Adding the old test file from Manish - won't build for now

vtkm/cont/BoundingIntervalHierarchyExec.h

@@ -0,0 +1,153 @@
+//============================================================================
+////  Copyright (c) Kitware, Inc.
+////  All rights reserved.
+////  See LICENSE.txt for details.
+////  This software is distributed WITHOUT ANY WARRANTY; without even
+////  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+////  PURPOSE.  See the above copyright notice for more information.
+////
+////  Copyright 2015 Sandia Corporation.
+////  Copyright 2015 UT-Battelle, LLC.
+////  Copyright 2015 Los Alamos National Security.
+////
+////  Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
+////  the U.S. Government retains certain rights in this software.
+////
+////  Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National
+////  Laboratory (LANL), the U.S. Government retains certain rights in
+////  this software.
+////============================================================================
+#ifndef vtk_m_cont_BoundingIntervalHierarchyExec_h
+#define vtk_m_cont_BoundingIntervalHierarchyExec_h
+
+#include <vtkm/VecFromPortalPermute.h>
+#include <vtkm/cont/ArrayHandle.h>
+#include <vtkm/cont/BoundingIntervalHierarchyNode.h>
+#include <vtkm/cont/CellLocator.h>
+#include <vtkm/exec/CellInside.h>
+#include <vtkm/exec/ParametricCoordinates.h>
+
+namespace vtkm
+{
+namespace exec
+{
+namespace
+{
+using NodeArrayHandle = vtkm::cont::ArrayHandle<vtkm::cont::BoundingIntervalHierarchyNode>;
+using CellIdArrayHandle = vtkm::cont::ArrayHandle<vtkm::Id>;
+
+} // namespace
+
+template <typename DeviceAdapter>
+class BoundingIntervalHierarchyExec : public vtkm::exec::CellLocator
+{
+public:
+  VTKM_CONT
+  BoundingIntervalHierarchyExec() {}
+
+  VTKM_CONT
+  BoundingIntervalHierarchyExec(const NodeArrayHandle& nodes,
+                                const CellIdArrayHandle& cellIds,
+                                DeviceAdapter)
+    : Nodes(nodes.PrepareForInput(DeviceAdapter()))
+    , CellIds(cellIds.PrepareForInput(DeviceAdapter()))
+  {
+  }
+
+  template <typename CellSetType, typename PointPortal>
+  VTKM_EXEC vtkm::Id Find(const vtkm::Vec<vtkm::Float64, 3>& point,
+                          const CellSetType& cellSet,
+                          const PointPortal& points,
+                          const vtkm::exec::FunctorBase& worklet) const
+  {
+    return Find(0, point, cellSet, points, worklet);
+  }
+
+private:
+  template <typename CellSetType, typename PointPortal>
+  VTKM_EXEC vtkm::Id Find(vtkm::Id index,
+                          const vtkm::Vec<vtkm::Float64, 3>& point,
+                          const CellSetType& cellSet,
+                          const PointPortal& points,
+                          const vtkm::exec::FunctorBase& worklet) const
+  {
+    const vtkm::cont::BoundingIntervalHierarchyNode& node = Nodes.Get(index);
+    if (node.ChildIndex < 0)
+    {
+      return FindInLeaf(point, node, cellSet, points, worklet);
+    }
+    else
+    {
+      const vtkm::Float64& c = point[node.Dimension];
+      vtkm::Id id1 = -1;
+      vtkm::Id id2 = -1;
+      if (c <= node.Node.LMax)
+      {
+        id1 = Find(node.ChildIndex, point, cellSet, points, worklet);
+      }
+      if (id1 == -1 && c >= node.Node.RMin)
+      {
+        id2 = Find(node.ChildIndex + 1, point, cellSet, points, worklet);
+      }
+      if (id1 == -1 && id2 == -1)
+      {
+        return -1;
+      }
+      else if (id1 == -1)
+      {
+        return id2;
+      }
+      else
+      {
+        return id1;
+      }
+    }
+  }
+
+  template <typename CellSetType, typename PointPortal>
+  VTKM_EXEC vtkm::Id FindInLeaf(const vtkm::Vec<vtkm::Float64, 3>& point,
+                                const vtkm::cont::BoundingIntervalHierarchyNode& node,
+                                const CellSetType& cellSet,
+                                const PointPortal& points,
+                                const vtkm::exec::FunctorBase& worklet) const
+  {
+    using IndicesType = typename CellSetType::IndicesType;
+    for (vtkm::Id i = node.Leaf.Start; i < node.Leaf.Start + node.Leaf.Size; ++i)
+    {
+      vtkm::Id cellId = CellIds.Get(i);
+      IndicesType cellPointIndices = cellSet.GetIndices(cellId);
+      vtkm::VecFromPortalPermute<IndicesType, PointPortal> cellPoints(&cellPointIndices, points);
+      if (IsPointInCell(point, cellSet.GetCellShape(cellId), cellPoints, worklet))
+      {
+        return cellId;
+      }
+    }
+    return -1;
+  }
+
+  template <typename CoordsType, typename CellShapeTag>
+  VTKM_EXEC static bool IsPointInCell(const vtkm::Vec<vtkm::Float64, 3>& point,
+                                      CellShapeTag cellShape,
+                                      const CoordsType& cellPoints,
+                                      const vtkm::exec::FunctorBase& worklet)
+  {
+    bool success = false;
+    vtkm::Vec<vtkm::Float64, 3> parametricCoords =
+      vtkm::exec::WorldCoordinatesToParametricCoordinates(
+        cellPoints, point, cellShape, success, worklet);
+    return success && vtkm::exec::CellInside(parametricCoords, cellShape);
+  }
+
+  using NodePortal = typename NodeArrayHandle::template ExecutionTypes<DeviceAdapter>::PortalConst;
+  using CellIdPortal =
+    typename CellIdArrayHandle::template ExecutionTypes<DeviceAdapter>::PortalConst;
+
+  NodePortal Nodes;
+  CellIdPortal CellIds;
+}; // class BoundingIntervalHierarchyExec
+
+} // namespace exec
+
+} // namespace vtkm
+
+#endif //vtk_m_cont_BoundingIntervalHierarchyExec_h

vtkm/cont/BoundingIntervalHierarchyNode.h

@@ -0,0 +1,52 @@
+//============================================================================
+//  Copyright (c) Kitware, Inc.
+//  All rights reserved.
+//  See LICENSE.txt for details.
+//  This software is distributed WITHOUT ANY WARRANTY; without even
+//  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+//  PURPOSE.  See the above copyright notice for more information.
+//
+//  Copyright 2014 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
+//  Copyright 2014 UT-Battelle, LLC.
+//  Copyright 2014 Los Alamos National Security.
+//
+//  Under the terms of Contract DE-NA0003525 with NTESS,
+//  the U.S. Government retains certain rights in this software.
+//
+//  Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National
+//  Laboratory (LANL), the U.S. Government retains certain rights in
+//  this software.
+//============================================================================
+
+#ifndef vtk_m_cont_BoundingIntervalHierarchyNode_h
+#define vtk_m_cont_BoundingIntervalHierarchyNode_h
+
+#include <vtkm/Types.h>
+
+namespace vtkm
+{
+namespace cont
+{
+
+struct BoundingIntervalHierarchyNode
+{
+  vtkm::IdComponent Dimension;
+  vtkm::Id ChildIndex;
+  union {
+    struct
+    {
+      vtkm::Float64 LMax;
+      vtkm::Float64 RMin;
+    } Node;
+    struct
+    {
+      vtkm::Id Start;
+      vtkm::Id Size;
+    } Leaf;
+  };
+}; // struct BoundingIntervalHierarchyNode
+
+} // namespace cont
+} // namespace vtkm
+
+#endif // vtk_m_cont_BoundingIntervalHierarchyNode_h

vtkm/cont/CMakeLists.txt

@@ -46,8 +46,12 @@ set(headers
   ArrayHandleConcatenate.h
   ArrayRangeCompute.h
   AssignerMultiBlock.h
+  BoundingIntervalHierarchyNode.h
+  BoundingIntervalHierarchyExec.h
+  BoundingIntervalHierarchy.h
   BoundsCompute.h
   BoundsGlobalCompute.h
+  CellLocator.h
   CellLocatorHelper.h
   CellLocatorTwoLevelUniformGrid.h
   CellSet.h
@@ -84,6 +88,7 @@ set(headers
   FieldRangeGlobalCompute.h
   ImplicitFunctionHandle.h
   MultiBlock.h
+  PointLocator.h
   PointLocatorUniformGrid.h
   RuntimeDeviceInformation.h
   RuntimeDeviceTracker.h

vtkm/cont/CellLocator.h

@@ -0,0 +1,89 @@
+#ifndef vtk_m_cont_CellLocator_h
+#define vtk_m_cont_CellLocator_h
+
+#include <vtkm/Types.h>
+#include <vtkm/cont/CoordinateSystem.h>
+#include <vtkm/cont/DeviceAdapter.h>
+#include <vtkm/cont/DynamicCellSet.h>
+#include <vtkm/cont/ExecutionObjectBase.h>
+
+namespace vtkm
+{
+
+namespace exec
+{
+// This will actually be used in the Execution Environment.
+// As this object is returned by the PrepareForExecution on
+// the CellLocator we need it to be covarient, and this acts
+// like a base class.
+
+class CellLocator
+{
+public:
+  VTKM_EXEC void FindCell(const vtkm::Vec<vtkm::FloatDefault, 3>& point,
+                          vtkm::Id& cellId,
+                          vtkm::Vec<vtkm::FloatDefault, 3>& parametric)
+  {
+  } //const = 0;
+};
+
+} // namespace exec
+
+namespace cont
+{
+
+class CellLocator : public ExecutionObjectBase
+{
+
+public:
+  CellLocator()
+    : Dirty(true)
+  {
+  }
+
+  vtkm::cont::DynamicCellSet GetCellSet() const { return CellSet; }
+
+  void SetCellSet(const vtkm::cont::DynamicCellSet& cellSet)
+  {
+    CellSet = cellSet;
+    Dirty = true;
+  }
+
+  vtkm::cont::CoordinateSystem GetCoords() const { return Coords; }
+
+  void SetCoords(const vtkm::cont::CoordinateSystem& coords)
+  {
+    Coords = coords;
+    Dirty = true;
+  }
+
+  //This is going to need a TryExecute
+  virtual void Build() = 0;
+
+  void Update()
+  {
+    if (Dirty)
+      Build();
+    Dirty = false;
+  }
+
+  template <typename DeviceAdapter>
+  VTKM_CONT std::unique_ptr<vtkm::exec::CellLocator> PrepareForExecution(DeviceAdapter device)
+  {
+    vtkm::cont::DeviceAdapterId deviceId = vtkm::cont::DeviceAdapterTraits<DeviceAdapter>::GetId();
+    return PrepareForExecution(deviceId);
+  }
+
+  VTKM_CONT virtual std::unique_ptr<vtkm::exec::CellLocator> PrepareForExecution(
+    vtkm::cont::DeviceAdapterId device) = 0;
+
+private:
+  vtkm::cont::DynamicCellSet CellSet;
+  vtkm::cont::CoordinateSystem Coords;
+  bool Dirty;
+};
+
+} // namespace cont
+} // namespace vtkm
+
+#endif // vtk_m_cont_CellLocator_h

vtkm/cont/PointLocator.h

@@ -0,0 +1,67 @@
+#include <vtkm/Types.h>
+#include <vtkm/cont/CoordinateSystem.h>
+#include <vtkm/cont/DeviceAdapter.h>
+#include <vtkm/cont/ExecutionObjectBase.h>
+
+namespace vtkm
+{
+
+namespace exec
+{
+
+class PointLocator
+{
+public:
+  VTKM_EXEC virtual void FindNearestNeighbor(const vtkm::Vec<vtkm::FloatDefault, 3>& queryPoint,
+                                             vtkm::Id& pointId,
+                                             vtkm::FloatDefault& distanceSquared) const = 0;
+};
+
+} // namespace exec
+
+namespace cont
+{
+
+class PointLocator : public ExecutionObjectBase
+{
+
+public:
+  PointLocator()
+    : Dirty(true)
+  {
+  }
+
+  vtkm::cont::CoordinateSystem GetCoords() const { return Coords; }
+
+  void SetCoords(const vtkm::cont::CoordinateSystem& coords)
+  {
+    Coords = coords;
+    Dirty = true;
+  }
+
+  virtual void Build() = 0;
+
+  void Update()
+  {
+    if (Dirty)
+      Build();
+    Dirty = false;
+  }
+
+  template <typename DeviceAdapter>
+  VTKM_CONT std::unique_ptr<vtkm::exec::PointLocator> PrepareForExecution(DeviceAdapter device)
+  {
+    vtkm::cont::DeviceAdapterId deviceId = vtkm::cont::DeviceAdapterTraits<DeviceAdapter>::GetId();
+    return PrepareForExecution(deviceId);
+  }
+
+  VTKM_CONT virtual std::unique_ptr<vtkm::exec::PointLocator> PrepareForExecution(
+    vtkm::cont::DeviceAdapterId device) = 0;
+
+private:
+  vtkm::cont::CoordinateSystem Coords;
+  bool Dirty;
+};
+
+} // namespace cont
+} // namespace vtkm

vtkm/cont/testing/UnitTestBoundingIntervalHierarchy.cxx

@@ -0,0 +1,153 @@
+//============================================================================
+//  Copyright (c) Kitware, Inc.
+//  All rights reserved.
+//  See LICENSE.txt for details.
+//  This software is distributed WITHOUT ANY WARRANTY; without even
+//  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+//  PURPOSE.  See the above copyright notice for more information.
+//
+//  Copyright 2014 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
+//  Copyright 2014 UT-Battelle, LLC.
+//  Copyright 2014 Los Alamos National Security.
+//
+//  Under the terms of Contract DE-NA0003525 with NTESS,
+//  the U.S. Government retains certain rights in this software.
+//
+//  Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National
+//  Laboratory (LANL), the U.S. Government retains certain rights in
+//  this software.
+//============================================================================
+
+#include <vtkm/cont/ArrayHandleConcatenate.h>
+#include <vtkm/cont/DataSetBuilderUniform.h>
+#include <vtkm/cont/Timer.h>
+#include <vtkm/cont/testing/Testing.h>
+#include <vtkm/io/reader/VTKDataSetReader.h>
+#include <vtkm/worklet/spatialstructure/BoundingIntervalHierarchy.h>
+#include <vtkm/worklet/spatialstructure/BoundingIntervalHierarchyBuilder.h>
+
+namespace
+{
+struct CellCentroidCalculator : public vtkm::worklet::WorkletMapPointToCell
+{
+  typedef void ControlSignature(CellSetIn, FieldInPoint<>, FieldOut<>);
+  typedef _3 ExecutionSignature(_1, PointCount, _2);
+
+  template <typename CellShape, typename InputPointField>
+  VTKM_EXEC typename InputPointField::ComponentType operator()(
+    CellShape shape,
+    vtkm::IdComponent numPoints,
+    const InputPointField& inputPointField) const
+  {
+    vtkm::Vec<vtkm::FloatDefault, 3> parametricCenter =
+      vtkm::exec::ParametricCoordinatesCenter(numPoints, shape, *this);
+    return vtkm::exec::CellInterpolate(inputPointField, parametricCenter, shape, *this);
+  }
+}; // struct CellCentroidCalculator
+
+struct BoundingIntervalHierarchyTester : public vtkm::worklet::WorkletMapField
+{
+  typedef void ControlSignature(FieldIn<>,
+                                ExecObject,
+                                WholeCellSetIn<>,
+                                WholeArrayIn<>,
+                                FieldIn<>,
+                                FieldOut<>);
+  typedef _6 ExecutionSignature(_1, _2, _3, _4, _5);
+
+  template <typename Point,
+            typename BoundingIntervalHierarchyExecObject,
+            typename CellSet,
+            typename CoordsPortal>
+  VTKM_EXEC vtkm::IdComponent operator()(const Point& point,
+                                         const BoundingIntervalHierarchyExecObject& bih,
+                                         const CellSet& cellSet,
+                                         const CoordsPortal& coords,
+                                         const vtkm::Id expectedId) const
+  {
+    vtkm::Id cellId = bih.Find(point, cellSet, coords, *this);
+    return (1 - static_cast<vtkm::IdComponent>(expectedId == cellId));
+  }
+}; // struct BoundingIntervalHierarchyTester
+
+vtkm::cont::DataSet ConstructDataSet(vtkm::Id size)
+{
+  return vtkm::cont::DataSetBuilderUniform().Create(vtkm::Id3(size, size, size));
+}
+
+void TestBoundingIntervalHierarchy(vtkm::cont::DataSet dataSet, vtkm::IdComponent numPlanes)
+{
+  using DeviceAdapter = VTKM_DEFAULT_DEVICE_ADAPTER_TAG;
+  using Algorithms = vtkm::cont::DeviceAdapterAlgorithm<DeviceAdapter>;
+  using Timer = vtkm::cont::Timer<DeviceAdapter>;
+  namespace spatial = vtkm::worklet::spatialstructure;
+
+  vtkm::cont::DynamicCellSet cellSet = dataSet.GetCellSet();
+  vtkm::cont::DynamicArrayHandleCoordinateSystem vertices = dataSet.GetCoordinateSystem().GetData();
+
+  std::cout << "Using numPlanes: " << numPlanes << "\n";
+  spatial::BoundingIntervalHierarchy bih = spatial::BoundingIntervalHierarchyBuilder(numPlanes, 5)
+                                             .Build(cellSet, vertices, DeviceAdapter());
+
+  Timer centroidsTimer;
+  vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::FloatDefault, 3>> centroids;
+  vtkm::worklet::DispatcherMapTopology<CellCentroidCalculator>().Invoke(
+    cellSet, vertices, centroids);
+  std::cout << "Centroids calculation time: " << centroidsTimer.GetElapsedTime() << "\n";
+
+  vtkm::worklet::spatialstructure::BoundingIntervalHierarchyExecutionObject<DeviceAdapter>
+    bihExecObject = bih.PrepareForInput<DeviceAdapter>();
+
+  vtkm::cont::ArrayHandleCounting<vtkm::Id> expectedCellIds(0, 1, cellSet.GetNumberOfCells());
+
+  Timer interpolationTimer;
+  vtkm::cont::ArrayHandle<vtkm::IdComponent> results;
+#if VTKM_DEVICE_ADAPTER == VTKM_DEVICE_ADAPTER_CUDA
+  //set up stack size for cuda envinroment
+  size_t stackSizeBackup;
+  cudaDeviceGetLimit(&stackSizeBackup, cudaLimitStackSize);
+  cudaDeviceSetLimit(cudaLimitStackSize, 1024 * 200);
+#endif
+  vtkm::worklet::DispatcherMapField<BoundingIntervalHierarchyTester>().Invoke(
+    centroids, bihExecObject, cellSet, vertices, expectedCellIds, results);
+#if VTKM_DEVICE_ADAPTER == VTKM_DEVICE_ADAPTER_CUDA
+  cudaDeviceSetLimit(cudaLimitStackSize, stackSizeBackup);
+#endif
+  vtkm::Id numDiffs = Algorithms::Reduce(results, 0, vtkm::Add());
+  vtkm::Float64 timeDiff = interpolationTimer.GetElapsedTime();
+  std::cout << "No of interpolations: " << results.GetNumberOfValues() << "\n";
+  std::cout << "Interpolation time: " << timeDiff << "\n";
+  std::cout << "Average interpolation rate: "
+            << (static_cast<vtkm::Float64>(results.GetNumberOfValues()) / timeDiff) << "\n";
+  std::cout << "No of diffs: " << numDiffs << "\n";
+}
+
+vtkm::cont::DataSet LoadFromFile(const char* file)
+{
+  vtkm::io::reader::VTKDataSetReader reader(file);
+  return reader.ReadDataSet();
+}
+
+void TestBoundingIntervalHierarchyFromFile(const char* file, vtkm::IdComponent numPlanes)
+{
+  TestBoundingIntervalHierarchy(LoadFromFile(file), numPlanes);
+}
+
+void RunTest()
+{
+  TestBoundingIntervalHierarchy(ConstructDataSet(145), 3);
+  TestBoundingIntervalHierarchy(ConstructDataSet(145), 4);
+  TestBoundingIntervalHierarchy(ConstructDataSet(145), 6);
+  TestBoundingIntervalHierarchy(ConstructDataSet(145), 9);
+  TestBoundingIntervalHierarchyFromFile("buoyancy.vtk", 3);
+  TestBoundingIntervalHierarchyFromFile("buoyancy.vtk", 4);
+  TestBoundingIntervalHierarchyFromFile("buoyancy.vtk", 6);
+  TestBoundingIntervalHierarchyFromFile("buoyancy.vtk", 9);
+}
+
+} // anonymous namespace
+
+int UnitTestBoundingIntervalHierarchy(int, char* [])
+{
+  return vtkm::cont::testing::Testing::Run(RunTest);
+}