1 //============================================================================
2 //  Copyright (c) Kitware, Inc.
3 //  All rights reserved.
4 //  See LICENSE.txt for details.
5 //
6 //  This software is distributed WITHOUT ANY WARRANTY; without even
7 //  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
8 //  PURPOSE.  See the above copyright notice for more information.
9 //============================================================================
10 #ifndef vtk_m_exec_CellLocatorBoundingIntervalHierarchy_h
11 #define vtk_m_exec_CellLocatorBoundingIntervalHierarchy_h
12 
13 #include <vtkm/TopologyElementTag.h>
14 #include <vtkm/VecFromPortalPermute.h>
15 #include <vtkm/cont/ArrayHandle.h>
16 #include <vtkm/cont/CoordinateSystem.h>
17 #include <vtkm/exec/CellInside.h>
18 #include <vtkm/exec/ParametricCoordinates.h>
19 
20 namespace vtkm
21 {
22 namespace exec
23 {
24 
25 
26 
27 
28 struct CellLocatorBoundingIntervalHierarchyNode
29 {
30 #if defined(VTKM_CLANG)
31 #pragma GCC diagnostic push
32 #pragma GCC diagnostic ignored "-Wnested-anon-types"
33 #endif // gcc || clang
34   vtkm::IdComponent Dimension;
35   vtkm::Id ParentIndex;
36   vtkm::Id ChildIndex;
37   union {
38     struct
39     {
40       vtkm::FloatDefault LMax;
41       vtkm::FloatDefault RMin;
42     } Node;
43     struct
44     {
45       vtkm::Id Start;
46       vtkm::Id Size;
47     } Leaf;
48   };
49 #if defined(VTKM_CLANG)
50 #pragma GCC diagnostic pop
51 #endif // gcc || clang
52 
53   VTKM_EXEC_CONT
CellLocatorBoundingIntervalHierarchyNodeCellLocatorBoundingIntervalHierarchyNode54   CellLocatorBoundingIntervalHierarchyNode()
55     : Dimension()
56     , ParentIndex()
57     , ChildIndex()
58     , Node{ 0, 0 }
59   {
60   }
61 }; // struct CellLocatorBoundingIntervalHierarchyNode
62 
63 template <typename CellSetType>
64 class VTKM_ALWAYS_EXPORT CellLocatorBoundingIntervalHierarchy
65 {
66   using NodeArrayHandle =
67     vtkm::cont::ArrayHandle<vtkm::exec::CellLocatorBoundingIntervalHierarchyNode>;
68   using CellIdArrayHandle = vtkm::cont::ArrayHandle<vtkm::Id>;
69 
70 public:
71   VTKM_CONT
CellLocatorBoundingIntervalHierarchy(const NodeArrayHandle & nodes,const CellIdArrayHandle & cellIds,const CellSetType & cellSet,const vtkm::cont::CoordinateSystem::MultiplexerArrayType & coords,vtkm::cont::DeviceAdapterId device,vtkm::cont::Token & token)72   CellLocatorBoundingIntervalHierarchy(
73     const NodeArrayHandle& nodes,
74     const CellIdArrayHandle& cellIds,
75     const CellSetType& cellSet,
76     const vtkm::cont::CoordinateSystem::MultiplexerArrayType& coords,
77     vtkm::cont::DeviceAdapterId device,
78     vtkm::cont::Token& token)
79     : Nodes(nodes.PrepareForInput(device, token))
80     , CellIds(cellIds.PrepareForInput(device, token))
81     , CellSet(cellSet.PrepareForInput(device, VisitType(), IncidentType(), token))
82     , Coords(coords.PrepareForInput(device, token))
83   {
84   }
85 
86 
87   VTKM_EXEC
FindCell(const vtkm::Vec3f & point,vtkm::Id & cellId,vtkm::Vec3f & parametric)88   vtkm::ErrorCode FindCell(const vtkm::Vec3f& point,
89                            vtkm::Id& cellId,
90                            vtkm::Vec3f& parametric) const
91   {
92     cellId = -1;
93     vtkm::Id nodeIndex = 0;
94     FindCellState state = FindCellState::EnterNode;
95 
96     while ((cellId < 0) && !((nodeIndex == 0) && (state == FindCellState::AscendFromNode)))
97     {
98       switch (state)
99       {
100         case FindCellState::EnterNode:
101           VTKM_RETURN_ON_ERROR(this->EnterNode(state, point, cellId, nodeIndex, parametric));
102           break;
103         case FindCellState::AscendFromNode:
104           this->AscendFromNode(state, nodeIndex);
105           break;
106         case FindCellState::DescendLeftChild:
107           this->DescendLeftChild(state, point, nodeIndex);
108           break;
109         case FindCellState::DescendRightChild:
110           this->DescendRightChild(state, point, nodeIndex);
111           break;
112       }
113     }
114 
115     if (cellId >= 0)
116     {
117       return vtkm::ErrorCode::Success;
118     }
119     else
120     {
121       return vtkm::ErrorCode::CellNotFound;
122     }
123   }
124 
125   VTKM_DEPRECATED(1.6, "Locators are no longer pointers. Use . operator.")
126   VTKM_EXEC CellLocatorBoundingIntervalHierarchy* operator->() { return this; }
127   VTKM_DEPRECATED(1.6, "Locators are no longer pointers. Use . operator.")
128   VTKM_EXEC const CellLocatorBoundingIntervalHierarchy* operator->() const { return this; }
129 
130 private:
131   enum struct FindCellState
132   {
133     EnterNode,
134     AscendFromNode,
135     DescendLeftChild,
136     DescendRightChild
137   };
138 
139   VTKM_EXEC
EnterNode(FindCellState & state,const vtkm::Vec3f & point,vtkm::Id & cellId,vtkm::Id nodeIndex,vtkm::Vec3f & parametric)140   vtkm::ErrorCode EnterNode(FindCellState& state,
141                             const vtkm::Vec3f& point,
142                             vtkm::Id& cellId,
143                             vtkm::Id nodeIndex,
144                             vtkm::Vec3f& parametric) const
145   {
146     VTKM_ASSERT(state == FindCellState::EnterNode);
147 
148     const vtkm::exec::CellLocatorBoundingIntervalHierarchyNode& node = this->Nodes.Get(nodeIndex);
149 
150     if (node.ChildIndex < 0)
151     {
152       // In a leaf node. Look for a containing cell.
153       VTKM_RETURN_ON_ERROR(this->FindInLeaf(point, parametric, node, cellId));
154       state = FindCellState::AscendFromNode;
155     }
156     else
157     {
158       state = FindCellState::DescendLeftChild;
159     }
160     return vtkm::ErrorCode::Success;
161   }
162 
163   VTKM_EXEC
AscendFromNode(FindCellState & state,vtkm::Id & nodeIndex)164   void AscendFromNode(FindCellState& state, vtkm::Id& nodeIndex) const
165   {
166     VTKM_ASSERT(state == FindCellState::AscendFromNode);
167 
168     vtkm::Id childNodeIndex = nodeIndex;
169     const vtkm::exec::CellLocatorBoundingIntervalHierarchyNode& childNode =
170       this->Nodes.Get(childNodeIndex);
171     nodeIndex = childNode.ParentIndex;
172     const vtkm::exec::CellLocatorBoundingIntervalHierarchyNode& parentNode =
173       this->Nodes.Get(nodeIndex);
174 
175     if (parentNode.ChildIndex == childNodeIndex)
176     {
177       // Ascending from left child. Descend into the right child.
178       state = FindCellState::DescendRightChild;
179     }
180     else
181     {
182       VTKM_ASSERT(parentNode.ChildIndex + 1 == childNodeIndex);
183       // Ascending from right child. Ascend again. (Don't need to change state.)
184     }
185   }
186 
187   VTKM_EXEC
DescendLeftChild(FindCellState & state,const vtkm::Vec3f & point,vtkm::Id & nodeIndex)188   void DescendLeftChild(FindCellState& state, const vtkm::Vec3f& point, vtkm::Id& nodeIndex) const
189   {
190     VTKM_ASSERT(state == FindCellState::DescendLeftChild);
191 
192     const vtkm::exec::CellLocatorBoundingIntervalHierarchyNode& node = this->Nodes.Get(nodeIndex);
193     const vtkm::FloatDefault& coordinate = point[node.Dimension];
194     if (coordinate <= node.Node.LMax)
195     {
196       // Left child does contain the point. Do the actual descent.
197       nodeIndex = node.ChildIndex;
198       state = FindCellState::EnterNode;
199     }
200     else
201     {
202       // Left child does not contain the point. Skip to the right child.
203       state = FindCellState::DescendRightChild;
204     }
205   }
206 
207   VTKM_EXEC
DescendRightChild(FindCellState & state,const vtkm::Vec3f & point,vtkm::Id & nodeIndex)208   void DescendRightChild(FindCellState& state, const vtkm::Vec3f& point, vtkm::Id& nodeIndex) const
209   {
210     VTKM_ASSERT(state == FindCellState::DescendRightChild);
211 
212     const vtkm::exec::CellLocatorBoundingIntervalHierarchyNode& node = this->Nodes.Get(nodeIndex);
213     const vtkm::FloatDefault& coordinate = point[node.Dimension];
214     if (coordinate >= node.Node.RMin)
215     {
216       // Right child does contain the point. Do the actual descent.
217       nodeIndex = node.ChildIndex + 1;
218       state = FindCellState::EnterNode;
219     }
220     else
221     {
222       // Right child does not contain the point. Skip to ascent
223       state = FindCellState::AscendFromNode;
224     }
225   }
226 
FindInLeaf(const vtkm::Vec3f & point,vtkm::Vec3f & parametric,const vtkm::exec::CellLocatorBoundingIntervalHierarchyNode & node,vtkm::Id & containingCellId)227   VTKM_EXEC vtkm::ErrorCode FindInLeaf(
228     const vtkm::Vec3f& point,
229     vtkm::Vec3f& parametric,
230     const vtkm::exec::CellLocatorBoundingIntervalHierarchyNode& node,
231     vtkm::Id& containingCellId) const
232   {
233     using IndicesType = typename CellSetPortal::IndicesType;
234     for (vtkm::Id i = node.Leaf.Start; i < node.Leaf.Start + node.Leaf.Size; ++i)
235     {
236       vtkm::Id cellId = this->CellIds.Get(i);
237       IndicesType cellPointIndices = this->CellSet.GetIndices(cellId);
238       vtkm::VecFromPortalPermute<IndicesType, CoordsPortal> cellPoints(&cellPointIndices,
239                                                                        this->Coords);
240       bool found;
241       VTKM_RETURN_ON_ERROR(this->IsPointInCell(
242         point, parametric, this->CellSet.GetCellShape(cellId), cellPoints, found));
243       if (found)
244       {
245         containingCellId = cellId;
246         return vtkm::ErrorCode::Success;
247       }
248     }
249     containingCellId = -1;
250     return vtkm::ErrorCode::Success;
251   }
252 
253   template <typename CoordsType, typename CellShapeTag>
IsPointInCell(const vtkm::Vec3f & point,vtkm::Vec3f & parametric,CellShapeTag cellShape,const CoordsType & cellPoints,bool & isInside)254   VTKM_EXEC static vtkm::ErrorCode IsPointInCell(const vtkm::Vec3f& point,
255                                                  vtkm::Vec3f& parametric,
256                                                  CellShapeTag cellShape,
257                                                  const CoordsType& cellPoints,
258                                                  bool& isInside)
259   {
260     isInside = false;
261     VTKM_RETURN_ON_ERROR(vtkm::exec::WorldCoordinatesToParametricCoordinates(
262       cellPoints, point, cellShape, parametric));
263     isInside = vtkm::exec::CellInside(parametric, cellShape);
264     return vtkm::ErrorCode::Success;
265   }
266 
267   using VisitType = vtkm::TopologyElementTagCell;
268   using IncidentType = vtkm::TopologyElementTagPoint;
269   using NodePortal = typename NodeArrayHandle::ReadPortalType;
270   using CellIdPortal = typename CellIdArrayHandle::ReadPortalType;
271   using CellSetPortal =
272     typename CellSetType::template ExecConnectivityType<VisitType, IncidentType>;
273   using CoordsPortal = typename vtkm::cont::CoordinateSystem::MultiplexerArrayType::ReadPortalType;
274 
275   NodePortal Nodes;
276   CellIdPortal CellIds;
277   CellSetPortal CellSet;
278   CoordsPortal Coords;
279 }; // class CellLocatorBoundingIntervalHierarchy
280 
281 } // namespace exec
282 
283 } // namespace vtkm
284 
285 #endif //vtk_m_exec_CellLocatorBoundingIntervalHierarchy_h
286