VTK-m  2.0
exec/CellLocatorTwoLevel.h
Go to the documentation of this file.
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_CellLocatorTwoLevel_h
11 #define vtk_m_exec_CellLocatorTwoLevel_h
12 
13 #include <vtkm/exec/CellInside.h>
15 
16 #include <vtkm/cont/ArrayHandle.h>
18 
19 #include <vtkm/Math.h>
21 #include <vtkm/Types.h>
23 #include <vtkm/VecTraits.h>
24 
25 namespace vtkm
26 {
27 namespace internal
28 {
29 namespace cl_uniform_bins
30 {
31 
32 using DimensionType = vtkm::Int16;
33 using DimVec3 = vtkm::Vec<DimensionType, 3>;
34 using FloatVec3 = vtkm::Vec3f;
35 
36 struct Grid
37 {
38  DimVec3 Dimensions;
39  // Bug in CUDA 9.2 where having this gap for alignment was for some reason setting garbage
40  // in a union with other cell locators (or perhaps not properly copying data). This appears
41  // to be fixed by CUDA 10.2.
42  DimensionType Padding;
43  FloatVec3 Origin;
44  FloatVec3 BinSize;
45 };
46 
47 struct Bounds
48 {
49  FloatVec3 Min;
50  FloatVec3 Max;
51 };
52 
53 VTKM_EXEC inline vtkm::Id ComputeFlatIndex(const DimVec3& idx, const DimVec3 dim)
54 {
55  return idx[0] + (dim[0] * (idx[1] + (dim[1] * idx[2])));
56 }
57 
58 VTKM_EXEC inline Grid ComputeLeafGrid(const DimVec3& idx, const DimVec3& dim, const Grid& l1Grid)
59 {
60  return { dim,
61  0,
62  l1Grid.Origin + (static_cast<FloatVec3>(idx) * l1Grid.BinSize),
63  l1Grid.BinSize / static_cast<FloatVec3>(dim) };
64 }
65 
66 template <typename PointsVecType>
67 VTKM_EXEC inline Bounds ComputeCellBounds(const PointsVecType& points)
68 {
69  using CoordsType = typename vtkm::VecTraits<PointsVecType>::ComponentType;
71 
72  CoordsType minp = points[0], maxp = points[0];
73  for (vtkm::IdComponent i = 1; i < numPoints; ++i)
74  {
75  minp = vtkm::Min(minp, points[i]);
76  maxp = vtkm::Max(maxp, points[i]);
77  }
78 
79  return { FloatVec3(minp), FloatVec3(maxp) };
80 }
81 }
82 }
83 } // vtkm::internal::cl_uniform_bins
84 
85 namespace vtkm
86 {
87 namespace exec
88 {
89 
90 //--------------------------------------------------------------------
91 template <typename CellStructureType>
93 {
94 private:
97 
98  template <typename T>
100 
101  using CoordsPortalType =
102  typename vtkm::cont::CoordinateSystem::MultiplexerArrayType::ReadPortalType;
103 
104  // TODO: This function may return false positives for non 3D cells as the
105  // tests are done on the projection of the point on the cell. Extra checks
106  // should be added to test if the point actually falls on the cell.
107  template <typename CellShapeTag, typename CoordsType>
109  CellShapeTag cellShape,
110  CoordsType cellPoints,
111  FloatVec3& parametricCoordinates,
112  bool& inside)
113  {
114  auto bounds = vtkm::internal::cl_uniform_bins::ComputeCellBounds(cellPoints);
115  if (point[0] >= bounds.Min[0] && point[0] <= bounds.Max[0] && point[1] >= bounds.Min[1] &&
116  point[1] <= bounds.Max[1] && point[2] >= bounds.Min[2] && point[2] <= bounds.Max[2])
117  {
118  VTKM_RETURN_ON_ERROR(vtkm::exec::WorldCoordinatesToParametricCoordinates(
119  cellPoints, point, cellShape, parametricCoordinates));
120  inside = vtkm::exec::CellInside(parametricCoordinates, cellShape);
121  }
122  else
123  {
124  inside = false;
125  }
126  // Return success error code even point is not inside this cell
128  }
129 
130 public:
131  template <typename CellSetType>
132  VTKM_CONT CellLocatorTwoLevel(const vtkm::internal::cl_uniform_bins::Grid& topLevelGrid,
133  const vtkm::cont::ArrayHandle<DimVec3>& leafDimensions,
134  const vtkm::cont::ArrayHandle<vtkm::Id>& leafStartIndex,
135  const vtkm::cont::ArrayHandle<vtkm::Id>& cellStartIndex,
136  const vtkm::cont::ArrayHandle<vtkm::Id>& cellCount,
137  const vtkm::cont::ArrayHandle<vtkm::Id>& cellIds,
138  const CellSetType& cellSet,
139  const vtkm::cont::CoordinateSystem& coords,
141  vtkm::cont::Token& token)
142  : TopLevel(topLevelGrid)
143  , LeafDimensions(leafDimensions.PrepareForInput(device, token))
144  , LeafStartIndex(leafStartIndex.PrepareForInput(device, token))
145  , CellStartIndex(cellStartIndex.PrepareForInput(device, token))
146  , CellCount(cellCount.PrepareForInput(device, token))
147  , CellIds(cellIds.PrepareForInput(device, token))
148  , CellSet(cellSet.PrepareForInput(device,
151  token))
152  , Coords(coords.GetDataAsMultiplexer().PrepareForInput(device, token))
153  {
154  }
155 
156  struct LastCell
157  {
158  vtkm::Id CellId = -1;
159  vtkm::Id LeafIdx = -1;
160  };
161 
162  VTKM_EXEC
163  vtkm::ErrorCode FindCell(const FloatVec3& point, vtkm::Id& cellId, FloatVec3& parametric) const
164  {
165  LastCell lastCell;
166  return this->FindCellImpl(point, cellId, parametric, lastCell);
167  }
168 
169  VTKM_EXEC
171  vtkm::Id& cellId,
172  FloatVec3& parametric,
173  LastCell& lastCell) const
174  {
175  vtkm::Vec3f pc;
176  //See if point is inside the last cell.
177  if ((lastCell.CellId >= 0) && (lastCell.CellId < this->CellSet.GetNumberOfElements()) &&
178  this->PointInCell(point, lastCell.CellId, pc) == vtkm::ErrorCode::Success)
179  {
180  parametric = pc;
181  cellId = lastCell.CellId;
183  }
184 
185  //See if it's in the last leaf.
186  if ((lastCell.LeafIdx >= 0) && (lastCell.LeafIdx < this->CellCount.GetNumberOfValues()) &&
187  this->PointInLeaf(point, lastCell.LeafIdx, cellId, pc) == vtkm::ErrorCode::Success)
188  {
189  parametric = pc;
190  lastCell.CellId = cellId;
192  }
193 
194  //Call the full point search.
195  return this->FindCellImpl(point, cellId, parametric, lastCell);
196  }
197 
198 private:
199  VTKM_EXEC
201  const vtkm::Id& cid,
202  vtkm::Vec3f& parametric) const
203  {
204  auto indices = this->CellSet.GetIndices(cid);
205  auto pts = vtkm::make_VecFromPortalPermute(&indices, this->Coords);
206  vtkm::Vec3f pc;
207  bool inside;
208  auto status = PointInsideCell(point, this->CellSet.GetCellShape(cid), pts, pc, inside);
209  if (status == vtkm::ErrorCode::Success && inside)
210  {
211  parametric = pc;
213  }
214 
216  }
217 
218  VTKM_EXEC
220  const vtkm::Id& leafIdx,
221  vtkm::Id& cellId,
222  FloatVec3& parametric) const
223  {
224  vtkm::Id start = this->CellStartIndex.Get(leafIdx);
225  vtkm::Id end = start + this->CellCount.Get(leafIdx);
226 
227  for (vtkm::Id i = start; i < end; ++i)
228  {
229  vtkm::Vec3f pc;
230 
231  vtkm::Id cid = this->CellIds.Get(i);
232  if (this->PointInCell(point, cid, pc) == vtkm::ErrorCode::Success)
233  {
234  cellId = cid;
235  parametric = pc;
237  }
238  }
239 
241  }
242 
243 
244  VTKM_EXEC
246  vtkm::Id& cellId,
247  FloatVec3& parametric,
248  LastCell& lastCell) const
249  {
250  using namespace vtkm::internal::cl_uniform_bins;
251 
252  cellId = -1;
253  lastCell.CellId = -1;
254  lastCell.LeafIdx = -1;
255 
256  DimVec3 binId3 = static_cast<DimVec3>((point - this->TopLevel.Origin) / this->TopLevel.BinSize);
257  if (binId3[0] >= 0 && binId3[0] < this->TopLevel.Dimensions[0] && binId3[1] >= 0 &&
258  binId3[1] < this->TopLevel.Dimensions[1] && binId3[2] >= 0 &&
259  binId3[2] < this->TopLevel.Dimensions[2])
260  {
261  vtkm::Id binId = ComputeFlatIndex(binId3, this->TopLevel.Dimensions);
262 
263  auto ldim = this->LeafDimensions.Get(binId);
264  if (!ldim[0] || !ldim[1] || !ldim[2])
265  {
267  }
268 
269  auto leafGrid = ComputeLeafGrid(binId3, ldim, this->TopLevel);
270 
271  DimVec3 leafId3 = static_cast<DimVec3>((point - leafGrid.Origin) / leafGrid.BinSize);
272  // precision issues may cause leafId3 to be out of range so clamp it
273  leafId3 = vtkm::Max(DimVec3(0), vtkm::Min(ldim - DimVec3(1), leafId3));
274 
275  vtkm::Id leafStart = this->LeafStartIndex.Get(binId);
276  vtkm::Id leafIdx = leafStart + ComputeFlatIndex(leafId3, leafGrid.Dimensions);
277 
278  if (this->PointInLeaf(point, leafIdx, cellId, parametric) == vtkm::ErrorCode::Success)
279  {
280  lastCell.CellId = cellId;
281  lastCell.LeafIdx = leafIdx;
283  }
284  }
285 
287  }
288 
289  vtkm::internal::cl_uniform_bins::Grid TopLevel;
290 
293 
297 
298  CellStructureType CellSet;
300 };
301 }
302 } // vtkm::exec
303 
304 #endif //vtk_m_exec_CellLocatorTwoLevel_h
vtkm::TopologyElementTagPoint
A tag used to identify the point elements in a topology.
Definition: TopologyElementTag.h:34
vtkm::exec::CellLocatorTwoLevel
Definition: exec/CellLocatorTwoLevel.h:92
vtkm::cont::ArrayHandle
Manages an array-worth of data.
Definition: ArrayHandle.h:283
ArrayHandle.h
vtkm::ErrorCode
ErrorCode
Definition: ErrorCode.h:19
vtkm::exec::CellLocatorTwoLevel::Coords
CoordsPortalType Coords
Definition: exec/CellLocatorTwoLevel.h:299
vtkm::exec::CellLocatorTwoLevel::CoordsPortalType
typename vtkm::cont::CoordinateSystem::MultiplexerArrayType::ReadPortalType CoordsPortalType
Definition: exec/CellLocatorTwoLevel.h:102
VTKM_EXEC
#define VTKM_EXEC
Definition: ExportMacros.h:51
vtkm
Groups connected points that have the same field value.
Definition: Atomic.h:19
vtkm::make_VecFromPortalPermute
VTKM_EXEC VecFromPortalPermute< IndexVecType, PortalType > make_VecFromPortalPermute(const IndexVecType *index, const PortalType &portal)
Definition: VecFromPortalPermute.h:166
Types.h
vtkm::exec::CellLocatorTwoLevel::FindCell
VTKM_EXEC vtkm::ErrorCode FindCell(const FloatVec3 &point, vtkm::Id &cellId, FloatVec3 &parametric) const
Definition: exec/CellLocatorTwoLevel.h:163
vtkm::exec::CellLocatorTwoLevel::CellSet
CellStructureType CellSet
Definition: exec/CellLocatorTwoLevel.h:298
vtkm::exec::CellLocatorTwoLevel::FindCellImpl
VTKM_EXEC vtkm::ErrorCode FindCellImpl(const FloatVec3 &point, vtkm::Id &cellId, FloatVec3 &parametric, LastCell &lastCell) const
Definition: exec/CellLocatorTwoLevel.h:245
vtkm::IdComponent
vtkm::Int32 IdComponent
Represents a component ID (index of component in a vector).
Definition: Types.h:168
vtkm::exec::CellLocatorTwoLevel::LeafDimensions
ReadPortal< DimVec3 > LeafDimensions
Definition: exec/CellLocatorTwoLevel.h:291
vtkm::exec::CellLocatorTwoLevel::CellStartIndex
ReadPortal< vtkm::Id > CellStartIndex
Definition: exec/CellLocatorTwoLevel.h:294
vtkm::ErrorCode::Success
@ Success
vtkm::exec::CellLocatorTwoLevel::PointInLeaf
VTKM_EXEC vtkm::ErrorCode PointInLeaf(const FloatVec3 &point, const vtkm::Id &leafIdx, vtkm::Id &cellId, FloatVec3 &parametric) const
Definition: exec/CellLocatorTwoLevel.h:219
vtkm::Int16
int16_t Int16
Definition: Types.h:158
vtkm::exec::CellLocatorTwoLevel::LastCell::CellId
vtkm::Id CellId
Definition: exec/CellLocatorTwoLevel.h:158
vtkm::exec::CellLocatorTwoLevel::LastCell::LeafIdx
vtkm::Id LeafIdx
Definition: exec/CellLocatorTwoLevel.h:159
CoordinateSystem.h
vtkm::cont::ArrayHandle::ReadPortalType
typename StorageType::ReadPortalType ReadPortalType
Definition: ArrayHandle.h:294
vtkm::Id
vtkm::Int32 Id
Represents an ID (index into arrays).
Definition: Types.h:191
vtkm::exec::CellLocatorTwoLevel::PointInsideCell
static VTKM_EXEC vtkm::ErrorCode PointInsideCell(FloatVec3 point, CellShapeTag cellShape, CoordsType cellPoints, FloatVec3 &parametricCoordinates, bool &inside)
Definition: exec/CellLocatorTwoLevel.h:108
VecFromPortalPermute.h
vtkm::exec::CellLocatorTwoLevel::ReadPortal
typename vtkm::cont::ArrayHandle< T >::ReadPortalType ReadPortal
Definition: exec/CellLocatorTwoLevel.h:99
vtkm::cont::CoordinateSystem
Definition: CoordinateSystem.h:25
vtkm::cont::Token
A token to hold the scope of an ArrayHandle or other object.
Definition: Token.h:35
CellInside.h
Math.h
VTKM_CONT
#define VTKM_CONT
Definition: ExportMacros.h:57
vtkm::exec::CellLocatorTwoLevel::CellCount
ReadPortal< vtkm::Id > CellCount
Definition: exec/CellLocatorTwoLevel.h:295
vtkm::ErrorCode::CellNotFound
@ CellNotFound
vtkm::Vec3f
vtkm::Vec< vtkm::FloatDefault, 3 > Vec3f
Vec3f corresponds to a 3-dimensional vector of floating point values.
Definition: Types.h:1014
vtkm::cont::DeviceAdapterId
Definition: DeviceAdapterTag.h:52
vtkm::Vec< DimensionType, 3 >
vtkm::exec::CellLocatorTwoLevel::FindCell
VTKM_EXEC vtkm::ErrorCode FindCell(const FloatVec3 &point, vtkm::Id &cellId, FloatVec3 &parametric, LastCell &lastCell) const
Definition: exec/CellLocatorTwoLevel.h:170
VTKM_RETURN_ON_ERROR
#define VTKM_RETURN_ON_ERROR(call)
Definition: ErrorCode.h:111
vtkm::VecTraits::ComponentType
typename VecType::ComponentType ComponentType
Type of the components in the vector.
Definition: VecTraits.h:73
vtkm::exec::CellLocatorTwoLevel::CellLocatorTwoLevel
VTKM_CONT CellLocatorTwoLevel(const vtkm::internal::cl_uniform_bins::Grid &topLevelGrid, const vtkm::cont::ArrayHandle< DimVec3 > &leafDimensions, const vtkm::cont::ArrayHandle< vtkm::Id > &leafStartIndex, const vtkm::cont::ArrayHandle< vtkm::Id > &cellStartIndex, const vtkm::cont::ArrayHandle< vtkm::Id > &cellCount, const vtkm::cont::ArrayHandle< vtkm::Id > &cellIds, const CellSetType &cellSet, const vtkm::cont::CoordinateSystem &coords, vtkm::cont::DeviceAdapterId device, vtkm::cont::Token &token)
Definition: exec/CellLocatorTwoLevel.h:132
vtkm::exec::CellLocatorTwoLevel::TopLevel
vtkm::internal::cl_uniform_bins::Grid TopLevel
Definition: exec/CellLocatorTwoLevel.h:289
vtkm::TopologyElementTagCell
A tag used to identify the cell elements in a topology.
Definition: TopologyElementTag.h:24
vtkm::exec::CellLocatorTwoLevel::LeafStartIndex
ReadPortal< vtkm::Id > LeafStartIndex
Definition: exec/CellLocatorTwoLevel.h:292
VTKM_ALWAYS_EXPORT
#define VTKM_ALWAYS_EXPORT
Definition: ExportMacros.h:92
vtkm::exec::CellLocatorTwoLevel::CellIds
ReadPortal< vtkm::Id > CellIds
Definition: exec/CellLocatorTwoLevel.h:296
vtkm::exec::CellLocatorTwoLevel::PointInCell
VTKM_EXEC vtkm::ErrorCode PointInCell(const vtkm::Vec3f &point, const vtkm::Id &cid, vtkm::Vec3f &parametric) const
Definition: exec/CellLocatorTwoLevel.h:200
vtkm::exec::CellLocatorTwoLevel::LastCell
Definition: exec/CellLocatorTwoLevel.h:156
ParametricCoordinates.h
vtkm::VecTraits::GetNumberOfComponents
static vtkm::IdComponent GetNumberOfComponents(const VecType &vec)
Number of components in the given vector.
VecTraits.h
TopologyElementTag.h