Shuenhoy/vtk-m
1
0
Fork 0
mirror of https://gitlab.kitware.com/vtk/vtk-m synced 2024-09-08 13:23:51 +00:00
Code Issues 2 Actions Packages Projects Releases Wiki Activity

Add a CellLocator

Browse Source 
Implements a two-level uniform grid cell locator
This commit is contained in:
Sujin Philip 2017-09-28 16:45:51 -04:00
parent 02f48cfaaa
commit 41679cb5f9
12 changed files with 1197 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
vtkm/cont/CMakeLists.txt
View File

@ -44,6 +44,7 @@ set(headers
ArrayHandleConcatenate.h
ArrayRangeCompute.h
ArrayRangeCompute.hxx
CellLocatorTwoLevelUniformGrid.h
CellSet.h
CellSetExplicit.h
CellSetListTag.h

747
vtkm/cont/CellLocatorTwoLevelUniformGrid.h Normal file
View File

@ -0,0 +1,747 @@
//============================================================================
// 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 <vtkm/cont/ArrayHandle.h>
#include <vtkm/cont/ArrayHandleConstant.h>
#include <vtkm/cont/ArrayHandleTransform.h>
#include <vtkm/cont/CoordinateSystem.h>
#include <vtkm/cont/DeviceAdapterAlgorithm.h>
#include <vtkm/cont/DynamicCellSet.h>
#include <vtkm/exec/CellInside.h>
#include <vtkm/exec/ParametricCoordinates.h>
#include <vtkm/worklet/DispatcherMapField.h>
#include <vtkm/worklet/DispatcherMapTopology.h>
#include <vtkm/worklet/WorkletMapField.h>
#include <vtkm/worklet/WorkletMapTopology.h>
#include <vtkm/Math.h>
#include <vtkm/Types.h>
#include <vtkm/VecFromPortalPermute.h>
#include <vtkm/VecTraits.h>
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<DimensionType, 3>;
using FloatVec3 = vtkm::Vec<vtkm::FloatDefault, 3>;
struct BinsBBox
{
DimVec3 Min;
DimVec3 Max;
};
struct Bounds
{
FloatVec3 Min;
FloatVec3 Max;
};
struct Grid
{
DimVec3 Dimensions;
FloatVec3 Origin;
FloatVec3 BinSize;
};
struct TwoLevelUniformGrid
{
Grid TopLevel;
vtkm::cont::ArrayHandle<DimVec3> LeafDimensions;
vtkm::cont::ArrayHandle<vtkm::Id> LeafStartIndex;
vtkm::cont::ArrayHandle<vtkm::Id> CellStartIndex;
vtkm::cont::ArrayHandle<vtkm::Id> CellCount;
vtkm::cont::ArrayHandle<vtkm::Id> CellIds;
};
VTKM_EXEC_CONT static DimVec3 ComputeGridDimension(vtkm::Id numberOfCells,
const FloatVec3& size,
vtkm::FloatDefault density)
{
vtkm::FloatDefault nsides = 1.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<vtkm::FloatDefault>(numberOfCells) / (volume * density)),
1.0f / nsides);
return vtkm::Max(DimVec3(1),
DimVec3(static_cast<DimensionType>(size[0] * r),
static_cast<DimensionType>(size[1] * r),
static_cast<DimensionType>(size[2] * r)));
}
template <typename PointsVecType>
VTKM_EXEC static 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) };
}
// 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 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;
}
VTKM_EXEC static BinsBBox ComputeIntersectingBins(const Bounds cellBounds, const Grid& grid)
{
auto minb = (cellBounds.Min - grid.Origin) / grid.BinSize;
auto maxb = (cellBounds.Max - grid.Origin) / grid.BinSize;
return { vtkm::Max(DimVec3(0), DimVec3(minb)),
vtkm::Min(grid.Dimensions - DimVec3(1), DimVec3(maxb)) };
}
VTKM_EXEC static vtkm::Id GetNumberOfBins(const BinsBBox& binsBBox)
{
return (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)
, Step(1,
dim[0] - (bbox.Max[0] - bbox.Min[0] + 1),
(dim[0] * dim[1]) - ((bbox.Max[1] - bbox.Min[1] + 1) * dim[0]))
, FlatIdx(bbox.Min[0] + (dim[0] * (bbox.Min[1] + (dim[1] * bbox.Min[2]))))
{
}
VTKM_EXEC_CONT void Init()
{
this->Idx = this->BBox.Min;
this->FlatIdx =
this->Idx[0] + (this->Dim[0] * (this->Idx[1] + (this->Dim[1] * this->Idx[2])));
}
VTKM_EXEC_CONT bool Done() const { return this->Idx[2] > this->BBox.Max[2]; }
VTKM_EXEC_CONT void Next()
{
++this->Idx[0];
this->FlatIdx += this->Step[0];
if (this->Idx[0] > this->BBox.Max[0])
{
this->Idx[0] = this->BBox.Min[0];
++this->Idx[1];
this->FlatIdx += this->Step[1];
if (this->Idx[1] > this->BBox.Max[1])
{
this->Idx[1] = this->BBox.Min[1];
++this->Idx[2];
this->FlatIdx += this->Step[2];
}
}
}
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::Id3 Step;
vtkm::Id FlatIdx;
};
public:
class CountBinsL1 : public vtkm::worklet::WorkletMapPointToCell
{
public:
typedef void ControlSignature(CellSetIn cellset,
FieldInPoint<Vec3> coords,
FieldOutCell<IdType> bincount);
typedef void ExecutionSignature(_2, _3);
CountBinsL1(const Grid& grid)
: L1Grid(grid)
{
}
template <typename PointsVecType>
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:
typedef void ControlSignature(CellSetIn cellset,
FieldInPoint<Vec3> coords,
FieldInCell<IdType> offsets,
WholeArrayOut<IdType> binIds);
typedef void ExecutionSignature(_2, _3, _4);
FindBinsL1(const Grid& grid)
: L1Grid(grid)
{
}
template <typename PointsVecType, typename BinIdsPortalType>
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:
typedef void ControlSignature(FieldIn<IdType> binIds,
FieldIn<IdType> cellCounts,
WholeArrayOut<vtkm::ListTagBase<DimVec3>> dimensions);
typedef void ExecutionSignature(_1, _2, _3);
using InputDomain = _1;
GenerateBinsL1(FloatVec3 size, vtkm::FloatDefault density)
: Size(size)
, Density(density)
{
}
template <typename OutputDimensionsPortal>
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:
typedef void ControlSignature(CellSetIn cellset,
FieldInPoint<Vec3> coords,
WholeArrayIn<vtkm::ListTagBase<DimVec3>> binDimensions,
FieldOutCell<IdType> bincount);
typedef void ExecutionSignature(_2, _3, _4);
CountBinsL2(const Grid& grid)
: L1Grid(grid)
{
}
template <typename PointsVecType, typename BinDimensionsPortalType>
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;
leaf.Dimensions = binDimensions.Get(i.GetFlatIdx());
leaf.Origin =
this->L1Grid.Origin + (static_cast<FloatVec3>(i.GetIdx()) * this->L1Grid.BinSize);
leaf.BinSize = this->L1Grid.BinSize / static_cast<FloatVec3>(leaf.Dimensions);
auto binsBBoxL2 = ComputeIntersectingBins(cellBounds, leaf);
numBins += GetNumberOfBins(binsBBoxL2);
}
}
private:
Grid L1Grid;
};
class FindBinsL2 : public vtkm::worklet::WorkletMapPointToCell
{
public:
typedef void ControlSignature(CellSetIn cellset,
FieldInPoint<Vec3> coords,
WholeArrayIn<vtkm::ListTagBase<DimVec3>> binDimensions,
WholeArrayIn<IdType> binStarts,
FieldInCell<IdType> offsets,
WholeArrayOut<IdType> binIds,
WholeArrayOut<IdType> cellIds);
typedef void ExecutionSignature(InputIndex, _2, _3, _4, _5, _6, _7);
FindBinsL2(const Grid& grid)
: L1Grid(grid)
{
}
template <typename PointsVecType,
typename BinDimensionsPortalType,
typename BinStartsPortalType,
typename BinIdsPortalType,
typename CellIdsPortalType>
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;
leaf.Dimensions = binDimensions.Get(i.GetFlatIdx());
leaf.Origin =
this->L1Grid.Origin + (static_cast<FloatVec3>(i.GetIdx()) * this->L1Grid.BinSize);
leaf.BinSize = this->L1Grid.BinSize / static_cast<FloatVec3>(leaf.Dimensions);
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:
typedef void ControlSignature(FieldIn<IdType> binIds,
FieldIn<IdType> startsIn,
FieldIn<IdType> countsIn,
WholeArrayOut<IdType> startsOut,
WholeArrayOut<IdType> countsOut);
typedef void ExecutionSignature(_1, _2, _3, _4, _5);
using InputDomain = _1;
template <typename CellStartsPortalType, typename CellCountsPortalType>
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 <typename DeviceAdapter,
typename CellSetList = VTKM_DEFAULT_CELL_SET_LIST_TAG,
typename CoordsTypeList = VTKM_DEFAULT_COORDINATE_SYSTEM_TYPE_LIST_TAG,
typename CoordsStorageList = VTKM_DEFAULT_COORDINATE_SYSTEM_STORAGE_LIST_TAG>
void Build(DeviceAdapter,
CellSetList cellSetTypes = CellSetList(),
CoordsTypeList coordsValueTypes = CoordsTypeList(),
CoordsStorageList coordsStorageType = CoordsStorageList())
{
using Algorithm = vtkm::cont::DeviceAdapterAlgorithm<DeviceAdapter>;
auto cellset = this->CellSet.ResetCellSetList(cellSetTypes);
auto points =
this->Coordinates.GetData().ResetTypeAndStorageLists(coordsValueTypes, coordsStorageType);
TwoLevelUniformGrid ls;
// 1: Compute the top level grid
auto bounds = this->Coordinates.GetBounds(coordsValueTypes, coordsStorageType);
FloatVec3 bmin(static_cast<vtkm::FloatDefault>(bounds.X.Min),
static_cast<vtkm::FloatDefault>(bounds.Y.Min),
static_cast<vtkm::FloatDefault>(bounds.Z.Min));
FloatVec3 bmax(static_cast<vtkm::FloatDefault>(bounds.X.Max),
static_cast<vtkm::FloatDefault>(bounds.Y.Max),
static_cast<vtkm::FloatDefault>(bounds.Z.Max));
auto size = bmax - bmin;
auto fudge = vtkm::Max(FloatVec3(1e-6f), size * 1e-3f);
size += 2.0f * fudge;
ls.TopLevel.Dimensions =
ComputeGridDimension(cellset.GetNumberOfCells(), size, this->DensityL1);
ls.TopLevel.Origin = bmin - fudge;
ls.TopLevel.BinSize = size / static_cast<FloatVec3>(ls.TopLevel.Dimensions);
// 2: For each cell, find the number of top level bins they intersect
vtkm::cont::ArrayHandle<vtkm::Id> binCounts;
CountBinsL1 countL1(ls.TopLevel);
vtkm::worklet::DispatcherMapTopology<CountBinsL1, DeviceAdapter>(countL1).Invoke(
cellset, points, binCounts);
// 3: Total number of unique (cell, bin) pairs (for pre-allocating arrays)
vtkm::Id numPairsL1 = Algorithm::ScanExclusive(binCounts, binCounts);
// 4: For each cell find the top level bins that intersect it
vtkm::cont::ArrayHandle<vtkm::Id> binIds;
binIds.Allocate(numPairsL1);
FindBinsL1 findL1(ls.TopLevel);
vtkm::worklet::DispatcherMapTopology<FindBinsL1, DeviceAdapter>(findL1).Invoke(
cellset, points, binCounts, binIds);
binCounts.ReleaseResources();
// 5: From above, find the number of cells that intersect each top level bin
Algorithm::Sort(binIds);
vtkm::cont::ArrayHandle<vtkm::Id> bins;
vtkm::cont::ArrayHandle<vtkm::Id> cellsPerBin;
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];
Algorithm::Copy(vtkm::cont::make_ArrayHandleConstant(DimVec3(0), numberOfBins),
ls.LeafDimensions);
GenerateBinsL1 generateL1(ls.TopLevel.BinSize, this->DensityL2);
vtkm::worklet::DispatcherMapField<GenerateBinsL1, DeviceAdapter>(generateL1)
.Invoke(bins, cellsPerBin, ls.LeafDimensions);
bins.ReleaseResources();
cellsPerBin.ReleaseResources();
// 7: Compute number of level-2 bins
vtkm::Id numberOfLeaves = 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);
vtkm::worklet::DispatcherMapTopology<CountBinsL2, DeviceAdapter>(countL2).Invoke(
cellset, points, ls.LeafDimensions, binCounts);
// 9: Total number of unique (cell, bin) pairs (for pre-allocating arrays)
vtkm::Id numPairsL2 = 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);
vtkm::worklet::DispatcherMapTopology<FindBinsL2, DeviceAdapter>(findL2).Invoke(
cellset, points, ls.LeafDimensions, ls.LeafStartIndex, binCounts, binIds, ls.CellIds);
binCounts.ReleaseResources();
// 11: From above, find the cells that each l2 bin intersects
Algorithm::SortByKey(binIds, ls.CellIds);
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<vtkm::Id> cellsStart;
Algorithm::ScanExclusive(cellsPerBin, cellsStart);
Algorithm::Copy(vtkm::cont::ArrayHandleConstant<vtkm::Id>(0, numberOfLeaves),
ls.CellStartIndex);
Algorithm::Copy(vtkm::cont::ArrayHandleConstant<vtkm::Id>(0, numberOfLeaves), ls.CellCount);
vtkm::worklet::DispatcherMapField<GenerateBinsL2, DeviceAdapter>().Invoke(
bins, cellsStart, cellsPerBin, ls.CellStartIndex, ls.CellCount);
std::swap(this->LookupStructure, ls);
}
template <typename DeviceAdapter>
struct TwoLevelUniformGridExecution : public vtkm::exec::ExecutionObjectBase
{
template <typename T>
using ArrayPortalConst =
typename vtkm::cont::ArrayHandle<T>::template ExecutionTypes<DeviceAdapter>::PortalConst;
Grid TopLevel;
ArrayPortalConst<DimVec3> LeafDimensions;
ArrayPortalConst<vtkm::Id> LeafStartIndex;
ArrayPortalConst<vtkm::Id> CellStartIndex;
ArrayPortalConst<vtkm::Id> CellCount;
ArrayPortalConst<vtkm::Id> CellIds;
};
class FindCellWorklet : public vtkm::worklet::WorkletMapField
{
public:
typedef void ControlSignature(FieldIn<Vec3> points,
WholeCellSetIn<> cellSet,
WholeArrayIn<Vec3> coordinates,
ExecObject lookupStruct,
FieldOut<IdType> cellIds,
FieldOut<Vec3> parametricCoordinates);
typedef void ExecutionSignature(_1, _2, _3, _4, _5, _6);
using InputDomain = _1;
template <typename PointType,
typename CellSetType,
typename CoordsPortalType,
typename LookupStructureType>
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<FloatDefault>(point[0]),
static_cast<FloatDefault>(point[1]),
static_cast<FloatDefault>(point[2]));
const Grid& topLevelGrid = lookupStruct.TopLevel;
DimVec3 binId3 = static_cast<DimVec3>((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 = binId3[0] +
(topLevelGrid.Dimensions[0] * (binId3[1] + (topLevelGrid.Dimensions[1] * binId3[2])));
Grid leafGrid;
leafGrid.Dimensions = lookupStruct.LeafDimensions.Get(binId);
if (leafGrid.Dimensions[0] + leafGrid.Dimensions[1] + leafGrid.Dimensions[2] == 0)
{
return;
}
leafGrid.Origin =
topLevelGrid.Origin + static_cast<FloatVec3>(binId3) * topLevelGrid.BinSize;
leafGrid.BinSize = topLevelGrid.BinSize / static_cast<FloatVec3>(leafGrid.Dimensions);
vtkm::Id leafStart = lookupStruct.LeafStartIndex.Get(binId);
DimVec3 leafId3 = static_cast<DimVec3>((p - leafGrid.Origin) / leafGrid.BinSize);
vtkm::Id leafId = leafStart +
(leafId3[0] +
(leafGrid.Dimensions[0] * (leafId3[1] + (leafGrid.Dimensions[1] * leafId3[2]))));
VTKM_ASSERT(leafId3[0] >= 0 && leafId3[0] < leafGrid.Dimensions[0] && leafId3[1] >= 0 &&
leafId3[1] < leafGrid.Dimensions[1] && leafId3[2] >= 0 &&
leafId3[2] < leafGrid.Dimensions[2]);
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<decltype(indices), CoordsPortalType> 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 <typename PointComponentType,
typename PointStorageType,
typename DeviceAdapter,
typename CellSetList = VTKM_DEFAULT_CELL_SET_LIST_TAG,
typename CoordsTypeList = VTKM_DEFAULT_COORDINATE_SYSTEM_TYPE_LIST_TAG,
typename CoordsStorageList = VTKM_DEFAULT_COORDINATE_SYSTEM_STORAGE_LIST_TAG>
void FindCells(
const vtkm::cont::ArrayHandle<vtkm::Vec<PointComponentType, 3>, PointStorageType>& points,
vtkm::cont::ArrayHandle<vtkm::Id>& cellIds,
vtkm::cont::ArrayHandle<FloatVec3>& parametricCoords,
DeviceAdapter device,
CellSetList cellSetTypes = CellSetList(),
CoordsTypeList coordsValueTypes = CoordsTypeList(),
CoordsStorageList coordsStorageType = CoordsStorageList()) const
{
vtkm::worklet::DispatcherMapField<FindCellWorklet, DeviceAdapter>().Invoke(
points,
this->CellSet.ResetCellSetList(cellSetTypes),
this->Coordinates.GetData().ResetTypeAndStorageLists(coordsValueTypes, coordsStorageType),
this->PrepareForDevice(device),
cellIds,
parametricCoords);
}
private:
template <typename DeviceAdapter>
TwoLevelUniformGridExecution<DeviceAdapter> PrepareForDevice(DeviceAdapter device) const
{
TwoLevelUniformGridExecution<DeviceAdapter> deviceObject;
deviceObject.TopLevel = this->LookupStructure.TopLevel;
deviceObject.LeafDimensions = this->LookupStructure.LeafDimensions.PrepareForInput(device);
deviceObject.LeafStartIndex = this->LookupStructure.LeafStartIndex.PrepareForInput(device);
deviceObject.CellStartIndex = this->LookupStructure.CellStartIndex.PrepareForInput(device);
deviceObject.CellCount = this->LookupStructure.CellCount.PrepareForInput(device);
deviceObject.CellIds = this->LookupStructure.CellIds.PrepareForInput(device);
return deviceObject;
}
vtkm::FloatDefault DensityL1, DensityL2;
vtkm::cont::DynamicCellSet CellSet;
vtkm::cont::CoordinateSystem Coordinates;
TwoLevelUniformGrid LookupStructure;
};
}
}
#endif // vtk_m_cont_CellLocatorTwoLevelUniformGrid_h

1
vtkm/cont/cuda/testing/CMakeLists.txt
View File

@ -20,6 +20,7 @@
set(unit_tests
UnitTestCudaArrayHandleFancy.cu
UnitTestCudaCellLocatorTwoLevelUniformGrid.cu
UnitTestCudaComputeRange.cu
UnitTestCudaDataSetExplicit.cu
UnitTestCudaDataSetSingleType.cu

26
vtkm/cont/cuda/testing/UnitTestCudaCellLocatorTwoLevelUniformGrid.cu Normal file
View File

@ -0,0 +1,26 @@
//============================================================================
// 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.
//============================================================================
#include <vtkm/cont/testing/TestingCellLocatorTwoLevelUniformGrid.h>
int UnitTestCudaCellLocatorTwoLevelUniformGrid(int, char* [])
{
return vtkm::cont::testing::Testing::Run(
TestingCellLocatorTwoLevelUniformGrid<vtkm::cont::DeviceAdapterTagCuda>);
}

1
vtkm/cont/serial/testing/CMakeLists.txt
View File

@ -21,6 +21,7 @@
set(unit_tests
UnitTestSerialArrayHandle.cxx
UnitTestSerialArrayHandleFancy.cxx
UnitTestSerialCellLocatorTwoLevelUniformGrid.cxx
UnitTestSerialComputeRange.cxx
UnitTestSerialDataSetExplicit.cxx
UnitTestSerialDataSetSingleType.cxx

26
vtkm/cont/serial/testing/UnitTestSerialCellLocatorTwoLevelUniformGrid.cxx Normal file
View File

@ -0,0 +1,26 @@
//============================================================================
// 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.
//============================================================================
#include <vtkm/cont/testing/TestingCellLocatorTwoLevelUniformGrid.h>
int UnitTestSerialCellLocatorTwoLevelUniformGrid(int, char* [])
{
return vtkm::cont::testing::Testing::Run(
TestingCellLocatorTwoLevelUniformGrid<vtkm::cont::DeviceAdapterTagSerial>);
}

1
vtkm/cont/tbb/testing/CMakeLists.txt
View File

@ -21,6 +21,7 @@
set(unit_tests
UnitTestTBBArrayHandle.cxx
UnitTestTBBArrayHandleFancy.cxx
UnitTestTBBCellLocatorTwoLevelUniformGrid.cxx
UnitTestTBBComputeRange.cxx
UnitTestTBBDataSetExplicit.cxx
UnitTestTBBDataSetSingleType.cxx

26
vtkm/cont/tbb/testing/UnitTestTBBCellLocatorTwoLevelUniformGrid.cxx Normal file
View File

@ -0,0 +1,26 @@
//============================================================================
// 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.
//============================================================================
#include <vtkm/cont/testing/TestingCellLocatorTwoLevelUniformGrid.h>
int UnitTestTBBCellLocatorTwoLevelUniformGrid(int, char* [])
{
return vtkm::cont::testing::Testing::Run(
TestingCellLocatorTwoLevelUniformGrid<vtkm::cont::DeviceAdapterTagTBB>);
}

1
vtkm/cont/testing/CMakeLists.txt
View File

@ -23,6 +23,7 @@ set(headers
MakeTestDataSet.h
Testing.h
TestingArrayHandles.h
TestingCellLocatorTwoLevelUniformGrid.h
TestingComputeRange.h
TestingDeviceAdapter.h
TestingDataSetExplicit.h

248
vtkm/cont/testing/TestingCellLocatorTwoLevelUniformGrid.h Normal file
View File

@ -0,0 +1,248 @@
//============================================================================
// 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_testing_TestingCellLocatorTwoLevelUniformGrid_h
#define vtk_m_cont_testing_TestingCellLocatorTwoLevelUniformGrid_h
#include <vtkm/cont/CellLocatorTwoLevelUniformGrid.h>
#include <vtkm/cont/DataSetBuilderUniform.h>
#include <vtkm/cont/testing/Testing.h>
#include <vtkm/exec/CellInterpolate.h>
#include <vtkm/worklet/DispatcherMapTopology.h>
#include <vtkm/worklet/ScatterPermutation.h>
#include <vtkm/worklet/Tetrahedralize.h>
#include <vtkm/worklet/Triangulate.h>
#include <vtkm/worklet/WorkletMapTopology.h>
#include <vtkm/CellShape.h>
#include <ctime>
#include <random>
namespace
{
using PointType = vtkm::Vec<vtkm::FloatDefault, 3>;
std::default_random_engine RandomGenerator;
class ParametricToWorldCoordinates : public vtkm::worklet::WorkletMapPointToCell
{
public:
typedef void ControlSignature(CellSetIn cellset,
FieldInPoint<Vec3> coords,
FieldInOutCell<Vec3> pcs,
FieldOutCell<Vec3> wcs);
typedef void ExecutionSignature(CellShape, _2, _3, _4);
using ScatterType = vtkm::worklet::ScatterPermutation<>;
explicit ParametricToWorldCoordinates(const vtkm::cont::ArrayHandle<vtkm::Id>& cellIds)
: Scatter(cellIds)
{
}
const ScatterType& GetScatter() const { return this->Scatter; }
template <typename CellShapeTagType, typename PointsVecType>
VTKM_EXEC void operator()(CellShapeTagType cellShape,
PointsVecType points,
const PointType& pc,
PointType& wc) const
{
wc = vtkm::exec::CellInterpolate(points, pc, cellShape, *this);
}
private:
ScatterType Scatter;
};
template <vtkm::IdComponent DIMENSIONS, typename DeviceAdapter>
vtkm::cont::DataSet MakeTestDataSet(const vtkm::Vec<vtkm::Id, DIMENSIONS>& dims,
DeviceAdapter device)
{
using Algorithm = vtkm::cont::DeviceAdapterAlgorithm<DeviceAdapter>;
using Connectivity = vtkm::internal::ConnectivityStructuredInternals<DIMENSIONS>;
const vtkm::IdComponent PointsPerCell = 1 << DIMENSIONS;
auto uniformDs =
vtkm::cont::DataSetBuilderUniform::Create(dims,
vtkm::Vec<vtkm::FloatDefault, DIMENSIONS>(0.0f),
vtkm::Vec<vtkm::FloatDefault, DIMENSIONS>(1.0f));
// copy points
vtkm::cont::ArrayHandle<PointType> points;
Algorithm::Copy(uniformDs.GetCoordinateSystem()
.GetData()
.template Cast<vtkm::cont::ArrayHandleUniformPointCoordinates>(),
points);
vtkm::Id numberOfCells = uniformDs.GetCellSet().GetNumberOfCells();
vtkm::Id numberOfIndices = numberOfCells * PointsPerCell;
Connectivity structured;
structured.SetPointDimensions(dims);
// copy connectivity
vtkm::cont::ArrayHandle<vtkm::Id> connectivity;
connectivity.Allocate(numberOfIndices);
for (vtkm::Id i = 0, idx = 0; i < numberOfCells; ++i)
{
auto ptids = structured.GetPointsOfCell(i);
for (vtkm::IdComponent j = 0; j < PointsPerCell; ++j, ++idx)
{
connectivity.GetPortalControl().Set(idx, ptids[j]);
}
}
auto uniformCs =
uniformDs.GetCellSet().template Cast<vtkm::cont::CellSetStructured<DIMENSIONS>>();
vtkm::cont::CellSetSingleType<> cellset;
// triangulate the cellset
switch (DIMENSIONS)
{
case 2:
cellset = vtkm::worklet::Triangulate().Run(uniformCs, device);
break;
case 3:
cellset = vtkm::worklet::Tetrahedralize().Run(uniformCs, device);
break;
default:
VTKM_ASSERT(false);
}
// It is posible that the warping will result in invalid cells. So use a
// local random generator with a known seed that does not create invalid cells.
std::default_random_engine rgen;
// Warp the coordinates
std::uniform_real_distribution<vtkm::FloatDefault> warpFactor(-0.25f, 0.25f);
auto pointsPortal = points.GetPortalControl();
for (vtkm::Id i = 0; i < pointsPortal.GetNumberOfValues(); ++i)
{
PointType warpVec(0);
for (vtkm::IdComponent c = 0; c < DIMENSIONS; ++c)
{
warpVec[c] = warpFactor(rgen);
}
pointsPortal.Set(i, pointsPortal.Get(i) + warpVec);
}
// build dataset
vtkm::cont::DataSet out;
out.AddCoordinateSystem(vtkm::cont::CoordinateSystem("coords", points));
out.AddCellSet(cellset);
return out;
}
template <vtkm::IdComponent DIMENSIONS>
void GenerateRandomInput(const vtkm::cont::DataSet& ds,
vtkm::Id count,
vtkm::cont::ArrayHandle<vtkm::Id>& cellIds,
vtkm::cont::ArrayHandle<PointType>& pcoords,
vtkm::cont::ArrayHandle<PointType>& wcoords)
{
vtkm::Id numberOfCells = ds.GetCellSet().GetNumberOfCells();
std::uniform_int_distribution<vtkm::Id> cellIdGen(0, numberOfCells - 1);
cellIds.Allocate(count);
pcoords.Allocate(count);
wcoords.Allocate(count);
for (vtkm::Id i = 0; i < count; ++i)
{
cellIds.GetPortalControl().Set(i, cellIdGen(RandomGenerator));
PointType pc(0.0f);
vtkm::FloatDefault minPc = 1e-3f;
vtkm::FloatDefault sum = 0.0f;
for (vtkm::IdComponent c = 0; c < DIMENSIONS; ++c)
{
vtkm::FloatDefault maxPc =
1.0f - (static_cast<vtkm::FloatDefault>(DIMENSIONS - c) * minPc) - sum;
std::uniform_real_distribution<vtkm::FloatDefault> pcoordGen(minPc, maxPc);
pc[c] = pcoordGen(RandomGenerator);
sum += pc[c];
}
pcoords.GetPortalControl().Set(i, pc);
}
ParametricToWorldCoordinates pc2wc(cellIds);
vtkm::worklet::DispatcherMapTopology<ParametricToWorldCoordinates>(pc2wc).Invoke(
ds.GetCellSet(), ds.GetCoordinateSystem().GetData(), pcoords, wcoords);
}
template <vtkm::IdComponent DIMENSIONS, typename DeviceAdapter>
void TestCellLocator(const vtkm::Vec<vtkm::Id, DIMENSIONS>& dim,
vtkm::Id numberOfPoints,
DeviceAdapter device)
{
auto ds = MakeTestDataSet(dim, device);
std::cout << "Testing " << DIMENSIONS << "D dataset with " << ds.GetCellSet().GetNumberOfCells()
<< " cells\n";
vtkm::cont::CellLocatorTwoLevelUniformGrid locator;
locator.SetDensityL1(64.0f);
locator.SetDensityL2(1.0f);
locator.SetCellSet(ds.GetCellSet());
locator.SetCoordinates(ds.GetCoordinateSystem());
locator.Build(device);
vtkm::cont::ArrayHandle<vtkm::Id> expCellIds;
vtkm::cont::ArrayHandle<PointType> expPCoords;
vtkm::cont::ArrayHandle<PointType> points;
GenerateRandomInput<DIMENSIONS>(ds, numberOfPoints, expCellIds, expPCoords, points);
std::cout << "Finding cells for " << numberOfPoints << " points\n";
vtkm::cont::ArrayHandle<vtkm::Id> cellIds;
vtkm::cont::ArrayHandle<PointType> pcoords;
locator.FindCells(points, cellIds, pcoords, device);
for (vtkm::Id i = 0; i < numberOfPoints; ++i)
{
VTKM_TEST_ASSERT(cellIds.GetPortalConstControl().Get(i) ==
expCellIds.GetPortalConstControl().Get(i),
"Incorrect cell ids");
VTKM_TEST_ASSERT(test_equal(pcoords.GetPortalConstControl().Get(i),
expPCoords.GetPortalConstControl().Get(i),
1e-3),
"Incorrect parameteric coordinates");
}
}
} // anonymous
template <typename DeviceAdapter>
void TestingCellLocatorTwoLevelUniformGrid()
{
vtkm::UInt32 seed = static_cast<vtkm::UInt32>(std::time(nullptr));
std::cout << "Seed: " << seed << std::endl;
RandomGenerator.seed(seed);
TestCellLocator(vtkm::Id3(8), 512, DeviceAdapter()); // 3D dataset
TestCellLocator(vtkm::Id2(18), 512, DeviceAdapter()); // 2D dataset
}
#endif // vtk_m_cont_testing_TestingCellLocatorTwoLevelUniformGrid_h

1
vtkm/exec/CMakeLists.txt
View File

@ -23,6 +23,7 @@ set(headers
CellDerivative.h
CellEdge.h
CellFace.h
CellInside.h
CellInterpolate.h
ConnectivityExplicit.h
ConnectivityPermuted.h

118
vtkm/exec/CellInside.h Normal file
View File

@ -0,0 +1,118 @@
//============================================================================
// 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_exec_CellInside_h
#define vtk_m_exec_CellInside_h
#include <vtkm/CellShape.h>
#include <vtkm/Types.h>
namespace vtkm
{
namespace exec
{
template <typename T>
static inline VTKM_EXEC bool CellInside(const vtkm::Vec<T, 3>&, vtkm::CellShapeTagEmpty)
{
return false;
}
template <typename T>
static inline VTKM_EXEC bool CellInside(const vtkm::Vec<T, 3>& pcoords, vtkm::CellShapeTagVertex)
{
return pcoords[0] == T(0) && pcoords[1] == T(0) && pcoords[2] == T(0);
}
template <typename T>
static inline VTKM_EXEC bool CellInside(const vtkm::Vec<T, 3>& pcoords, vtkm::CellShapeTagLine)
{
return pcoords[0] >= T(0) && pcoords[0] <= T(1);
}
template <typename T>
static inline VTKM_EXEC bool CellInside(const vtkm::Vec<T, 3>& pcoords, vtkm::CellShapeTagTriangle)
{
return pcoords[0] >= T(0) && pcoords[1] >= T(0) && (pcoords[0] + pcoords[1] <= T(1));
}
template <typename T>
static inline VTKM_EXEC bool CellInside(const vtkm::Vec<T, 3>& pcoords, vtkm::CellShapeTagPolygon)
{
return ((pcoords[0] - T(0.5)) * (pcoords[0] - T(0.5))) +
((pcoords[1] - T(0.5)) * (pcoords[1] - T(0.5))) <=
T(0.25);
}
template <typename T>
static inline VTKM_EXEC bool CellInside(const vtkm::Vec<T, 3>& pcoords, vtkm::CellShapeTagQuad)
{
return pcoords[0] >= T(0) && pcoords[0] <= T(1) && pcoords[1] >= T(0) && pcoords[1] <= T(1);
}
template <typename T>
static inline VTKM_EXEC bool CellInside(const vtkm::Vec<T, 3>& pcoords, vtkm::CellShapeTagTetra)
{
return pcoords[0] >= T(0) && pcoords[1] >= T(0) && pcoords[2] >= T(0) &&
(pcoords[0] + pcoords[1] + pcoords[2] <= T(1));
}
template <typename T>
static inline VTKM_EXEC bool CellInside(const vtkm::Vec<T, 3>& pcoords,
vtkm::CellShapeTagHexahedron)
{
return pcoords[0] >= T(0) && pcoords[0] <= T(1) && pcoords[1] >= T(0) && pcoords[1] <= T(1) &&
pcoords[2] >= T(0) && pcoords[2] <= T(1);
}
template <typename T>
static inline VTKM_EXEC bool CellInside(const vtkm::Vec<T, 3>& pcoords, vtkm::CellShapeTagWedge)
{
return pcoords[0] >= T(0) && pcoords[1] >= T(0) && pcoords[2] >= T(0) && pcoords[2] <= T(1) &&
(pcoords[0] + pcoords[1] <= T(1));
}
template <typename T>
static inline VTKM_EXEC bool CellInside(const vtkm::Vec<T, 3>& pcoords, vtkm::CellShapeTagPyramid)
{
return pcoords[0] >= T(0) && pcoords[0] <= T(1) && pcoords[1] >= T(0) && pcoords[1] <= T(1) &&
pcoords[2] >= T(0) && pcoords[2] <= T(1);
}
/// Checks if the paramteric coordinates `pcoords` are on the inside for the
/// specified cell type.
///
template <typename T>
static inline VTKM_EXEC bool CellInside(const vtkm::Vec<T, 3>& pcoords,
vtkm::CellShapeTagGeneric shape)
{
bool result = false;
switch (shape.Id)
{
vtkmGenericCellShapeMacro(result = CellInside(pcoords, CellShapeTag()));
default:
break;
}
return result;
}
}
} // vtkm::exec
#endif // vtk_m_exec_CellInside_h