diff --git a/vtkm/cont/CMakeLists.txt b/vtkm/cont/CMakeLists.txt index e3c52f99f..c827a03d3 100644 --- a/vtkm/cont/CMakeLists.txt +++ b/vtkm/cont/CMakeLists.txt @@ -55,7 +55,7 @@ set(headers CellLocatorBoundingIntervalHierarchy.h CellLocatorHelper.h CellLocatorRectilinearGrid.h - CellLocatorTwoLevelUniformGrid.h + CellLocatorUniformBins.h CellLocatorUniformGrid.h CellSet.h CellSetExplicit.h @@ -156,6 +156,7 @@ set(device_sources BoundsGlobalCompute.cxx CellLocatorBoundingIntervalHierarchy.cxx CellLocatorRectilinearGrid.cxx + CellLocatorUniformBins.cxx CellLocatorUniformGrid.cxx CellSet.cxx CellSetExplicit.cxx diff --git a/vtkm/cont/CellLocatorHelper.h b/vtkm/cont/CellLocatorHelper.h index aa5c0a80e..170d15ac4 100644 --- a/vtkm/cont/CellLocatorHelper.h +++ b/vtkm/cont/CellLocatorHelper.h @@ -20,7 +20,7 @@ #ifndef vtk_m_cont_CellLocatorHelper_h #define vtk_m_cont_CellLocatorHelper_h -#include +#include #include namespace vtkm @@ -53,11 +53,8 @@ public: /// Builds the cell locator lookup structure /// - template - void Build(CellSetList cellSetTypes = CellSetList()) + void Build() { - VTKM_IS_LIST_TAG(CellSetList); - if (IsUniformGrid(this->CellSet, this->Coordinates)) { // nothing to build for uniform grid @@ -66,12 +63,31 @@ public: { this->Locator.SetCellSet(this->CellSet); this->Locator.SetCoordinates(this->Coordinates); - this->Locator.Build(cellSetTypes); + this->Locator.Update(); } } class FindCellWorklet : public vtkm::worklet::WorkletMapField { + public: + using ControlSignature = void(FieldIn points, + ExecObject locator, + FieldOut cellIds, + FieldOut pcoords); + using ExecutionSignature = void(_1, _2, _3, _4); + + template + VTKM_EXEC void operator()(const vtkm::Vec& point, + const LocatorType& locator, + vtkm::Id& cellId, + vtkm::Vec& pcoords) const + { + locator->FindCell(point, cellId, pcoords, *this); + } + }; + + class FindCellWorkletUG : public vtkm::worklet::WorkletMapField + { private: template using ConnectivityType = vtkm::exec::ConnectivityStructured>& parametricCoords, CellSetList cellSetTypes = CellSetList()) const { + (void)cellSetTypes; // unused for now + if (IsUniformGrid(this->CellSet, this->Coordinates)) { auto coordinates = this->Coordinates.GetData().Cast(); auto cellset = this->CellSet.ResetCellSetList(StructuredCellSetList()); - vtkm::worklet::DispatcherMapField dispatcher; + vtkm::worklet::DispatcherMapField dispatcher; dispatcher.Invoke(points, cellset, coordinates, cellIds, parametricCoords); } else { - this->Locator.FindCells(points, cellIds, parametricCoords, cellSetTypes); + vtkm::worklet::DispatcherMapField dispatcher; + dispatcher.Invoke(points, this->Locator, cellIds, parametricCoords); } } private: vtkm::cont::DynamicCellSet CellSet; vtkm::cont::CoordinateSystem Coordinates; - vtkm::cont::CellLocatorTwoLevelUniformGrid Locator; + vtkm::cont::CellLocatorUniformBins Locator; }; } } // vtkm::cont diff --git a/vtkm/cont/CellLocatorTwoLevelUniformGrid.h b/vtkm/cont/CellLocatorTwoLevelUniformGrid.h deleted file mode 100644 index e6d641126..000000000 --- a/vtkm/cont/CellLocatorTwoLevelUniformGrid.h +++ /dev/null @@ -1,724 +0,0 @@ -//============================================================================ -// Copyright (c) Kitware, Inc. -// All rights reserved. -// See LICENSE.txt for details. -// 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. -// -// Copyright 2017 National Technology & Engineering Solutions of Sandia, LLC (NTESS). -// Copyright 2017 UT-Battelle, LLC. -// Copyright 2017 Los Alamos National Security. -// -// Under the terms of Contract DE-NA0003525 with NTESS, -// the U.S. Government retains certain rights in this software. -// -// Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National -// Laboratory (LANL), the U.S. Government retains certain rights in -// this software. -//============================================================================ -#ifndef vtk_m_cont_CellLocatorTwoLevelUniformGrid_h -#define vtk_m_cont_CellLocatorTwoLevelUniformGrid_h - -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include - -#include -#include -#include - -#include -#include -#include -#include - -namespace vtkm -{ -namespace cont -{ - -class CellLocatorTwoLevelUniformGrid -{ -public: - CellLocatorTwoLevelUniformGrid() - : DensityL1(32.0f) - , DensityL2(2.0f) - { - } - - /// Get/Set the desired approximate number of cells per level 1 bin - /// - void SetDensityL1(vtkm::FloatDefault val) { this->DensityL1 = val; } - vtkm::FloatDefault GetDensityL1() const { return this->DensityL1; } - - /// Get/Set the desired approximate number of cells per level 1 bin - /// - void SetDensityL2(vtkm::FloatDefault val) { this->DensityL2 = val; } - vtkm::FloatDefault GetDensityL2() const { return this->DensityL2; } - - void SetCellSet(const vtkm::cont::DynamicCellSet& cellset) { this->CellSet = cellset; } - - const vtkm::cont::DynamicCellSet& GetCellSet() const { return this->CellSet; } - - void SetCoordinates(const vtkm::cont::CoordinateSystem& coords) { this->Coordinates = coords; } - - const vtkm::cont::CoordinateSystem& GetCoordinates() const { return this->Coordinates; } - - void PrintSummary(std::ostream& out) const - { - out << "DensityL1: " << this->DensityL1 << "\n"; - out << "DensityL2: " << this->DensityL2 << "\n"; - out << "Input CellSet: \n"; - this->CellSet.PrintSummary(out); - out << "Input Coordinates: \n"; - this->Coordinates.PrintSummary(out); - out << "LookupStructure:\n"; - out << " TopLevelGrid\n"; - out << " Dimensions: " << this->LookupStructure.TopLevel.Dimensions << "\n"; - out << " Origin: " << this->LookupStructure.TopLevel.Origin << "\n"; - out << " BinSize: " << this->LookupStructure.TopLevel.BinSize << "\n"; - out << " LeafDimensions:\n"; - vtkm::cont::printSummary_ArrayHandle(this->LookupStructure.LeafDimensions, out); - out << " LeafStartIndex:\n"; - vtkm::cont::printSummary_ArrayHandle(this->LookupStructure.LeafStartIndex, out); - out << " CellStartIndex:\n"; - vtkm::cont::printSummary_ArrayHandle(this->LookupStructure.CellStartIndex, out); - out << " CellCount:\n"; - vtkm::cont::printSummary_ArrayHandle(this->LookupStructure.CellCount, out); - out << " CellIds:\n"; - vtkm::cont::printSummary_ArrayHandle(this->LookupStructure.CellIds, out); - } - -private: - using DimensionType = vtkm::Int16; - using DimVec3 = vtkm::Vec; - using FloatVec3 = vtkm::Vec; - - struct BinsBBox - { - DimVec3 Min; - DimVec3 Max; - - VTKM_EXEC - bool Empty() const - { - return (this->Max[0] < this->Min[0]) || (this->Max[1] < this->Min[1]) || - (this->Max[2] < this->Min[2]); - } - }; - - struct Bounds - { - FloatVec3 Min; - FloatVec3 Max; - }; - - VTKM_EXEC_CONT static DimVec3 ComputeGridDimension(vtkm::Id numberOfCells, - const FloatVec3& size, - vtkm::FloatDefault density) - { - vtkm::FloatDefault nsides = 0.0f; - vtkm::FloatDefault volume = 1.0f; - vtkm::FloatDefault maxside = vtkm::Max(size[0], vtkm::Max(size[1], size[2])); - for (int i = 0; i < 3; ++i) - { - if (size[i] / maxside >= 1e-4f) - { - nsides += 1.0f; - volume *= size[i]; - } - } - - auto r = vtkm::Pow((static_cast(numberOfCells) / (volume * density)), - 1.0f / nsides); - return vtkm::Max(DimVec3(1), - DimVec3(static_cast(size[0] * r), - static_cast(size[1] * r), - static_cast(size[2] * r))); - } - - VTKM_EXEC static vtkm::Id ComputeFlatIndex(const DimVec3& idx, const DimVec3 dim) - { - return idx[0] + (dim[0] * (idx[1] + (dim[1] * idx[2]))); - } - - VTKM_EXEC static vtkm::exec::twolevelgrid::Grid ComputeLeafGrid( - const DimVec3& idx, - const DimVec3& dim, - const vtkm::exec::twolevelgrid::Grid& l1Grid) - { - return { dim, - l1Grid.Origin + (static_cast(idx) * l1Grid.BinSize), - l1Grid.BinSize / static_cast(dim) }; - } - - template - VTKM_EXEC static Bounds ComputeCellBounds(const PointsVecType& points) - { - using CoordsType = typename vtkm::VecTraits::ComponentType; - auto numPoints = vtkm::VecTraits::GetNumberOfComponents(points); - - CoordsType minp = points[0], maxp = points[0]; - for (vtkm::IdComponent i = 1; i < numPoints; ++i) - { - minp = vtkm::Min(minp, points[i]); - maxp = vtkm::Max(maxp, points[i]); - } - - return { FloatVec3(minp), FloatVec3(maxp) }; - } - - VTKM_EXEC static BinsBBox ComputeIntersectingBins(const Bounds cellBounds, - const vtkm::exec::twolevelgrid::Grid& grid) - { - auto minb = static_cast((cellBounds.Min - grid.Origin) / grid.BinSize); - auto maxb = static_cast((cellBounds.Max - grid.Origin) / grid.BinSize); - - return { vtkm::Max(DimVec3(0), minb), vtkm::Min(grid.Dimensions - DimVec3(1), maxb) }; - } - - VTKM_EXEC static vtkm::Id GetNumberOfBins(const BinsBBox& binsBBox) - { - return binsBBox.Empty() ? 0 : ((binsBBox.Max[0] - binsBBox.Min[0] + 1) * - (binsBBox.Max[1] - binsBBox.Min[1] + 1) * - (binsBBox.Max[2] - binsBBox.Min[2] + 1)); - } - - class BBoxIterator - { - public: - VTKM_EXEC_CONT BBoxIterator(const BinsBBox& bbox, const DimVec3& dim) - : BBox(bbox) - , Dim(dim) - , Idx(bbox.Min) - , StepY(dim[0] - (bbox.Max[0] - bbox.Min[0] + 1)) - , StepZ((dim[0] * dim[1]) - ((bbox.Max[1] - bbox.Min[1] + 1) * dim[0])) - , FlatIdx(ComputeFlatIndex(this->Idx, dim)) - , DoneFlag(bbox.Empty()) - { - } - - VTKM_EXEC_CONT void Init() - { - this->Idx = this->BBox.Min; - this->FlatIdx = ComputeFlatIndex(this->Idx, this->Dim); - this->DoneFlag = this->BBox.Empty(); - } - - VTKM_EXEC_CONT bool Done() const { return this->DoneFlag; } - - VTKM_EXEC_CONT void Next() - { - if (!this->DoneFlag) - { - ++this->Idx[0]; - this->FlatIdx += 1; - if (this->Idx[0] > this->BBox.Max[0]) - { - this->Idx[0] = this->BBox.Min[0]; - ++this->Idx[1]; - this->FlatIdx += this->StepY; - if (this->Idx[1] > this->BBox.Max[1]) - { - this->Idx[1] = this->BBox.Min[1]; - ++this->Idx[2]; - this->FlatIdx += this->StepZ; - if (this->Idx[2] > this->BBox.Max[2]) - { - this->DoneFlag = true; - } - } - } - } - } - - VTKM_EXEC_CONT const DimVec3& GetIdx() const { return this->Idx; } - - VTKM_EXEC_CONT vtkm::Id GetFlatIdx() const { return this->FlatIdx; } - - private: - BinsBBox BBox; - DimVec3 Dim; - DimVec3 Idx; - vtkm::Id StepY, StepZ; - vtkm::Id FlatIdx; - bool DoneFlag; - }; - - // TODO: This function may return false positives for non 3D cells as the - // tests are done on the projection of the point on the cell. Extra checks - // should be added to test if the point actually falls on the cell. - template - VTKM_EXEC static bool PointInsideCell(FloatVec3 point, - CellShapeTag cellShape, - CoordsType cellPoints, - const vtkm::exec::FunctorBase& worklet, - FloatVec3& parametricCoordinates) - { - auto bounds = ComputeCellBounds(cellPoints); - if (point[0] >= bounds.Min[0] && point[0] <= bounds.Max[0] && point[1] >= bounds.Min[1] && - point[1] <= bounds.Max[1] && point[2] >= bounds.Min[2] && point[2] <= bounds.Max[2]) - { - bool success = false; - parametricCoordinates = vtkm::exec::WorldCoordinatesToParametricCoordinates( - cellPoints, point, cellShape, success, worklet); - return success && vtkm::exec::CellInside(parametricCoordinates, cellShape); - } - return false; - } - -public: - class CountBinsL1 : public vtkm::worklet::WorkletMapPointToCell - { - public: - using ControlSignature = void(CellSetIn cellset, FieldInPoint coords, FieldOutCell bincount); - using ExecutionSignature = void(_2, _3); - - CountBinsL1(const vtkm::exec::twolevelgrid::Grid& grid) - : L1Grid(grid) - { - } - - template - VTKM_EXEC void operator()(const PointsVecType& points, vtkm::Id& numBins) const - { - auto cellBounds = ComputeCellBounds(points); - auto binsBBox = ComputeIntersectingBins(cellBounds, this->L1Grid); - numBins = GetNumberOfBins(binsBBox); - } - - private: - vtkm::exec::twolevelgrid::Grid L1Grid; - }; - - class FindBinsL1 : public vtkm::worklet::WorkletMapPointToCell - { - public: - using ControlSignature = void(CellSetIn cellset, - FieldInPoint coords, - FieldInCell offsets, - WholeArrayOut binIds); - using ExecutionSignature = void(_2, _3, _4); - - FindBinsL1(const vtkm::exec::twolevelgrid::Grid& grid) - : L1Grid(grid) - { - } - - template - VTKM_EXEC void operator()(const PointsVecType& points, - vtkm::Id offset, - BinIdsPortalType& binIds) const - { - auto cellBounds = ComputeCellBounds(points); - auto binsBBox = ComputeIntersectingBins(cellBounds, this->L1Grid); - - for (BBoxIterator i(binsBBox, this->L1Grid.Dimensions); !i.Done(); i.Next()) - { - binIds.Set(offset, i.GetFlatIdx()); - ++offset; - } - } - - private: - vtkm::exec::twolevelgrid::Grid L1Grid; - }; - - class GenerateBinsL1 : public vtkm::worklet::WorkletMapField - { - public: - using ControlSignature = void(FieldIn binIds, FieldIn cellCounts, WholeArrayOut dimensions); - using ExecutionSignature = void(_1, _2, _3); - - using InputDomain = _1; - - GenerateBinsL1(FloatVec3 size, vtkm::FloatDefault density) - : Size(size) - , Density(density) - { - } - - template - VTKM_EXEC void operator()(vtkm::Id binId, - vtkm::Id numCells, - OutputDimensionsPortal& dimensions) const - { - dimensions.Set(binId, ComputeGridDimension(numCells, this->Size, this->Density)); - } - - private: - FloatVec3 Size; - vtkm::FloatDefault Density; - }; - - class CountBinsL2 : public vtkm::worklet::WorkletMapPointToCell - { - public: - using ControlSignature = void(CellSetIn cellset, - FieldInPoint coords, - WholeArrayIn binDimensions, - FieldOutCell bincount); - using ExecutionSignature = void(_2, _3, _4); - - CountBinsL2(const vtkm::exec::twolevelgrid::Grid& grid) - : L1Grid(grid) - { - } - - template - VTKM_EXEC void operator()(const PointsVecType& points, - const BinDimensionsPortalType& binDimensions, - vtkm::Id& numBins) const - { - auto cellBounds = ComputeCellBounds(points); - auto binsBBox = ComputeIntersectingBins(cellBounds, this->L1Grid); - - numBins = 0; - for (BBoxIterator i(binsBBox, this->L1Grid.Dimensions); !i.Done(); i.Next()) - { - vtkm::exec::twolevelgrid::Grid leaf = - ComputeLeafGrid(i.GetIdx(), binDimensions.Get(i.GetFlatIdx()), this->L1Grid); - auto binsBBoxL2 = ComputeIntersectingBins(cellBounds, leaf); - numBins += GetNumberOfBins(binsBBoxL2); - } - } - - private: - vtkm::exec::twolevelgrid::Grid L1Grid; - }; - - class FindBinsL2 : public vtkm::worklet::WorkletMapPointToCell - { - public: - using ControlSignature = void(CellSetIn cellset, - FieldInPoint coords, - WholeArrayIn binDimensions, - WholeArrayIn binStarts, - FieldInCell offsets, - WholeArrayOut binIds, - WholeArrayOut cellIds); - using ExecutionSignature = void(InputIndex, _2, _3, _4, _5, _6, _7); - - FindBinsL2(const vtkm::exec::twolevelgrid::Grid& grid) - : L1Grid(grid) - { - } - - template - VTKM_EXEC void operator()(vtkm::Id cellId, - const PointsVecType& points, - const BinDimensionsPortalType& binDimensions, - const BinStartsPortalType& binStarts, - vtkm::Id offset, - BinIdsPortalType binIds, - CellIdsPortalType cellIds) const - { - auto cellBounds = ComputeCellBounds(points); - auto binsBBox = ComputeIntersectingBins(cellBounds, this->L1Grid); - - for (BBoxIterator i(binsBBox, this->L1Grid.Dimensions); !i.Done(); i.Next()) - { - vtkm::exec::twolevelgrid::Grid leaf = - ComputeLeafGrid(i.GetIdx(), binDimensions.Get(i.GetFlatIdx()), this->L1Grid); - auto binsBBoxL2 = ComputeIntersectingBins(cellBounds, leaf); - vtkm::Id leafStart = binStarts.Get(i.GetFlatIdx()); - - for (BBoxIterator j(binsBBoxL2, leaf.Dimensions); !j.Done(); j.Next()) - { - binIds.Set(offset, leafStart + j.GetFlatIdx()); - cellIds.Set(offset, cellId); - ++offset; - } - } - } - - private: - vtkm::exec::twolevelgrid::Grid L1Grid; - }; - - class GenerateBinsL2 : public vtkm::worklet::WorkletMapField - { - public: - using ControlSignature = void(FieldIn binIds, - FieldIn startsIn, - FieldIn countsIn, - WholeArrayOut startsOut, - WholeArrayOut countsOut); - using ExecutionSignature = void(_1, _2, _3, _4, _5); - - using InputDomain = _1; - - template - VTKM_EXEC void operator()(vtkm::Id binIndex, - vtkm::Id start, - vtkm::Id count, - CellStartsPortalType& cellStarts, - CellCountsPortalType& cellCounts) const - { - cellStarts.Set(binIndex, start); - cellCounts.Set(binIndex, count); - } - }; - - struct DimensionsToCount - { - VTKM_EXEC vtkm::Id operator()(const DimVec3& dim) const { return dim[0] * dim[1] * dim[2]; } - }; - - /// Builds the cell locator lookup structure - /// - template - void Build(CellSetList cellSetTypes = CellSetList()) - { - VTKM_IS_LIST_TAG(CellSetList); - - vtkm::worklet::Invoker invoke; - - auto cellset = this->CellSet.ResetCellSetList(cellSetTypes); - const auto& coords = this->Coordinates; - TwoLevelUniformGrid ls; - - // 1: Compute the top level grid - auto bounds = this->Coordinates.GetBounds(); - FloatVec3 bmin(static_cast(bounds.X.Min), - static_cast(bounds.Y.Min), - static_cast(bounds.Z.Min)); - FloatVec3 bmax(static_cast(bounds.X.Max), - static_cast(bounds.Y.Max), - static_cast(bounds.Z.Max)); - auto size = bmax - bmin; - auto fudge = vtkm::Max(FloatVec3(1e-6f), size * 1e-4f); - size += 2.0f * fudge; - - ls.TopLevel.Dimensions = - ComputeGridDimension(cellset.GetNumberOfCells(), size, this->DensityL1); - ls.TopLevel.Origin = bmin - fudge; - ls.TopLevel.BinSize = size / static_cast(ls.TopLevel.Dimensions); - - // 2: For each cell, find the number of top level bins they intersect - vtkm::cont::ArrayHandle binCounts; - CountBinsL1 countL1(ls.TopLevel); - invoke(countL1, cellset, coords, binCounts); - - // 3: Total number of unique (cell, bin) pairs (for pre-allocating arrays) - vtkm::Id numPairsL1 = vtkm::cont::Algorithm::ScanExclusive(binCounts, binCounts); - - // 4: For each cell find the top level bins that intersect it - vtkm::cont::ArrayHandle binIds; - binIds.Allocate(numPairsL1); - FindBinsL1 findL1(ls.TopLevel); - invoke(findL1, cellset, coords, binCounts, binIds); - binCounts.ReleaseResources(); - - // 5: From above, find the number of cells that intersect each top level bin - vtkm::cont::Algorithm::Sort(binIds); - vtkm::cont::ArrayHandle bins; - vtkm::cont::ArrayHandle cellsPerBin; - vtkm::cont::Algorithm::ReduceByKey( - binIds, - vtkm::cont::make_ArrayHandleConstant(vtkm::Id(1), numPairsL1), - bins, - cellsPerBin, - vtkm::Sum()); - binIds.ReleaseResources(); - - // 6: Compute level-2 dimensions - vtkm::Id numberOfBins = - ls.TopLevel.Dimensions[0] * ls.TopLevel.Dimensions[1] * ls.TopLevel.Dimensions[2]; - vtkm::cont::ArrayCopy(vtkm::cont::make_ArrayHandleConstant(DimVec3(0), numberOfBins), - ls.LeafDimensions); - GenerateBinsL1 generateL1(ls.TopLevel.BinSize, this->DensityL2); - invoke(generateL1, bins, cellsPerBin, ls.LeafDimensions); - bins.ReleaseResources(); - cellsPerBin.ReleaseResources(); - - // 7: Compute number of level-2 bins - vtkm::Id numberOfLeaves = vtkm::cont::Algorithm::ScanExclusive( - vtkm::cont::make_ArrayHandleTransform(ls.LeafDimensions, DimensionsToCount()), - ls.LeafStartIndex); - - - // 8: For each cell, find the number of l2 bins they intersect - CountBinsL2 countL2(ls.TopLevel); - invoke(countL2, cellset, coords, ls.LeafDimensions, binCounts); - - // 9: Total number of unique (cell, bin) pairs (for pre-allocating arrays) - vtkm::Id numPairsL2 = vtkm::cont::Algorithm::ScanExclusive(binCounts, binCounts); - - // 10: For each cell, find the l2 bins they intersect - binIds.Allocate(numPairsL2); - ls.CellIds.Allocate(numPairsL2); - FindBinsL2 findL2(ls.TopLevel); - invoke( - findL2, cellset, coords, ls.LeafDimensions, ls.LeafStartIndex, binCounts, binIds, ls.CellIds); - binCounts.ReleaseResources(); - - // 11: From above, find the cells that each l2 bin intersects - vtkm::cont::Algorithm::SortByKey(binIds, ls.CellIds); - vtkm::cont::Algorithm::ReduceByKey( - binIds, - vtkm::cont::make_ArrayHandleConstant(vtkm::Id(1), numPairsL2), - bins, - cellsPerBin, - vtkm::Sum()); - binIds.ReleaseResources(); - - // 12: Generate the leaf bin arrays - vtkm::cont::ArrayHandle cellsStart; - vtkm::cont::Algorithm::ScanExclusive(cellsPerBin, cellsStart); - - vtkm::cont::ArrayCopy(vtkm::cont::ArrayHandleConstant(0, numberOfLeaves), - ls.CellStartIndex); - vtkm::cont::ArrayCopy(vtkm::cont::ArrayHandleConstant(0, numberOfLeaves), - ls.CellCount); - invoke(GenerateBinsL2{}, bins, cellsStart, cellsPerBin, ls.CellStartIndex, ls.CellCount); - - std::swap(this->LookupStructure, ls); - } - - struct TwoLevelUniformGrid : public vtkm::cont::ExecutionObjectBase - { - template - VTKM_CONT vtkm::exec::twolevelgrid::TwoLevelUniformGridExecutionObject - PrepareForExecution(DeviceAdapter device) const - { - vtkm::exec::twolevelgrid::TwoLevelUniformGridExecutionObject deviceObject; - deviceObject.TopLevel = this->TopLevel; - deviceObject.LeafDimensions = this->LeafDimensions.PrepareForInput(device); - deviceObject.LeafStartIndex = this->LeafStartIndex.PrepareForInput(device); - deviceObject.CellStartIndex = this->CellStartIndex.PrepareForInput(device); - deviceObject.CellCount = this->CellCount.PrepareForInput(device); - deviceObject.CellIds = this->CellIds.PrepareForInput(device); - return deviceObject; - } - - vtkm::exec::twolevelgrid::Grid TopLevel; - - vtkm::cont::ArrayHandle LeafDimensions; - vtkm::cont::ArrayHandle LeafStartIndex; - - vtkm::cont::ArrayHandle CellStartIndex; - vtkm::cont::ArrayHandle CellCount; - vtkm::cont::ArrayHandle CellIds; - }; - - class FindCellWorklet : public vtkm::worklet::WorkletMapField - { - public: - using ControlSignature = void(FieldIn points, - WholeCellSetIn<> cellSet, - WholeArrayIn coordinates, - ExecObject lookupStruct, - FieldOut cellIds, - FieldOut parametricCoordinates); - using ExecutionSignature = void(_1, _2, _3, _4, _5, _6); - - using InputDomain = _1; - - template - VTKM_EXEC void operator()(const PointType& point, - const CellSetType& cellSet, - const CoordsPortalType& coords, - const LookupStructureType& lookupStruct, - vtkm::Id& cellId, - FloatVec3& parametricCoordinates) const - { - cellId = -1; - FloatVec3 p(static_cast(point[0]), - static_cast(point[1]), - static_cast(point[2])); - - const vtkm::exec::twolevelgrid::Grid& topLevelGrid = lookupStruct.TopLevel; - - DimVec3 binId3 = static_cast((p - topLevelGrid.Origin) / topLevelGrid.BinSize); - if (binId3[0] >= 0 && binId3[0] < topLevelGrid.Dimensions[0] && binId3[1] >= 0 && - binId3[1] < topLevelGrid.Dimensions[1] && binId3[2] >= 0 && - binId3[2] < topLevelGrid.Dimensions[2]) - { - vtkm::Id binId = ComputeFlatIndex(binId3, topLevelGrid.Dimensions); - - auto ldim = lookupStruct.LeafDimensions.Get(binId); - if (!ldim[0] || !ldim[1] || !ldim[2]) - { - return; - } - - auto leafGrid = ComputeLeafGrid(binId3, ldim, topLevelGrid); - - DimVec3 leafId3 = static_cast((p - leafGrid.Origin) / leafGrid.BinSize); - // precision issues may cause leafId3 to be out of range so clamp it - leafId3 = vtkm::Max(DimVec3(0), vtkm::Min(ldim - DimVec3(1), leafId3)); - - vtkm::Id leafStart = lookupStruct.LeafStartIndex.Get(binId); - vtkm::Id leafId = leafStart + ComputeFlatIndex(leafId3, leafGrid.Dimensions); - - vtkm::Id start = lookupStruct.CellStartIndex.Get(leafId); - vtkm::Id end = start + lookupStruct.CellCount.Get(leafId); - for (vtkm::Id i = start; i < end; ++i) - { - vtkm::Id cid = lookupStruct.CellIds.Get(i); - auto indices = cellSet.GetIndices(cid); - vtkm::VecFromPortalPermute pts(&indices, coords); - FloatVec3 pc; - if (PointInsideCell(p, cellSet.GetCellShape(cid), pts, *this, pc)) - { - cellId = cid; - parametricCoordinates = pc; - break; - } - } - } - } - }; - - /// Finds the containing cells for the given array of points. Returns the cell ids - /// in the `cellIds` arrays. If a cell could not be found due to the point being - /// outside all the cells or due to numerical errors, the cell id is set to -1. - /// Parametric coordinates of the point inside the cell is returned in the - /// `parametricCoords` array. - /// - template - void FindCells( - const vtkm::cont::ArrayHandle, PointStorageType>& points, - vtkm::cont::ArrayHandle& cellIds, - vtkm::cont::ArrayHandle& parametricCoords, - CellSetList cellSetTypes = CellSetList()) const - { - vtkm::worklet::DispatcherMapField dispatcher; - dispatcher.Invoke(points, - this->CellSet.ResetCellSetList(cellSetTypes), - this->Coordinates, - this->LookupStructure, - cellIds, - parametricCoords); - } - -private: - vtkm::FloatDefault DensityL1, DensityL2; - - vtkm::cont::DynamicCellSet CellSet; - vtkm::cont::CoordinateSystem Coordinates; - TwoLevelUniformGrid LookupStructure; -}; -} -} - -#endif // vtk_m_cont_CellLocatorTwoLevelUniformGrid_h diff --git a/vtkm/cont/CellLocatorUniformBins.cxx b/vtkm/cont/CellLocatorUniformBins.cxx new file mode 100644 index 000000000..00e1a447e --- /dev/null +++ b/vtkm/cont/CellLocatorUniformBins.cxx @@ -0,0 +1,531 @@ +//============================================================================ +// Copyright (c) Kitware, Inc. +// All rights reserved. +// See LICENSE.txt for details. +// 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. +// +// Copyright 2019 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +// Copyright 2019 UT-Battelle, LLC. +// Copyright 2019 Los Alamos National Security. +// +// Under the terms of Contract DE-NA0003525 with NTESS, +// the U.S. Government retains certain rights in this software. +// +// Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National +// Laboratory (LANL), the U.S. Government retains certain rights in +// this software. +//============================================================================ +#include + +#include +#include +#include +#include + +#include +#include +#include + +using namespace vtkm::internal::cl_uniform_bins; + +namespace +{ + +struct BinsBBox +{ + DimVec3 Min; + DimVec3 Max; + + VTKM_EXEC + bool Empty() const + { + return (this->Max[0] < this->Min[0]) || (this->Max[1] < this->Min[1]) || + (this->Max[2] < this->Min[2]); + } +}; + +VTKM_EXEC_CONT static DimVec3 ComputeGridDimension(vtkm::Id numberOfCells, + const FloatVec3& size, + vtkm::FloatDefault density) +{ + vtkm::FloatDefault nsides = 0.0f; + vtkm::FloatDefault volume = 1.0f; + vtkm::FloatDefault maxside = vtkm::Max(size[0], vtkm::Max(size[1], size[2])); + for (int i = 0; i < 3; ++i) + { + if (size[i] / maxside >= 1e-4f) + { + nsides += 1.0f; + volume *= size[i]; + } + } + + auto r = + vtkm::Pow((static_cast(numberOfCells) / (volume * density)), 1.0f / nsides); + return vtkm::Max(DimVec3(1), + DimVec3(static_cast(size[0] * r), + static_cast(size[1] * r), + static_cast(size[2] * r))); +} + +VTKM_EXEC static BinsBBox ComputeIntersectingBins(const Bounds cellBounds, const Grid& grid) +{ + auto minb = static_cast((cellBounds.Min - grid.Origin) / grid.BinSize); + auto maxb = static_cast((cellBounds.Max - grid.Origin) / grid.BinSize); + + return { vtkm::Max(DimVec3(0), minb), vtkm::Min(grid.Dimensions - DimVec3(1), maxb) }; +} + +VTKM_EXEC static vtkm::Id GetNumberOfBins(const BinsBBox& binsBBox) +{ + return binsBBox.Empty() ? 0 : ((binsBBox.Max[0] - binsBBox.Min[0] + 1) * + (binsBBox.Max[1] - binsBBox.Min[1] + 1) * + (binsBBox.Max[2] - binsBBox.Min[2] + 1)); +} + +class BBoxIterator +{ +public: + VTKM_EXEC_CONT BBoxIterator(const BinsBBox& bbox, const DimVec3& dim) + : BBox(bbox) + , Dim(dim) + , StepY(dim[0] - (bbox.Max[0] - bbox.Min[0] + 1)) + , StepZ((dim[0] * dim[1]) - ((bbox.Max[1] - bbox.Min[1] + 1) * dim[0])) + { + this->Init(); + } + + VTKM_EXEC_CONT void Init() + { + this->Idx = this->BBox.Min; + this->FlatIdx = ComputeFlatIndex(this->Idx, this->Dim); + this->DoneFlag = this->BBox.Empty(); + } + + VTKM_EXEC_CONT bool Done() const { return this->DoneFlag; } + + VTKM_EXEC_CONT void Next() + { + if (!this->DoneFlag) + { + ++this->Idx[0]; + this->FlatIdx += 1; + if (this->Idx[0] > this->BBox.Max[0]) + { + this->Idx[0] = this->BBox.Min[0]; + ++this->Idx[1]; + this->FlatIdx += this->StepY; + if (this->Idx[1] > this->BBox.Max[1]) + { + this->Idx[1] = this->BBox.Min[1]; + ++this->Idx[2]; + this->FlatIdx += this->StepZ; + if (this->Idx[2] > this->BBox.Max[2]) + { + this->DoneFlag = true; + } + } + } + } + } + + VTKM_EXEC_CONT const DimVec3& GetIdx() const { return this->Idx; } + + VTKM_EXEC_CONT vtkm::Id GetFlatIdx() const { return this->FlatIdx; } + +private: + const BinsBBox BBox; + const DimVec3 Dim; + const vtkm::Id StepY, StepZ; + + DimVec3 Idx; + vtkm::Id FlatIdx; + bool DoneFlag; +}; + +class CountBinsL1 : public vtkm::worklet::WorkletMapPointToCell +{ +public: + using ControlSignature = void(CellSetIn cellset, FieldInPoint coords, FieldOutCell bincount); + using ExecutionSignature = void(_2, _3); + + CountBinsL1(const Grid& grid) + : L1Grid(grid) + { + } + + template + VTKM_EXEC void operator()(const PointsVecType& points, vtkm::Id& numBins) const + { + auto cellBounds = ComputeCellBounds(points); + auto binsBBox = ComputeIntersectingBins(cellBounds, this->L1Grid); + numBins = GetNumberOfBins(binsBBox); + } + +private: + Grid L1Grid; +}; + +class FindBinsL1 : public vtkm::worklet::WorkletMapPointToCell +{ +public: + using ControlSignature = void(CellSetIn cellset, + FieldInPoint coords, + FieldInCell offsets, + WholeArrayOut binIds); + using ExecutionSignature = void(_2, _3, _4); + + FindBinsL1(const Grid& grid) + : L1Grid(grid) + { + } + + template + VTKM_EXEC void operator()(const PointsVecType& points, + vtkm::Id offset, + BinIdsPortalType& binIds) const + { + auto cellBounds = ComputeCellBounds(points); + auto binsBBox = ComputeIntersectingBins(cellBounds, this->L1Grid); + + for (BBoxIterator i(binsBBox, this->L1Grid.Dimensions); !i.Done(); i.Next()) + { + binIds.Set(offset, i.GetFlatIdx()); + ++offset; + } + } + +private: + Grid L1Grid; +}; + +class GenerateBinsL1 : public vtkm::worklet::WorkletMapField +{ +public: + using ControlSignature = void(FieldIn binIds, FieldIn cellCounts, WholeArrayOut dimensions); + using ExecutionSignature = void(_1, _2, _3); + + using InputDomain = _1; + + GenerateBinsL1(FloatVec3 size, vtkm::FloatDefault density) + : Size(size) + , Density(density) + { + } + + template + VTKM_EXEC void operator()(vtkm::Id binId, + vtkm::Id numCells, + OutputDimensionsPortal& dimensions) const + { + dimensions.Set(binId, ComputeGridDimension(numCells, this->Size, this->Density)); + } + +private: + FloatVec3 Size; + vtkm::FloatDefault Density; +}; + +class CountBinsL2 : public vtkm::worklet::WorkletMapPointToCell +{ +public: + using ControlSignature = void(CellSetIn cellset, + FieldInPoint coords, + WholeArrayIn binDimensions, + FieldOutCell bincount); + using ExecutionSignature = void(_2, _3, _4); + + CountBinsL2(const Grid& grid) + : L1Grid(grid) + { + } + + template + VTKM_EXEC void operator()(const PointsVecType& points, + const BinDimensionsPortalType& binDimensions, + vtkm::Id& numBins) const + { + auto cellBounds = ComputeCellBounds(points); + auto binsBBox = ComputeIntersectingBins(cellBounds, this->L1Grid); + + numBins = 0; + for (BBoxIterator i(binsBBox, this->L1Grid.Dimensions); !i.Done(); i.Next()) + { + Grid leaf = ComputeLeafGrid(i.GetIdx(), binDimensions.Get(i.GetFlatIdx()), this->L1Grid); + auto binsBBoxL2 = ComputeIntersectingBins(cellBounds, leaf); + numBins += GetNumberOfBins(binsBBoxL2); + } + } + +private: + Grid L1Grid; +}; + +class FindBinsL2 : public vtkm::worklet::WorkletMapPointToCell +{ +public: + using ControlSignature = void(CellSetIn cellset, + FieldInPoint coords, + WholeArrayIn binDimensions, + WholeArrayIn binStarts, + FieldInCell offsets, + WholeArrayOut binIds, + WholeArrayOut cellIds); + using ExecutionSignature = void(InputIndex, _2, _3, _4, _5, _6, _7); + + FindBinsL2(const Grid& grid) + : L1Grid(grid) + { + } + + template + VTKM_EXEC void operator()(vtkm::Id cellId, + const PointsVecType& points, + const BinDimensionsPortalType& binDimensions, + const BinStartsPortalType& binStarts, + vtkm::Id offset, + BinIdsPortalType binIds, + CellIdsPortalType cellIds) const + { + auto cellBounds = ComputeCellBounds(points); + auto binsBBox = ComputeIntersectingBins(cellBounds, this->L1Grid); + + for (BBoxIterator i(binsBBox, this->L1Grid.Dimensions); !i.Done(); i.Next()) + { + Grid leaf = ComputeLeafGrid(i.GetIdx(), binDimensions.Get(i.GetFlatIdx()), this->L1Grid); + auto binsBBoxL2 = ComputeIntersectingBins(cellBounds, leaf); + vtkm::Id leafStart = binStarts.Get(i.GetFlatIdx()); + + for (BBoxIterator j(binsBBoxL2, leaf.Dimensions); !j.Done(); j.Next()) + { + binIds.Set(offset, leafStart + j.GetFlatIdx()); + cellIds.Set(offset, cellId); + ++offset; + } + } + } + +private: + Grid L1Grid; +}; + +class GenerateBinsL2 : public vtkm::worklet::WorkletMapField +{ +public: + using ControlSignature = void(FieldIn binIds, + FieldIn startsIn, + FieldIn countsIn, + WholeArrayOut startsOut, + WholeArrayOut countsOut); + using ExecutionSignature = void(_1, _2, _3, _4, _5); + + using InputDomain = _1; + + template + VTKM_EXEC void operator()(vtkm::Id binIndex, + vtkm::Id start, + vtkm::Id count, + CellStartsPortalType& cellStarts, + CellCountsPortalType& cellCounts) const + { + cellStarts.Set(binIndex, start); + cellCounts.Set(binIndex, count); + } +}; + +struct DimensionsToCount +{ + VTKM_EXEC vtkm::Id operator()(const DimVec3& dim) const { return dim[0] * dim[1] * dim[2]; } +}; + +} // anonymous namespace + +namespace vtkm +{ +namespace cont +{ + +//---------------------------------------------------------------------------- +/// Builds the cell locator lookup structure +/// +VTKM_CONT void CellLocatorUniformBins::Build() +{ + vtkm::worklet::Invoker invoke; + + auto cellset = this->GetCellSet(); + const auto& coords = this->GetCoordinates(); + + // 1: Compute the top level grid + auto bounds = coords.GetBounds(); + FloatVec3 bmin(static_cast(bounds.X.Min), + static_cast(bounds.Y.Min), + static_cast(bounds.Z.Min)); + FloatVec3 bmax(static_cast(bounds.X.Max), + static_cast(bounds.Y.Max), + static_cast(bounds.Z.Max)); + auto size = bmax - bmin; + auto fudge = vtkm::Max(FloatVec3(1e-6f), size * 1e-4f); + size += 2.0f * fudge; + + this->TopLevel.Dimensions = + ComputeGridDimension(cellset.GetNumberOfCells(), size, this->DensityL1); + this->TopLevel.Origin = bmin - fudge; + this->TopLevel.BinSize = size / static_cast(this->TopLevel.Dimensions); + + // 2: For each cell, find the number of top level bins they intersect + vtkm::cont::ArrayHandle binCounts; + CountBinsL1 countL1(this->TopLevel); + invoke(countL1, cellset, coords, binCounts); + + // 3: Total number of unique (cell, bin) pairs (for pre-allocating arrays) + vtkm::Id numPairsL1 = vtkm::cont::Algorithm::ScanExclusive(binCounts, binCounts); + + // 4: For each cell find the top level bins that intersect it + vtkm::cont::ArrayHandle binIds; + binIds.Allocate(numPairsL1); + FindBinsL1 findL1(this->TopLevel); + invoke(findL1, cellset, coords, binCounts, binIds); + binCounts.ReleaseResources(); + + // 5: From above, find the number of cells that intersect each top level bin + vtkm::cont::Algorithm::Sort(binIds); + vtkm::cont::ArrayHandle bins; + vtkm::cont::ArrayHandle cellsPerBin; + vtkm::cont::Algorithm::ReduceByKey(binIds, + vtkm::cont::make_ArrayHandleConstant(vtkm::Id(1), numPairsL1), + bins, + cellsPerBin, + vtkm::Sum()); + binIds.ReleaseResources(); + + // 6: Compute level-2 dimensions + vtkm::Id numberOfBins = + this->TopLevel.Dimensions[0] * this->TopLevel.Dimensions[1] * this->TopLevel.Dimensions[2]; + vtkm::cont::ArrayCopy(vtkm::cont::make_ArrayHandleConstant(DimVec3(0), numberOfBins), + this->LeafDimensions); + GenerateBinsL1 generateL1(this->TopLevel.BinSize, this->DensityL2); + invoke(generateL1, bins, cellsPerBin, this->LeafDimensions); + bins.ReleaseResources(); + cellsPerBin.ReleaseResources(); + + // 7: Compute number of level-2 bins + vtkm::Id numberOfLeaves = vtkm::cont::Algorithm::ScanExclusive( + vtkm::cont::make_ArrayHandleTransform(this->LeafDimensions, DimensionsToCount()), + this->LeafStartIndex); + + + // 8: For each cell, find the number of l2 bins they intersect + CountBinsL2 countL2(this->TopLevel); + invoke(countL2, cellset, coords, this->LeafDimensions, binCounts); + + // 9: Total number of unique (cell, bin) pairs (for pre-allocating arrays) + vtkm::Id numPairsL2 = vtkm::cont::Algorithm::ScanExclusive(binCounts, binCounts); + + // 10: For each cell, find the l2 bins they intersect + binIds.Allocate(numPairsL2); + this->CellIds.Allocate(numPairsL2); + FindBinsL2 findL2(this->TopLevel); + invoke(findL2, + cellset, + coords, + this->LeafDimensions, + this->LeafStartIndex, + binCounts, + binIds, + this->CellIds); + binCounts.ReleaseResources(); + + // 11: From above, find the cells that each l2 bin intersects + vtkm::cont::Algorithm::SortByKey(binIds, this->CellIds); + vtkm::cont::Algorithm::ReduceByKey(binIds, + vtkm::cont::make_ArrayHandleConstant(vtkm::Id(1), numPairsL2), + bins, + cellsPerBin, + vtkm::Sum()); + binIds.ReleaseResources(); + + // 12: Generate the leaf bin arrays + vtkm::cont::ArrayHandle cellsStart; + vtkm::cont::Algorithm::ScanExclusive(cellsPerBin, cellsStart); + + vtkm::cont::ArrayCopy(vtkm::cont::ArrayHandleConstant(0, numberOfLeaves), + this->CellStartIndex); + vtkm::cont::ArrayCopy(vtkm::cont::ArrayHandleConstant(0, numberOfLeaves), + this->CellCount); + invoke(GenerateBinsL2{}, bins, cellsStart, cellsPerBin, this->CellStartIndex, this->CellCount); +} + +//---------------------------------------------------------------------------- +struct CellLocatorUniformBins::MakeExecObject +{ + template + VTKM_CONT void operator()(const CellSetType& cellSet, + DeviceAdapter, + const CellLocatorUniformBins& self) const + { + auto execObject = + new vtkm::exec::CellLocatorUniformBins(self.TopLevel, + self.LeafDimensions, + self.LeafStartIndex, + self.CellStartIndex, + self.CellCount, + self.CellIds, + cellSet, + self.GetCoordinates()); + self.ExecutionObjectHandle.Reset(execObject); + } +}; + +struct CellLocatorUniformBins::PrepareForExecutionFunctor +{ + template + VTKM_CONT bool operator()(DeviceAdapter, const CellLocatorUniformBins& self) const + { + self.GetCellSet().CastAndCall(MakeExecObject{}, DeviceAdapter{}, self); + return true; + } +}; + +VTKM_CONT const vtkm::exec::CellLocator* CellLocatorUniformBins::PrepareForExecution( + vtkm::cont::DeviceAdapterId device) const +{ + if (!vtkm::cont::TryExecuteOnDevice(device, PrepareForExecutionFunctor(), *this)) + { + throwFailedRuntimeDeviceTransfer("CellLocatorUniformBins", device); + } + return this->ExecutionObjectHandle.PrepareForExecution(device); +} + +//---------------------------------------------------------------------------- +void CellLocatorUniformBins::PrintSummary(std::ostream& out) const +{ + out << "DensityL1: " << this->DensityL1 << "\n"; + out << "DensityL2: " << this->DensityL2 << "\n"; + out << "Input CellSet: \n"; + this->GetCellSet().PrintSummary(out); + out << "Input Coordinates: \n"; + this->GetCoordinates().PrintSummary(out); + out << "LookupStructure:\n"; + out << " TopLevelGrid\n"; + out << " Dimensions: " << this->TopLevel.Dimensions << "\n"; + out << " Origin: " << this->TopLevel.Origin << "\n"; + out << " BinSize: " << this->TopLevel.BinSize << "\n"; + out << " LeafDimensions:\n"; + vtkm::cont::printSummary_ArrayHandle(this->LeafDimensions, out); + out << " LeafStartIndex:\n"; + vtkm::cont::printSummary_ArrayHandle(this->LeafStartIndex, out); + out << " CellStartIndex:\n"; + vtkm::cont::printSummary_ArrayHandle(this->CellStartIndex, out); + out << " CellCount:\n"; + vtkm::cont::printSummary_ArrayHandle(this->CellCount, out); + out << " CellIds:\n"; + vtkm::cont::printSummary_ArrayHandle(this->CellIds, out); +} +} +} // vtkm::cont diff --git a/vtkm/cont/CellLocatorUniformBins.h b/vtkm/cont/CellLocatorUniformBins.h new file mode 100644 index 000000000..1a3c3589a --- /dev/null +++ b/vtkm/cont/CellLocatorUniformBins.h @@ -0,0 +1,284 @@ +//============================================================================ +// Copyright (c) Kitware, Inc. +// All rights reserved. +// See LICENSE.txt for details. +// 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. +// +// Copyright 2017 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +// Copyright 2017 UT-Battelle, LLC. +// Copyright 2017 Los Alamos National Security. +// +// Under the terms of Contract DE-NA0003525 with NTESS, +// the U.S. Government retains certain rights in this software. +// +// Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National +// Laboratory (LANL), the U.S. Government retains certain rights in +// this software. +//============================================================================ +#ifndef vtk_m_cont_CellLocatorUniformBins_h +#define vtk_m_cont_CellLocatorUniformBins_h + +#include +#include +#include + +#include +#include + +#include +#include +#include +#include + +namespace vtkm +{ +namespace internal +{ +namespace cl_uniform_bins +{ + +using DimensionType = vtkm::Int16; +using DimVec3 = vtkm::Vec; +using FloatVec3 = vtkm::Vec; + +struct Grid +{ + DimVec3 Dimensions; + FloatVec3 Origin; + FloatVec3 BinSize; +}; + +struct Bounds +{ + FloatVec3 Min; + FloatVec3 Max; +}; + +VTKM_EXEC inline vtkm::Id ComputeFlatIndex(const DimVec3& idx, const DimVec3 dim) +{ + return idx[0] + (dim[0] * (idx[1] + (dim[1] * idx[2]))); +} + +VTKM_EXEC inline Grid ComputeLeafGrid(const DimVec3& idx, const DimVec3& dim, const Grid& l1Grid) +{ + return { dim, + l1Grid.Origin + (static_cast(idx) * l1Grid.BinSize), + l1Grid.BinSize / static_cast(dim) }; +} + +template +VTKM_EXEC inline Bounds ComputeCellBounds(const PointsVecType& points) +{ + using CoordsType = typename vtkm::VecTraits::ComponentType; + auto numPoints = vtkm::VecTraits::GetNumberOfComponents(points); + + CoordsType minp = points[0], maxp = points[0]; + for (vtkm::IdComponent i = 1; i < numPoints; ++i) + { + minp = vtkm::Min(minp, points[i]); + maxp = vtkm::Max(maxp, points[i]); + } + + return { FloatVec3(minp), FloatVec3(maxp) }; +} +} +} +} // vtkm::internal::cl_uniform_bins + +namespace vtkm +{ +namespace exec +{ + +//-------------------------------------------------------------------- +template +class CellLocatorUniformBins : public vtkm::exec::CellLocator +{ +private: + using DimVec3 = vtkm::internal::cl_uniform_bins::DimVec3; + using FloatVec3 = vtkm::internal::cl_uniform_bins::FloatVec3; + + template + using ArrayPortalConst = + typename vtkm::cont::ArrayHandle::template ExecutionTypes::PortalConst; + + using CoordsPortalType = + decltype(vtkm::cont::ArrayHandleVirtualCoordinates{}.PrepareForInput(DeviceAdapter{})); + + using CellSetP2CExecType = + decltype(std::declval().PrepareForInput(DeviceAdapter{}, + vtkm::TopologyElementTagPoint{}, + vtkm::TopologyElementTagCell{})); + + // TODO: This function may return false positives for non 3D cells as the + // tests are done on the projection of the point on the cell. Extra checks + // should be added to test if the point actually falls on the cell. + template + VTKM_EXEC static bool PointInsideCell(FloatVec3 point, + CellShapeTag cellShape, + CoordsType cellPoints, + const vtkm::exec::FunctorBase& worklet, + FloatVec3& parametricCoordinates) + { + auto bounds = vtkm::internal::cl_uniform_bins::ComputeCellBounds(cellPoints); + if (point[0] >= bounds.Min[0] && point[0] <= bounds.Max[0] && point[1] >= bounds.Min[1] && + point[1] <= bounds.Max[1] && point[2] >= bounds.Min[2] && point[2] <= bounds.Max[2]) + { + bool success = false; + parametricCoordinates = vtkm::exec::WorldCoordinatesToParametricCoordinates( + cellPoints, point, cellShape, success, worklet); + return success && vtkm::exec::CellInside(parametricCoordinates, cellShape); + } + return false; + } + +public: + VTKM_CONT CellLocatorUniformBins(const vtkm::internal::cl_uniform_bins::Grid& topLevelGrid, + const vtkm::cont::ArrayHandle& leafDimensions, + const vtkm::cont::ArrayHandle& leafStartIndex, + const vtkm::cont::ArrayHandle& cellStartIndex, + const vtkm::cont::ArrayHandle& cellCount, + const vtkm::cont::ArrayHandle& cellIds, + const CellSetType& cellSet, + const vtkm::cont::CoordinateSystem& coords) + : TopLevel(topLevelGrid) + , LeafDimensions(leafDimensions.PrepareForInput(DeviceAdapter{})) + , LeafStartIndex(leafStartIndex.PrepareForInput(DeviceAdapter{})) + , CellStartIndex(cellStartIndex.PrepareForInput(DeviceAdapter{})) + , CellCount(cellCount.PrepareForInput(DeviceAdapter{})) + , CellIds(cellIds.PrepareForInput(DeviceAdapter{})) + , CellSet(cellSet.PrepareForInput(DeviceAdapter{}, + vtkm::TopologyElementTagPoint{}, + vtkm::TopologyElementTagCell{})) + , Coords(coords.GetData().PrepareForInput(DeviceAdapter{})) + { + } + + VTKM_EXEC + void FindCell(const FloatVec3& point, + vtkm::Id& cellId, + FloatVec3& parametric, + const vtkm::exec::FunctorBase& worklet) const override + { + using namespace vtkm::internal::cl_uniform_bins; + + cellId = -1; + + DimVec3 binId3 = static_cast((point - this->TopLevel.Origin) / this->TopLevel.BinSize); + if (binId3[0] >= 0 && binId3[0] < this->TopLevel.Dimensions[0] && binId3[1] >= 0 && + binId3[1] < this->TopLevel.Dimensions[1] && binId3[2] >= 0 && + binId3[2] < this->TopLevel.Dimensions[2]) + { + vtkm::Id binId = ComputeFlatIndex(binId3, this->TopLevel.Dimensions); + + auto ldim = this->LeafDimensions.Get(binId); + if (!ldim[0] || !ldim[1] || !ldim[2]) + { + return; + } + + auto leafGrid = ComputeLeafGrid(binId3, ldim, this->TopLevel); + + DimVec3 leafId3 = static_cast((point - leafGrid.Origin) / leafGrid.BinSize); + // precision issues may cause leafId3 to be out of range so clamp it + leafId3 = vtkm::Max(DimVec3(0), vtkm::Min(ldim - DimVec3(1), leafId3)); + + vtkm::Id leafStart = this->LeafStartIndex.Get(binId); + vtkm::Id leafId = leafStart + ComputeFlatIndex(leafId3, leafGrid.Dimensions); + + vtkm::Id start = this->CellStartIndex.Get(leafId); + vtkm::Id end = start + this->CellCount.Get(leafId); + for (vtkm::Id i = start; i < end; ++i) + { + vtkm::Id cid = this->CellIds.Get(i); + auto indices = this->CellSet.GetIndices(cid); + auto pts = vtkm::make_VecFromPortalPermute(&indices, this->Coords); + FloatVec3 pc; + if (PointInsideCell(point, this->CellSet.GetCellShape(cid), pts, worklet, pc)) + { + cellId = cid; + parametric = pc; + break; + } + } + } + } + +private: + vtkm::internal::cl_uniform_bins::Grid TopLevel; + + ArrayPortalConst LeafDimensions; + ArrayPortalConst LeafStartIndex; + + ArrayPortalConst CellStartIndex; + ArrayPortalConst CellCount; + ArrayPortalConst CellIds; + + CellSetP2CExecType CellSet; + CoordsPortalType Coords; +}; +} +} // vtkm::exec + + +namespace vtkm +{ +namespace cont +{ + +//---------------------------------------------------------------------------- +class VTKM_CONT_EXPORT CellLocatorUniformBins : public vtkm::cont::CellLocator +{ +public: + CellLocatorUniformBins() + : DensityL1(32.0f) + , DensityL2(2.0f) + { + } + + /// Get/Set the desired approximate number of cells per level 1 bin + /// + void SetDensityL1(vtkm::FloatDefault val) + { + this->DensityL1 = val; + this->SetModified(); + } + vtkm::FloatDefault GetDensityL1() const { return this->DensityL1; } + + /// Get/Set the desired approximate number of cells per level 1 bin + /// + void SetDensityL2(vtkm::FloatDefault val) + { + this->DensityL2 = val; + this->SetModified(); + } + vtkm::FloatDefault GetDensityL2() const { return this->DensityL2; } + + void PrintSummary(std::ostream& out) const; + + const vtkm::exec::CellLocator* PrepareForExecution( + vtkm::cont::DeviceAdapterId device) const override; + +private: + VTKM_CONT void Build() override; + + vtkm::FloatDefault DensityL1, DensityL2; + + vtkm::internal::cl_uniform_bins::Grid TopLevel; + vtkm::cont::ArrayHandle LeafDimensions; + vtkm::cont::ArrayHandle LeafStartIndex; + vtkm::cont::ArrayHandle CellStartIndex; + vtkm::cont::ArrayHandle CellCount; + vtkm::cont::ArrayHandle CellIds; + + mutable vtkm::cont::VirtualObjectHandle ExecutionObjectHandle; + + struct MakeExecObject; + struct PrepareForExecutionFunctor; +}; +} +} // vtkm::cont + +#endif // vtk_m_cont_CellLocatorUniformBins_h diff --git a/vtkm/cont/cuda/testing/CMakeLists.txt b/vtkm/cont/cuda/testing/CMakeLists.txt index cd274f387..49cff5989 100644 --- a/vtkm/cont/cuda/testing/CMakeLists.txt +++ b/vtkm/cont/cuda/testing/CMakeLists.txt @@ -23,7 +23,7 @@ set(unit_tests UnitTestCudaArrayHandleFancy.cu UnitTestCudaArrayHandleVirtualCoordinates.cu UnitTestCudaCellLocatorRectilinearGrid.cu - UnitTestCudaCellLocatorTwoLevelUniformGrid.cu + UnitTestCudaCellLocatorUniformBins.cu UnitTestCudaCellLocatorUniformGrid.cu UnitTestCudaComputeRange.cu UnitTestCudaColorTable.cu diff --git a/vtkm/cont/cuda/testing/UnitTestCudaCellLocatorTwoLevelUniformGrid.cu b/vtkm/cont/cuda/testing/UnitTestCudaCellLocatorUniformBins.cu similarity index 85% rename from vtkm/cont/cuda/testing/UnitTestCudaCellLocatorTwoLevelUniformGrid.cu rename to vtkm/cont/cuda/testing/UnitTestCudaCellLocatorUniformBins.cu index fc52b2027..41a1365fc 100644 --- a/vtkm/cont/cuda/testing/UnitTestCudaCellLocatorTwoLevelUniformGrid.cu +++ b/vtkm/cont/cuda/testing/UnitTestCudaCellLocatorUniformBins.cu @@ -23,12 +23,12 @@ // for a part of an operation where the TBB device was specified. #define VTKM_DEVICE_ADAPTER VTKM_DEVICE_ADAPTER_ERROR -#include +#include -int UnitTestCudaCellLocatorTwoLevelUniformGrid(int argc, char* argv[]) +int UnitTestCudaCellLocatorUniformBins(int argc, char* argv[]) { auto tracker = vtkm::cont::GetRuntimeDeviceTracker(); tracker.ForceDevice(vtkm::cont::DeviceAdapterTagCuda{}); return vtkm::cont::testing::Testing::Run( - TestingCellLocatorTwoLevelUniformGrid, argc, argv); + TestingCellLocatorUniformBins, argc, argv); } diff --git a/vtkm/cont/openmp/testing/CMakeLists.txt b/vtkm/cont/openmp/testing/CMakeLists.txt index 9659ea7ac..0a10b09bb 100644 --- a/vtkm/cont/openmp/testing/CMakeLists.txt +++ b/vtkm/cont/openmp/testing/CMakeLists.txt @@ -23,7 +23,7 @@ set(unit_tests UnitTestOpenMPArrayHandleFancy.cxx UnitTestOpenMPArrayHandleVirtualCoordinates.cxx UnitTestOpenMPCellLocatorRectilinearGrid.cxx - UnitTestOpenMPCellLocatorTwoLevelUniformGrid.cxx + UnitTestOpenMPCellLocatorUniformBins.cxx UnitTestOpenMPCellLocatorUniformGrid.cxx UnitTestOpenMPColorTable.cxx UnitTestOpenMPComputeRange.cxx diff --git a/vtkm/cont/openmp/testing/UnitTestOpenMPCellLocatorTwoLevelUniformGrid.cxx b/vtkm/cont/openmp/testing/UnitTestOpenMPCellLocatorUniformBins.cxx similarity index 83% rename from vtkm/cont/openmp/testing/UnitTestOpenMPCellLocatorTwoLevelUniformGrid.cxx rename to vtkm/cont/openmp/testing/UnitTestOpenMPCellLocatorUniformBins.cxx index 2817f3431..6f8c6b4a4 100644 --- a/vtkm/cont/openmp/testing/UnitTestOpenMPCellLocatorTwoLevelUniformGrid.cxx +++ b/vtkm/cont/openmp/testing/UnitTestOpenMPCellLocatorUniformBins.cxx @@ -21,12 +21,12 @@ #define VTKM_DEVICE_ADAPTER VTKM_DEVICE_ADAPTER_ERROR #include -#include +#include -int UnitTestOpenMPCellLocatorTwoLevelUniformGrid(int argc, char* argv[]) +int UnitTestOpenMPCellLocatorUniformBins(int argc, char* argv[]) { auto tracker = vtkm::cont::GetRuntimeDeviceTracker(); tracker.ForceDevice(vtkm::cont::DeviceAdapterTagOpenMP{}); return vtkm::cont::testing::Testing::Run( - TestingCellLocatorTwoLevelUniformGrid, argc, argv); + TestingCellLocatorUniformBins, argc, argv); } diff --git a/vtkm/cont/serial/testing/CMakeLists.txt b/vtkm/cont/serial/testing/CMakeLists.txt index fcea48e6b..9ed2ebf7b 100644 --- a/vtkm/cont/serial/testing/CMakeLists.txt +++ b/vtkm/cont/serial/testing/CMakeLists.txt @@ -23,7 +23,7 @@ set(unit_tests UnitTestSerialArrayHandleFancy.cxx UnitTestSerialArrayHandleVirtualCoordinates.cxx UnitTestSerialCellLocatorRectilinearGrid.cxx - UnitTestSerialCellLocatorTwoLevelUniformGrid.cxx + UnitTestSerialCellLocatorUniformBins.cxx UnitTestSerialCellLocatorUniformGrid.cxx UnitTestSerialComputeRange.cxx UnitTestSerialColorTable.cxx diff --git a/vtkm/cont/serial/testing/UnitTestSerialCellLocatorTwoLevelUniformGrid.cxx b/vtkm/cont/serial/testing/UnitTestSerialCellLocatorUniformBins.cxx similarity index 85% rename from vtkm/cont/serial/testing/UnitTestSerialCellLocatorTwoLevelUniformGrid.cxx rename to vtkm/cont/serial/testing/UnitTestSerialCellLocatorUniformBins.cxx index f3568f00e..a499b930b 100644 --- a/vtkm/cont/serial/testing/UnitTestSerialCellLocatorTwoLevelUniformGrid.cxx +++ b/vtkm/cont/serial/testing/UnitTestSerialCellLocatorUniformBins.cxx @@ -24,12 +24,12 @@ #define VTKM_DEVICE_ADAPTER VTKM_DEVICE_ADAPTER_ERROR #include -#include +#include -int UnitTestSerialCellLocatorTwoLevelUniformGrid(int argc, char* argv[]) +int UnitTestSerialCellLocatorUniformBins(int argc, char* argv[]) { auto tracker = vtkm::cont::GetRuntimeDeviceTracker(); tracker.ForceDevice(vtkm::cont::DeviceAdapterTagSerial{}); return vtkm::cont::testing::Testing::Run( - TestingCellLocatorTwoLevelUniformGrid, argc, argv); + TestingCellLocatorUniformBins, argc, argv); } diff --git a/vtkm/cont/tbb/testing/CMakeLists.txt b/vtkm/cont/tbb/testing/CMakeLists.txt index 9a498f6fb..9aea9163a 100644 --- a/vtkm/cont/tbb/testing/CMakeLists.txt +++ b/vtkm/cont/tbb/testing/CMakeLists.txt @@ -23,7 +23,7 @@ set(unit_tests UnitTestTBBArrayHandleFancy.cxx UnitTestTBBArrayHandleVirtualCoordinates.cxx UnitTestTBBCellLocatorRectilinearGrid.cxx - UnitTestTBBCellLocatorTwoLevelUniformGrid.cxx + UnitTestTBBCellLocatorUniformBins.cxx UnitTestTBBCellLocatorUniformGrid.cxx UnitTestTBBColorTable.cxx UnitTestTBBComputeRange.cxx diff --git a/vtkm/cont/tbb/testing/UnitTestTBBCellLocatorTwoLevelUniformGrid.cxx b/vtkm/cont/tbb/testing/UnitTestTBBCellLocatorUniformBins.cxx similarity index 85% rename from vtkm/cont/tbb/testing/UnitTestTBBCellLocatorTwoLevelUniformGrid.cxx rename to vtkm/cont/tbb/testing/UnitTestTBBCellLocatorUniformBins.cxx index cca641eb3..03b76756c 100644 --- a/vtkm/cont/tbb/testing/UnitTestTBBCellLocatorTwoLevelUniformGrid.cxx +++ b/vtkm/cont/tbb/testing/UnitTestTBBCellLocatorUniformBins.cxx @@ -23,12 +23,12 @@ // for a part of an operation where the TBB device was specified. #define VTKM_DEVICE_ADAPTER VTKM_DEVICE_ADAPTER_ERROR -#include +#include -int UnitTestTBBCellLocatorTwoLevelUniformGrid(int argc, char* argv[]) +int UnitTestTBBCellLocatorUniformBins(int argc, char* argv[]) { auto tracker = vtkm::cont::GetRuntimeDeviceTracker(); tracker.ForceDevice(vtkm::cont::DeviceAdapterTagTBB{}); return vtkm::cont::testing::Testing::Run( - TestingCellLocatorTwoLevelUniformGrid, argc, argv); + TestingCellLocatorUniformBins, argc, argv); } diff --git a/vtkm/cont/testing/CMakeLists.txt b/vtkm/cont/testing/CMakeLists.txt index ccaebf2ee..9ce50eecb 100644 --- a/vtkm/cont/testing/CMakeLists.txt +++ b/vtkm/cont/testing/CMakeLists.txt @@ -25,7 +25,7 @@ set(headers TestingArrayHandles.h TestingArrayHandleVirtualCoordinates.h TestingCellLocatorRectilinearGrid.h - TestingCellLocatorTwoLevelUniformGrid.h + TestingCellLocatorUniformBins.h TestingCellLocatorUniformGrid.h TestingColorTable.h TestingComputeRange.h diff --git a/vtkm/cont/testing/TestingCellLocatorTwoLevelUniformGrid.h b/vtkm/cont/testing/TestingCellLocatorUniformBins.h similarity index 86% rename from vtkm/cont/testing/TestingCellLocatorTwoLevelUniformGrid.h rename to vtkm/cont/testing/TestingCellLocatorUniformBins.h index 9e9096d12..6eaa36851 100644 --- a/vtkm/cont/testing/TestingCellLocatorTwoLevelUniformGrid.h +++ b/vtkm/cont/testing/TestingCellLocatorUniformBins.h @@ -17,20 +17,22 @@ // Laboratory (LANL), the U.S. Government retains certain rights in // this software. //============================================================================ -#ifndef vtk_m_cont_testing_TestingCellLocatorTwoLevelUniformGrid_h -#define vtk_m_cont_testing_TestingCellLocatorTwoLevelUniformGrid_h +#ifndef vtk_m_cont_testing_TestingCellLocatorUniformBins_h +#define vtk_m_cont_testing_TestingCellLocatorUniformBins_h #include -#include +#include #include #include #include +#include #include #include #include #include +#include #include #include @@ -185,6 +187,25 @@ void GenerateRandomInput(const vtkm::cont::DataSet& ds, dispatcher.Invoke(ds.GetCellSet(), ds.GetCoordinateSystem().GetData(), pcoords, wcoords); } +class FindCellWorklet : public vtkm::worklet::WorkletMapField +{ +public: + using ControlSignature = void(FieldIn points, + ExecObject locator, + FieldOut cellIds, + FieldOut pcoords); + using ExecutionSignature = void(_1, _2, _3, _4); + + template + VTKM_EXEC void operator()(const vtkm::Vec& point, + const LocatorType& locator, + vtkm::Id& cellId, + vtkm::Vec& pcoords) const + { + locator->FindCell(point, cellId, pcoords, *this); + } +}; + template void TestCellLocator(const vtkm::Vec& dim, vtkm::Id numberOfPoints) { @@ -193,12 +214,12 @@ void TestCellLocator(const vtkm::Vec& dim, vtkm::Id number std::cout << "Testing " << DIMENSIONS << "D dataset with " << ds.GetCellSet().GetNumberOfCells() << " cells\n"; - vtkm::cont::CellLocatorTwoLevelUniformGrid locator; + vtkm::cont::CellLocatorUniformBins locator; locator.SetDensityL1(64.0f); locator.SetDensityL2(1.0f); locator.SetCellSet(ds.GetCellSet()); locator.SetCoordinates(ds.GetCoordinateSystem()); - locator.Build(); + locator.Update(); vtkm::cont::ArrayHandle expCellIds; vtkm::cont::ArrayHandle expPCoords; @@ -208,7 +229,9 @@ void TestCellLocator(const vtkm::Vec& dim, vtkm::Id number std::cout << "Finding cells for " << numberOfPoints << " points\n"; vtkm::cont::ArrayHandle cellIds; vtkm::cont::ArrayHandle pcoords; - locator.FindCells(points, cellIds, pcoords); + + vtkm::worklet::DispatcherMapField dispatcher; + dispatcher.Invoke(points, locator, cellIds, pcoords); for (vtkm::Id i = 0; i < numberOfPoints; ++i) { @@ -225,7 +248,7 @@ void TestCellLocator(const vtkm::Vec& dim, vtkm::Id number } // anonymous template -void TestingCellLocatorTwoLevelUniformGrid() +void TestingCellLocatorUniformBins() { vtkm::cont::GetRuntimeDeviceTracker().ForceDevice(DeviceAdapter()); @@ -237,4 +260,4 @@ void TestingCellLocatorTwoLevelUniformGrid() TestCellLocator(vtkm::Id2(18), 512); // 2D dataset } -#endif // vtk_m_cont_testing_TestingCellLocatorTwoLevelUniformGrid_h +#endif // vtk_m_cont_testing_TestingCellLocatorUniformBins_h