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 vtkm_exec_celllocatorrectilineargrid_h 11 #define vtkm_exec_celllocatorrectilineargrid_h 12 13 #include <vtkm/Bounds.h> 14 #include <vtkm/TopologyElementTag.h> 15 #include <vtkm/Types.h> 16 #include <vtkm/VecFromPortalPermute.h> 17 18 #include <vtkm/cont/ArrayHandleCartesianProduct.h> 19 #include <vtkm/cont/CellSetStructured.h> 20 21 #include <vtkm/exec/CellInside.h> 22 #include <vtkm/exec/ConnectivityStructured.h> 23 #include <vtkm/exec/ParametricCoordinates.h> 24 25 namespace vtkm 26 { 27 28 namespace exec 29 { 30 31 class VTKM_ALWAYS_EXPORT CellLocatorRectilinearGrid 32 { 33 private: 34 using AxisHandle = vtkm::cont::ArrayHandle<vtkm::FloatDefault>; 35 using RectilinearType = 36 vtkm::cont::ArrayHandleCartesianProduct<AxisHandle, AxisHandle, AxisHandle>; 37 using AxisPortalType = typename AxisHandle::ReadPortalType; 38 using RectilinearPortalType = typename RectilinearType::ReadPortalType; 39 ToId3(vtkm::Id3 && src)40 VTKM_CONT static vtkm::Id3&& ToId3(vtkm::Id3&& src) { return std::move(src); } ToId3(vtkm::Id2 && src)41 VTKM_CONT static vtkm::Id3 ToId3(vtkm::Id2&& src) { return vtkm::Id3(src[0], src[1], 1); } ToId3(vtkm::Id && src)42 VTKM_CONT static vtkm::Id3 ToId3(vtkm::Id&& src) { return vtkm::Id3(src, 1, 1); } 43 44 public: 45 template <vtkm::IdComponent dimensions> CellLocatorRectilinearGrid(const vtkm::Id planeSize,const vtkm::Id rowSize,const vtkm::cont::CellSetStructured<dimensions> & cellSet,const RectilinearType & coords,vtkm::cont::DeviceAdapterId device,vtkm::cont::Token & token)46 VTKM_CONT CellLocatorRectilinearGrid(const vtkm::Id planeSize, 47 const vtkm::Id rowSize, 48 const vtkm::cont::CellSetStructured<dimensions>& cellSet, 49 const RectilinearType& coords, 50 vtkm::cont::DeviceAdapterId device, 51 vtkm::cont::Token& token) 52 : PlaneSize(planeSize) 53 , RowSize(rowSize) 54 , PointDimensions(ToId3(cellSet.GetPointDimensions())) 55 , Dimensions(dimensions) 56 { 57 auto coordsContPortal = coords.ReadPortal(); 58 RectilinearPortalType coordsExecPortal = coords.PrepareForInput(device, token); 59 this->AxisPortals[0] = coordsExecPortal.GetFirstPortal(); 60 this->MinPoint[0] = coordsContPortal.GetFirstPortal().Get(0); 61 this->MaxPoint[0] = coordsContPortal.GetFirstPortal().Get(this->PointDimensions[0] - 1); 62 63 this->AxisPortals[1] = coordsExecPortal.GetSecondPortal(); 64 this->MinPoint[1] = coordsContPortal.GetSecondPortal().Get(0); 65 this->MaxPoint[1] = coordsContPortal.GetSecondPortal().Get(this->PointDimensions[1] - 1); 66 if (dimensions == 3) 67 { 68 this->AxisPortals[2] = coordsExecPortal.GetThirdPortal(); 69 this->MinPoint[2] = coordsContPortal.GetThirdPortal().Get(0); 70 this->MaxPoint[2] = coordsContPortal.GetThirdPortal().Get(this->PointDimensions[2] - 1); 71 } 72 } 73 74 VTKM_EXEC IsInside(const vtkm::Vec3f & point)75 inline bool IsInside(const vtkm::Vec3f& point) const 76 { 77 bool inside = true; 78 if (point[0] < this->MinPoint[0] || point[0] > this->MaxPoint[0]) 79 inside = false; 80 if (point[1] < this->MinPoint[1] || point[1] > this->MaxPoint[1]) 81 inside = false; 82 if (this->Dimensions == 3) 83 { 84 if (point[2] < this->MinPoint[2] || point[2] > this->MaxPoint[2]) 85 inside = false; 86 } 87 return inside; 88 } 89 90 VTKM_EXEC FindCell(const vtkm::Vec3f & point,vtkm::Id & cellId,vtkm::Vec3f & parametric)91 vtkm::ErrorCode FindCell(const vtkm::Vec3f& point, 92 vtkm::Id& cellId, 93 vtkm::Vec3f& parametric) const 94 { 95 if (!this->IsInside(point)) 96 { 97 cellId = -1; 98 return vtkm::ErrorCode::CellNotFound; 99 } 100 101 // Get the Cell Id from the point. 102 vtkm::Id3 logicalCell(0, 0, 0); 103 for (vtkm::Int32 dim = 0; dim < this->Dimensions; ++dim) 104 { 105 // 106 // When searching for points, we consider the max value of the cell 107 // to be apart of the next cell. If the point falls on the boundary of the 108 // data set, then it is technically inside a cell. This checks for that case 109 // 110 if (point[dim] == MaxPoint[dim]) 111 { 112 logicalCell[dim] = this->PointDimensions[dim] - 2; 113 parametric[dim] = static_cast<vtkm::FloatDefault>(1); 114 continue; 115 } 116 117 vtkm::Id minIndex = 0; 118 vtkm::Id maxIndex = this->PointDimensions[dim] - 1; 119 vtkm::FloatDefault minVal; 120 vtkm::FloatDefault maxVal; 121 minVal = this->AxisPortals[dim].Get(minIndex); 122 maxVal = this->AxisPortals[dim].Get(maxIndex); 123 while (maxIndex > minIndex + 1) 124 { 125 vtkm::Id midIndex = (minIndex + maxIndex) / 2; 126 vtkm::FloatDefault midVal = this->AxisPortals[dim].Get(midIndex); 127 if (point[dim] <= midVal) 128 { 129 maxIndex = midIndex; 130 maxVal = midVal; 131 } 132 else 133 { 134 minIndex = midIndex; 135 minVal = midVal; 136 } 137 } 138 logicalCell[dim] = minIndex; 139 parametric[dim] = (point[dim] - minVal) / (maxVal - minVal); 140 } 141 // Get the actual cellId, from the logical cell index of the cell 142 cellId = logicalCell[2] * this->PlaneSize + logicalCell[1] * this->RowSize + logicalCell[0]; 143 144 return vtkm::ErrorCode::Success; 145 } 146 147 VTKM_DEPRECATED(1.6, "Locators are no longer pointers. Use . operator.") 148 VTKM_EXEC CellLocatorRectilinearGrid* operator->() { return this; } 149 VTKM_DEPRECATED(1.6, "Locators are no longer pointers. Use . operator.") 150 VTKM_EXEC const CellLocatorRectilinearGrid* operator->() const { return this; } 151 152 private: 153 vtkm::Id PlaneSize; 154 vtkm::Id RowSize; 155 156 AxisPortalType AxisPortals[3]; 157 vtkm::Id3 PointDimensions; 158 vtkm::Vec3f MinPoint; 159 vtkm::Vec3f MaxPoint; 160 vtkm::Id Dimensions; 161 }; 162 } //namespace exec 163 } //namespace vtkm 164 165 #endif //vtkm_exec_celllocatorrectilineargrid_h 166