CellLocatorBoundingIntervalHierarchyExec.h 8.79 KB
Newer Older
1
//============================================================================
2 3 4
//  Copyright (c) Kitware, Inc.
//  All rights reserved.
//  See LICENSE.txt for details.
5
//
6 7 8 9
//  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.
//============================================================================
10 11
#ifndef vtk_m_cont_CellLocatorBoundingIntervalHierarchyExec_h
#define vtk_m_cont_CellLocatorBoundingIntervalHierarchyExec_h
12

13
#include <vtkm/TopologyElementTag.h>
14 15
#include <vtkm/VecFromPortalPermute.h>
#include <vtkm/cont/ArrayHandle.h>
16
#include <vtkm/cont/ArrayHandleVirtualCoordinates.h>
17
#include <vtkm/exec/CellInside.h>
18
#include <vtkm/exec/CellLocator.h>
19 20 21 22 23 24
#include <vtkm/exec/ParametricCoordinates.h>

namespace vtkm
{
namespace exec
{
25 26

struct CellLocatorBoundingIntervalHierarchyNode
27
{
28 29 30 31 32 33 34 35 36 37 38 39 40 41 42
  vtkm::IdComponent Dimension;
  vtkm::Id ParentIndex;
  vtkm::Id ChildIndex;
  union {
    struct
    {
      vtkm::FloatDefault LMax;
      vtkm::FloatDefault RMin;
    } Node;
    struct
    {
      vtkm::Id Start;
      vtkm::Id Size;
    } Leaf;
  };
43

44 45 46 47 48 49 50 51 52
  VTKM_EXEC_CONT
  CellLocatorBoundingIntervalHierarchyNode()
    : Dimension()
    , ParentIndex()
    , ChildIndex()
    , Node{ 0, 0 }
  {
  }
}; // struct CellLocatorBoundingIntervalHierarchyNode
53

54
template <typename DeviceAdapter, typename CellSetType>
55
class CellLocatorBoundingIntervalHierarchyExec : public vtkm::exec::CellLocator
56
{
57 58 59 60
  using NodeArrayHandle =
    vtkm::cont::ArrayHandle<vtkm::exec::CellLocatorBoundingIntervalHierarchyNode>;
  using CellIdArrayHandle = vtkm::cont::ArrayHandle<vtkm::Id>;

61 62
public:
  VTKM_CONT
63
  CellLocatorBoundingIntervalHierarchyExec() {}
64 65

  VTKM_CONT
66 67 68 69 70
  CellLocatorBoundingIntervalHierarchyExec(const NodeArrayHandle& nodes,
                                           const CellIdArrayHandle& cellIds,
                                           const CellSetType& cellSet,
                                           const vtkm::cont::ArrayHandleVirtualCoordinates& coords,
                                           DeviceAdapter)
71 72
    : Nodes(nodes.PrepareForInput(DeviceAdapter()))
    , CellIds(cellIds.PrepareForInput(DeviceAdapter()))
73 74
    , CellSet(cellSet.PrepareForInput(DeviceAdapter(), FromType(), ToType()))
    , Coords(coords.PrepareForInput(DeviceAdapter()))
75 76 77
  {
  }

78 79 80 81 82
  VTKM_EXEC
  void FindCell(const vtkm::Vec<vtkm::FloatDefault, 3>& point,
                vtkm::Id& cellId,
                vtkm::Vec<vtkm::FloatDefault, 3>& parametric,
                const vtkm::exec::FunctorBase& worklet) const override
83
  {
84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105
    cellId = -1;
    vtkm::Id nodeIndex = 0;
    FindCellState state = FindCellState::EnterNode;

    while ((cellId < 0) && !((nodeIndex == 0) && (state == FindCellState::AscendFromNode)))
    {
      switch (state)
      {
        case FindCellState::EnterNode:
          this->EnterNode(state, point, cellId, nodeIndex, parametric, worklet);
          break;
        case FindCellState::AscendFromNode:
          this->AscendFromNode(state, nodeIndex);
          break;
        case FindCellState::DescendLeftChild:
          this->DescendLeftChild(state, point, nodeIndex);
          break;
        case FindCellState::DescendRightChild:
          this->DescendRightChild(state, point, nodeIndex);
          break;
      }
    }
106 107 108
  }

private:
109
  enum struct FindCellState
110
  {
111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126
    EnterNode,
    AscendFromNode,
    DescendLeftChild,
    DescendRightChild
  };

  VTKM_EXEC
  void EnterNode(FindCellState& state,
                 const vtkm::Vec<vtkm::FloatDefault, 3>& point,
                 vtkm::Id& cellId,
                 vtkm::Id nodeIndex,
                 vtkm::Vec<vtkm::FloatDefault, 3>& parametric,
                 const vtkm::exec::FunctorBase& worklet) const
  {
    VTKM_ASSERT(state == FindCellState::EnterNode);

127
    const vtkm::exec::CellLocatorBoundingIntervalHierarchyNode& node = this->Nodes.Get(nodeIndex);
128

129 130
    if (node.ChildIndex < 0)
    {
131 132 133
      // In a leaf node. Look for a containing cell.
      cellId = this->FindInLeaf(point, parametric, node, worklet);
      state = FindCellState::AscendFromNode;
134 135 136
    }
    else
    {
137 138 139 140 141 142 143 144 145 146
      state = FindCellState::DescendLeftChild;
    }
  }

  VTKM_EXEC
  void AscendFromNode(FindCellState& state, vtkm::Id& nodeIndex) const
  {
    VTKM_ASSERT(state == FindCellState::AscendFromNode);

    vtkm::Id childNodeIndex = nodeIndex;
147 148
    const vtkm::exec::CellLocatorBoundingIntervalHierarchyNode& childNode =
      this->Nodes.Get(childNodeIndex);
149
    nodeIndex = childNode.ParentIndex;
150 151
    const vtkm::exec::CellLocatorBoundingIntervalHierarchyNode& parentNode =
      this->Nodes.Get(nodeIndex);
152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171

    if (parentNode.ChildIndex == childNodeIndex)
    {
      // Ascending from left child. Descend into the right child.
      state = FindCellState::DescendRightChild;
    }
    else
    {
      VTKM_ASSERT(parentNode.ChildIndex + 1 == childNodeIndex);
      // Ascending from right child. Ascend again. (Don't need to change state.)
    }
  }

  VTKM_EXEC
  void DescendLeftChild(FindCellState& state,
                        const vtkm::Vec<vtkm::FloatDefault, 3>& point,
                        vtkm::Id& nodeIndex) const
  {
    VTKM_ASSERT(state == FindCellState::DescendLeftChild);

172
    const vtkm::exec::CellLocatorBoundingIntervalHierarchyNode& node = this->Nodes.Get(nodeIndex);
173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193
    const vtkm::FloatDefault& coordinate = point[node.Dimension];
    if (coordinate <= node.Node.LMax)
    {
      // Left child does contain the point. Do the actual descent.
      nodeIndex = node.ChildIndex;
      state = FindCellState::EnterNode;
    }
    else
    {
      // Left child does not contain the point. Skip to the right child.
      state = FindCellState::DescendRightChild;
    }
  }

  VTKM_EXEC
  void DescendRightChild(FindCellState& state,
                         const vtkm::Vec<vtkm::FloatDefault, 3>& point,
                         vtkm::Id& nodeIndex) const
  {
    VTKM_ASSERT(state == FindCellState::DescendRightChild);

194
    const vtkm::exec::CellLocatorBoundingIntervalHierarchyNode& node = this->Nodes.Get(nodeIndex);
195 196 197 198 199 200 201 202 203 204 205
    const vtkm::FloatDefault& coordinate = point[node.Dimension];
    if (coordinate >= node.Node.RMin)
    {
      // Right child does contain the point. Do the actual descent.
      nodeIndex = node.ChildIndex + 1;
      state = FindCellState::EnterNode;
    }
    else
    {
      // Right child does not contain the point. Skip to ascent
      state = FindCellState::AscendFromNode;
206 207 208
    }
  }

209 210
  VTKM_EXEC vtkm::Id FindInLeaf(const vtkm::Vec<vtkm::FloatDefault, 3>& point,
                                vtkm::Vec<vtkm::FloatDefault, 3>& parametric,
211
                                const vtkm::exec::CellLocatorBoundingIntervalHierarchyNode& node,
212 213
                                const vtkm::exec::FunctorBase& worklet) const
  {
214
    using IndicesType = typename CellSetPortal::IndicesType;
215 216
    for (vtkm::Id i = node.Leaf.Start; i < node.Leaf.Start + node.Leaf.Size; ++i)
    {
217 218 219 220 221
      vtkm::Id cellId = this->CellIds.Get(i);
      IndicesType cellPointIndices = this->CellSet.GetIndices(cellId);
      vtkm::VecFromPortalPermute<IndicesType, CoordsPortal> cellPoints(&cellPointIndices,
                                                                       this->Coords);
      if (IsPointInCell(point, parametric, this->CellSet.GetCellShape(cellId), cellPoints, worklet))
222 223 224 225 226 227 228 229
      {
        return cellId;
      }
    }
    return -1;
  }

  template <typename CoordsType, typename CellShapeTag>
230 231
  VTKM_EXEC static bool IsPointInCell(const vtkm::Vec<vtkm::FloatDefault, 3>& point,
                                      vtkm::Vec<vtkm::FloatDefault, 3>& parametric,
232 233 234 235 236
                                      CellShapeTag cellShape,
                                      const CoordsType& cellPoints,
                                      const vtkm::exec::FunctorBase& worklet)
  {
    bool success = false;
237 238 239
    parametric = vtkm::exec::WorldCoordinatesToParametricCoordinates(
      cellPoints, point, cellShape, success, worklet);
    return success && vtkm::exec::CellInside(parametric, cellShape);
240 241
  }

242 243
  using FromType = vtkm::TopologyElementTagPoint;
  using ToType = vtkm::TopologyElementTagCell;
244 245 246
  using NodePortal = typename NodeArrayHandle::template ExecutionTypes<DeviceAdapter>::PortalConst;
  using CellIdPortal =
    typename CellIdArrayHandle::template ExecutionTypes<DeviceAdapter>::PortalConst;
247 248 249 250
  using CellSetPortal =
    typename CellSetType::template ExecutionTypes<DeviceAdapter, FromType, ToType>::ExecObjectType;
  using CoordsPortal = typename vtkm::cont::ArrayHandleVirtualCoordinates::template ExecutionTypes<
    DeviceAdapter>::PortalConst;
251 252 253

  NodePortal Nodes;
  CellIdPortal CellIds;
254 255
  CellSetPortal CellSet;
  CoordsPortal Coords;
256
}; // class CellLocatorBoundingIntervalHierarchyExec
257 258 259 260 261

} // namespace exec

} // namespace vtkm

262
#endif //vtk_m_cont_CellLocatorBoundingIntervalHierarchyExec_h