vtk-m/vtkm/exec/CellLocatorTwoLevel.h
Kenneth Moreland 3e1339f9a7 Remove deprecated features from VTK-m
With the major revision 2.0 of VTK-m, many items previously marked as
deprecated were removed. If updating to a new version of VTK-m, it is
recommended to first update to VTK-m 1.9, which will include the deprecated
features but provide warnings (with the right compiler) that will point to
the replacement code. Once the deprecations have been fixed, updating to
2.0 should be smoother.
2022-11-17 07:12:31 -06:00

305 lines
10 KiB
C++

//============================================================================
// 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.
//============================================================================
#ifndef vtk_m_exec_CellLocatorTwoLevel_h
#define vtk_m_exec_CellLocatorTwoLevel_h
#include <vtkm/exec/CellInside.h>
#include <vtkm/exec/ParametricCoordinates.h>
#include <vtkm/cont/ArrayHandle.h>
#include <vtkm/cont/CoordinateSystem.h>
#include <vtkm/Math.h>
#include <vtkm/TopologyElementTag.h>
#include <vtkm/Types.h>
#include <vtkm/VecFromPortalPermute.h>
#include <vtkm/VecTraits.h>
namespace vtkm
{
namespace internal
{
namespace cl_uniform_bins
{
using DimensionType = vtkm::Int16;
using DimVec3 = vtkm::Vec<DimensionType, 3>;
using FloatVec3 = vtkm::Vec3f;
struct Grid
{
DimVec3 Dimensions;
// Bug in CUDA 9.2 where having this gap for alignment was for some reason setting garbage
// in a union with other cell locators (or perhaps not properly copying data). This appears
// to be fixed by CUDA 10.2.
DimensionType Padding;
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,
0,
l1Grid.Origin + (static_cast<FloatVec3>(idx) * l1Grid.BinSize),
l1Grid.BinSize / static_cast<FloatVec3>(dim) };
}
template <typename PointsVecType>
VTKM_EXEC inline Bounds ComputeCellBounds(const PointsVecType& points)
{
using CoordsType = typename vtkm::VecTraits<PointsVecType>::ComponentType;
auto numPoints = vtkm::VecTraits<PointsVecType>::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 <typename CellStructureType>
class VTKM_ALWAYS_EXPORT CellLocatorTwoLevel
{
private:
using DimVec3 = vtkm::internal::cl_uniform_bins::DimVec3;
using FloatVec3 = vtkm::internal::cl_uniform_bins::FloatVec3;
template <typename T>
using ReadPortal = typename vtkm::cont::ArrayHandle<T>::ReadPortalType;
using CoordsPortalType =
typename vtkm::cont::CoordinateSystem::MultiplexerArrayType::ReadPortalType;
// 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 <typename CellShapeTag, typename CoordsType>
VTKM_EXEC static vtkm::ErrorCode PointInsideCell(FloatVec3 point,
CellShapeTag cellShape,
CoordsType cellPoints,
FloatVec3& parametricCoordinates,
bool& inside)
{
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])
{
VTKM_RETURN_ON_ERROR(vtkm::exec::WorldCoordinatesToParametricCoordinates(
cellPoints, point, cellShape, parametricCoordinates));
inside = vtkm::exec::CellInside(parametricCoordinates, cellShape);
}
else
{
inside = false;
}
// Return success error code even point is not inside this cell
return vtkm::ErrorCode::Success;
}
public:
template <typename CellSetType>
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)
: TopLevel(topLevelGrid)
, LeafDimensions(leafDimensions.PrepareForInput(device, token))
, LeafStartIndex(leafStartIndex.PrepareForInput(device, token))
, CellStartIndex(cellStartIndex.PrepareForInput(device, token))
, CellCount(cellCount.PrepareForInput(device, token))
, CellIds(cellIds.PrepareForInput(device, token))
, CellSet(cellSet.PrepareForInput(device,
vtkm::TopologyElementTagCell{},
vtkm::TopologyElementTagPoint{},
token))
, Coords(coords.GetDataAsMultiplexer().PrepareForInput(device, token))
{
}
struct LastCell
{
vtkm::Id CellId = -1;
vtkm::Id LeafIdx = -1;
};
VTKM_EXEC
vtkm::ErrorCode FindCell(const FloatVec3& point, vtkm::Id& cellId, FloatVec3& parametric) const
{
LastCell lastCell;
return this->FindCellImpl(point, cellId, parametric, lastCell);
}
VTKM_EXEC
vtkm::ErrorCode FindCell(const FloatVec3& point,
vtkm::Id& cellId,
FloatVec3& parametric,
LastCell& lastCell) const
{
vtkm::Vec3f pc;
//See if point is inside the last cell.
if ((lastCell.CellId >= 0) && (lastCell.CellId < this->CellSet.GetNumberOfElements()) &&
this->PointInCell(point, lastCell.CellId, pc) == vtkm::ErrorCode::Success)
{
parametric = pc;
cellId = lastCell.CellId;
return vtkm::ErrorCode::Success;
}
//See if it's in the last leaf.
if ((lastCell.LeafIdx >= 0) && (lastCell.LeafIdx < this->CellCount.GetNumberOfValues()) &&
this->PointInLeaf(point, lastCell.LeafIdx, cellId, pc) == vtkm::ErrorCode::Success)
{
parametric = pc;
lastCell.CellId = cellId;
return vtkm::ErrorCode::Success;
}
//Call the full point search.
return this->FindCellImpl(point, cellId, parametric, lastCell);
}
private:
VTKM_EXEC
vtkm::ErrorCode PointInCell(const vtkm::Vec3f& point,
const vtkm::Id& cid,
vtkm::Vec3f& parametric) const
{
auto indices = this->CellSet.GetIndices(cid);
auto pts = vtkm::make_VecFromPortalPermute(&indices, this->Coords);
vtkm::Vec3f pc;
bool inside;
auto status = PointInsideCell(point, this->CellSet.GetCellShape(cid), pts, pc, inside);
if (status == vtkm::ErrorCode::Success && inside)
{
parametric = pc;
return vtkm::ErrorCode::Success;
}
return vtkm::ErrorCode::CellNotFound;
}
VTKM_EXEC
vtkm::ErrorCode PointInLeaf(const FloatVec3& point,
const vtkm::Id& leafIdx,
vtkm::Id& cellId,
FloatVec3& parametric) const
{
vtkm::Id start = this->CellStartIndex.Get(leafIdx);
vtkm::Id end = start + this->CellCount.Get(leafIdx);
for (vtkm::Id i = start; i < end; ++i)
{
vtkm::Vec3f pc;
vtkm::Id cid = this->CellIds.Get(i);
if (this->PointInCell(point, cid, pc) == vtkm::ErrorCode::Success)
{
cellId = cid;
parametric = pc;
return vtkm::ErrorCode::Success;
}
}
return vtkm::ErrorCode::CellNotFound;
}
VTKM_EXEC
vtkm::ErrorCode FindCellImpl(const FloatVec3& point,
vtkm::Id& cellId,
FloatVec3& parametric,
LastCell& lastCell) const
{
using namespace vtkm::internal::cl_uniform_bins;
cellId = -1;
lastCell.CellId = -1;
lastCell.LeafIdx = -1;
DimVec3 binId3 = static_cast<DimVec3>((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 vtkm::ErrorCode::CellNotFound;
}
auto leafGrid = ComputeLeafGrid(binId3, ldim, this->TopLevel);
DimVec3 leafId3 = static_cast<DimVec3>((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 leafIdx = leafStart + ComputeFlatIndex(leafId3, leafGrid.Dimensions);
if (this->PointInLeaf(point, leafIdx, cellId, parametric) == vtkm::ErrorCode::Success)
{
lastCell.CellId = cellId;
lastCell.LeafIdx = leafIdx;
return vtkm::ErrorCode::Success;
}
}
return vtkm::ErrorCode::CellNotFound;
}
vtkm::internal::cl_uniform_bins::Grid TopLevel;
ReadPortal<DimVec3> LeafDimensions;
ReadPortal<vtkm::Id> LeafStartIndex;
ReadPortal<vtkm::Id> CellStartIndex;
ReadPortal<vtkm::Id> CellCount;
ReadPortal<vtkm::Id> CellIds;
CellStructureType CellSet;
CoordsPortalType Coords;
};
}
} // vtkm::exec
#endif //vtk_m_exec_CellLocatorTwoLevel_h