diff --git a/vtkm/cont/CellLocatorUniformGrid.cxx b/vtkm/cont/CellLocatorUniformGrid.cxx index 44b0bcc55..22ef55e01 100644 --- a/vtkm/cont/CellLocatorUniformGrid.cxx +++ b/vtkm/cont/CellLocatorUniformGrid.cxx @@ -38,32 +38,38 @@ void CellLocatorUniformGrid::Build() if (cellSet.IsSameType(Structured2DType())) { this->Is3D = false; - this->Bounds = coords.GetBounds(); - vtkm::Id2 cellDims = - cellSet.Cast().GetSchedulingRange(vtkm::TopologyElementTagCell()); - this->CellDims = vtkm::Id3(cellDims[0], cellDims[1], 0); - this->RangeTransform[0] = static_cast(this->CellDims[0]) / - static_cast(this->Bounds.X.Length()); - this->RangeTransform[1] = static_cast(this->CellDims[1]) / - static_cast(this->Bounds.Y.Length()); + Structured2DType structuredCellSet = cellSet.Cast(); + vtkm::Id2 pointDims = structuredCellSet.GetSchedulingRange(vtkm::TopologyElementTagPoint()); + this->PointDims = vtkm::Id3(pointDims[0], pointDims[1], 1); } else if (cellSet.IsSameType(Structured3DType())) { this->Is3D = true; - this->Bounds = coords.GetBounds(); - this->CellDims = - cellSet.Cast().GetSchedulingRange(vtkm::TopologyElementTagCell()); - this->RangeTransform[0] = static_cast(this->CellDims[0]) / - static_cast(this->Bounds.X.Length()); - this->RangeTransform[1] = static_cast(this->CellDims[1]) / - static_cast(this->Bounds.Y.Length()); - this->RangeTransform[2] = static_cast(this->CellDims[2]) / - static_cast(this->Bounds.Z.Length()); + Structured3DType structuredCellSet = cellSet.Cast(); + this->PointDims = structuredCellSet.GetSchedulingRange(vtkm::TopologyElementTagPoint()); } else { throw vtkm::cont::ErrorInternal("Cells are not structured."); } + + UniformType uniformCoords = coords.GetData().Cast(); + this->Origin = uniformCoords.GetPortalConstControl().GetOrigin(); + + vtkm::Vec spacing = uniformCoords.GetPortalConstControl().GetSpacing(); + vtkm::Vec unitLength; + unitLength[0] = static_cast(PointDims[0] - 1); + unitLength[1] = static_cast(PointDims[1] - 1); + unitLength[2] = static_cast(PointDims[2] - 1); + + this->MaxPoint = Origin + spacing * unitLength; + this->InvSpacing[0] = 1.f / spacing[0]; + this->InvSpacing[1] = 1.f / spacing[1]; + this->InvSpacing[2] = 1.f / spacing[2]; + + this->CellDims[0] = PointDims[0] - 1; + this->CellDims[1] = PointDims[1] - 1; + this->CellDims[2] = PointDims[2] - 1; } namespace @@ -93,10 +99,11 @@ const vtkm::exec::CellLocator* CellLocatorUniformGrid::PrepareForExecution( success = vtkm::cont::TryExecuteOnDevice(device, CellLocatorUniformGridPrepareForExecutionFunctor<3>(), this->ExecutionObjectHandle, - this->Bounds, - this->RangeTransform, this->CellDims, - this->GetCellSet().template Cast(), + this->PointDims, + this->Origin, + this->InvSpacing, + this->MaxPoint, this->GetCoordinates().GetData()); } else @@ -104,10 +111,11 @@ const vtkm::exec::CellLocator* CellLocatorUniformGrid::PrepareForExecution( success = vtkm::cont::TryExecuteOnDevice(device, CellLocatorUniformGridPrepareForExecutionFunctor<2>(), this->ExecutionObjectHandle, - this->Bounds, - this->RangeTransform, this->CellDims, - this->GetCellSet().template Cast(), + this->PointDims, + this->Origin, + this->InvSpacing, + this->MaxPoint, this->GetCoordinates().GetData()); } if (!success) diff --git a/vtkm/cont/CellLocatorUniformGrid.h b/vtkm/cont/CellLocatorUniformGrid.h index f24a3360b..3b7eb8a40 100644 --- a/vtkm/cont/CellLocatorUniformGrid.h +++ b/vtkm/cont/CellLocatorUniformGrid.h @@ -32,9 +32,11 @@ protected: VTKM_CONT void Build() override; private: - vtkm::Bounds Bounds; - vtkm::Vec3f RangeTransform; vtkm::Id3 CellDims; + vtkm::Id3 PointDims; + vtkm::Vec Origin; + vtkm::Vec InvSpacing; + vtkm::Vec MaxPoint; bool Is3D = true; mutable vtkm::cont::VirtualObjectHandle ExecutionObjectHandle; diff --git a/vtkm/exec/CellLocatorUniformGrid.h b/vtkm/exec/CellLocatorUniformGrid.h index 6a2e95a83..2cbc0514c 100644 --- a/vtkm/exec/CellLocatorUniformGrid.h +++ b/vtkm/exec/CellLocatorUniformGrid.h @@ -41,18 +41,18 @@ private: public: VTKM_CONT - CellLocatorUniformGrid(const vtkm::Bounds& bounds, - const vtkm::Vec3f rangeTransform, - const vtkm::Id3 cellDims, - const vtkm::cont::CellSetStructured& cellSet, + CellLocatorUniformGrid(const vtkm::Id3 cellDims, + const vtkm::Id3 pointDims, + const vtkm::Vec origin, + const vtkm::Vec invSpacing, + const vtkm::Vec maxPoint, const vtkm::cont::ArrayHandleVirtualCoordinates& coords, DeviceAdapter) - : Bounds(bounds) - , RangeTransform(rangeTransform) - , CellDims(cellDims) - , PlaneSize(cellDims[0] * cellDims[1]) - , RowSize(cellDims[0]) - , CellSet(cellSet.PrepareForInput(DeviceAdapter(), FromType(), ToType())) + : CellDims(cellDims) + , PointDims(pointDims) + , Origin(origin) + , InvSpacing(invSpacing) + , MaxPoint(maxPoint) , Coords(coords.PrepareForInput(DeviceAdapter())) { } @@ -63,51 +63,65 @@ public: // troublesome with CUDA __host__ __device__ markup. } + template + VTKM_EXEC inline bool IsInside(const vtkm::Vec& point) const + { + bool inside = true; + if (point[0] < this->Origin[0] || point[0] > this->MaxPoint[0]) + inside = false; + if (point[1] < this->Origin[1] || point[1] > this->MaxPoint[1]) + inside = false; + if (point[2] < this->Origin[2] || point[2] > this->MaxPoint[2]) + inside = false; + return inside; + } + VTKM_EXEC void FindCell(const vtkm::Vec3f& point, vtkm::Id& cellId, vtkm::Vec3f& parametric, const vtkm::exec::FunctorBase& worklet) const override { - if (!Bounds.Contains(point)) + (void)worklet; //suppress unused warning + if (!IsInside(point)) { cellId = -1; return; } // Get the Cell Id from the point. vtkm::Id3 logicalCell(0, 0, 0); - logicalCell[0] = (point[0] == Bounds.X.Max) - ? CellDims[0] - 1 - : static_cast(vtkm::Floor((point[0] - Bounds.X.Min) * RangeTransform[0])); - logicalCell[1] = (point[1] == Bounds.Y.Max) - ? CellDims[1] - 1 - : static_cast(vtkm::Floor((point[1] - Bounds.Y.Min) * RangeTransform[1])); - if (dimensions == 3) - { - logicalCell[2] = (point[2] == Bounds.Z.Max) - ? CellDims[2] - 1 - : static_cast(vtkm::Floor((point[2] - Bounds.Z.Min) * RangeTransform[2])); - } - // Get the actual cellId, from the logical cell index of the cell - cellId = logicalCell[2] * PlaneSize + logicalCell[1] * RowSize + logicalCell[0]; - bool success = false; - using IndicesType = typename CellSetPortal::IndicesType; - IndicesType cellPointIndices = CellSet.GetIndices(cellId); - vtkm::VecFromPortalPermute cellPoints(&cellPointIndices, Coords); - auto cellShape = CellSet.GetCellShape(cellId); - // Get Parametric Coordinates from the cell, for the point. - parametric = vtkm::exec::WorldCoordinatesToParametricCoordinates( - cellPoints, point, cellShape, success, worklet); + vtkm::Vec temp; + temp = point - Origin; + temp = temp * InvSpacing; + + //make sure that if we border the upper edge, we sample the correct cell + logicalCell = temp; + if (logicalCell[0] == this->CellDims[0]) + { + logicalCell[0]--; + } + if (logicalCell[1] == this->CellDims[1]) + { + logicalCell[1]--; + } + if (logicalCell[2] == this->CellDims[2]) + { + logicalCell[2]--; + } + if (dimensions == 2) + logicalCell[2] = 0; + cellId = + (logicalCell[2] * this->CellDims[1] + logicalCell[1]) * this->CellDims[0] + logicalCell[0]; + parametric = temp - logicalCell; } private: - vtkm::Bounds Bounds; - vtkm::Vec3f RangeTransform; vtkm::Id3 CellDims; - vtkm::Id PlaneSize; - vtkm::Id RowSize; - CellSetPortal CellSet; + vtkm::Id3 PointDims; + vtkm::Vec Origin; + vtkm::Vec InvSpacing; + vtkm::Vec MaxPoint; CoordsPortal Coords; }; }