Fixes for uniform cell locator for points right on the boundary. Also added some tests for boundary points.

This commit is contained in:
Dave Pugmire 2019-02-01 08:43:54 -05:00
parent 71f000465b
commit 5a63b19c2b
3 changed files with 63 additions and 47 deletions

@ -51,18 +51,14 @@ public:
throw vtkm::cont::ErrorBadType("Cell set is not 3D structured type");
Bounds = coords.GetBounds();
vtkm::Vec<vtkm::Id, 3> celldims =
cellSet.Cast<StructuredType>().GetSchedulingRange(vtkm::TopologyElementTagCell());
CellDims = cellSet.Cast<StructuredType>().GetSchedulingRange(vtkm::TopologyElementTagCell());
RangeTransform[0] = static_cast<vtkm::FloatDefault>(celldims[0]) /
RangeTransform[0] = static_cast<vtkm::FloatDefault>(CellDims[0]) /
static_cast<vtkm::FloatDefault>(Bounds.X.Length());
RangeTransform[1] = static_cast<vtkm::FloatDefault>(celldims[1]) /
RangeTransform[1] = static_cast<vtkm::FloatDefault>(CellDims[1]) /
static_cast<vtkm::FloatDefault>(Bounds.Y.Length());
RangeTransform[2] = static_cast<vtkm::FloatDefault>(celldims[2]) /
RangeTransform[2] = static_cast<vtkm::FloatDefault>(CellDims[2]) /
static_cast<vtkm::FloatDefault>(Bounds.Z.Length());
PlaneSize = celldims[0] * celldims[1];
RowSize = celldims[0];
}
struct PrepareForExecutionFunctor
@ -76,8 +72,7 @@ public:
ExecutionType* execObject =
new ExecutionType(contLocator.Bounds,
contLocator.RangeTransform,
contLocator.PlaneSize,
contLocator.RowSize,
contLocator.CellDims,
contLocator.GetCellSet().template Cast<StructuredType>(),
contLocator.GetCoordinates().GetData(),
DeviceAdapter());
@ -105,8 +100,7 @@ private:
vtkm::Bounds Bounds;
vtkm::Vec<vtkm::FloatDefault, 3> RangeTransform;
vtkm::Id PlaneSize;
vtkm::Id RowSize;
vtkm::Vec<vtkm::Id, 3> CellDims;
mutable HandleType ExecHandle;
};
}

@ -37,9 +37,9 @@
class LocatorWorklet : public vtkm::worklet::WorkletMapField
{
public:
LocatorWorklet(vtkm::Bounds& bounds_, vtkm::Vec<vtkm::Id, 3>& dims_)
: bounds(bounds_)
, dims(dims_)
LocatorWorklet(vtkm::Bounds& bounds, vtkm::Vec<vtkm::Id, 3>& cellDims)
: Bounds(bounds)
, CellDims(cellDims)
{
}
@ -51,16 +51,24 @@ public:
template <typename PointType>
VTKM_EXEC vtkm::Id CalculateCellId(const PointType& point) const
{
if (!bounds.Contains(point))
if (!Bounds.Contains(point))
return -1;
vtkm::Vec<vtkm::Id, 3> logical;
logical[0] = static_cast<vtkm::Id>(
vtkm::Floor((point[0] / bounds.X.Length()) * static_cast<vtkm::FloatDefault>(dims[0] - 1)));
logical[1] = static_cast<vtkm::Id>(
vtkm::Floor((point[1] / bounds.Y.Length()) * static_cast<vtkm::FloatDefault>(dims[1] - 1)));
logical[2] = static_cast<vtkm::Id>(
vtkm::Floor((point[2] / bounds.Z.Length()) * static_cast<vtkm::FloatDefault>(dims[2] - 1)));
return logical[2] * (dims[0] - 1) * (dims[1] - 1) + logical[1] * (dims[0] - 1) + logical[0];
logical[0] = (point[0] == Bounds.X.Max)
? CellDims[0] - 1
: static_cast<vtkm::Id>(vtkm::Floor((point[0] / Bounds.X.Length()) *
static_cast<vtkm::FloatDefault>(CellDims[0])));
logical[1] = (point[1] == Bounds.Y.Max)
? CellDims[1] - 1
: static_cast<vtkm::Id>(vtkm::Floor((point[1] / Bounds.Y.Length()) *
static_cast<vtkm::FloatDefault>(CellDims[1])));
logical[2] = (point[2] == Bounds.Z.Max)
? CellDims[2] - 1
: static_cast<vtkm::Id>(vtkm::Floor((point[2] / Bounds.Z.Length()) *
static_cast<vtkm::FloatDefault>(CellDims[2])));
return logical[2] * CellDims[0] * CellDims[1] + logical[1] * CellDims[0] + logical[0];
}
template <typename PointType, typename LocatorType>
@ -76,8 +84,8 @@ public:
}
private:
vtkm::Bounds bounds;
vtkm::Vec<vtkm::Id, 3> dims;
vtkm::Bounds Bounds;
vtkm::Vec<vtkm::Id, 3> CellDims;
};
template <typename DeviceAdapter>
@ -98,9 +106,9 @@ public:
std::cout << "Z bounds : " << bounds.Z.Min << " to " << bounds.Z.Max << std::endl;
using StructuredType = vtkm::cont::CellSetStructured<3>;
vtkm::Vec<vtkm::Id, 3> dims =
cellSet.Cast<StructuredType>().GetSchedulingRange(vtkm::TopologyElementTagPoint());
std::cout << "Dimensions of dataset : " << dims << std::endl;
vtkm::Vec<vtkm::Id, 3> cellDims =
cellSet.Cast<StructuredType>().GetSchedulingRange(vtkm::TopologyElementTagCell());
std::cout << "Dimensions of dataset : " << cellDims << std::endl;
vtkm::cont::CellLocatorUniformGrid locator;
locator.SetCoordinates(coords);
@ -124,13 +132,29 @@ public:
PointType point = vtkm::make_Vec(outBounds(dre), outBounds(dre), outBounds(dre));
pointsVec.push_back(point);
}
std::uniform_real_distribution<vtkm::Float32> outBounds2(-1.0f, 0.0f);
for (size_t i = 0; i < 5; i++)
{
PointType point = vtkm::make_Vec(outBounds2(dre), outBounds2(dre), outBounds2(dre));
pointsVec.push_back(point);
}
// Add points right on the boundary.
pointsVec.push_back(vtkm::make_Vec(0, 0, 0));
pointsVec.push_back(vtkm::make_Vec(4, 4, 4));
pointsVec.push_back(vtkm::make_Vec(4, 0, 0));
pointsVec.push_back(vtkm::make_Vec(0, 4, 0));
pointsVec.push_back(vtkm::make_Vec(0, 0, 4));
pointsVec.push_back(vtkm::make_Vec(4, 4, 0));
pointsVec.push_back(vtkm::make_Vec(0, 4, 4));
pointsVec.push_back(vtkm::make_Vec(4, 0, 4));
vtkm::cont::ArrayHandle<PointType> points = vtkm::cont::make_ArrayHandle(pointsVec);
// Query the points using the locators.
vtkm::cont::ArrayHandle<vtkm::Id> cellIds;
vtkm::cont::ArrayHandle<PointType> parametric;
vtkm::cont::ArrayHandle<bool> match;
LocatorWorklet worklet(bounds, dims);
LocatorWorklet worklet(bounds, cellDims);
vtkm::worklet::DispatcherMapField<LocatorWorklet> dispatcher(worklet);
dispatcher.SetDevice(DeviceAdapter());
dispatcher.Invoke(points, locator, cellIds, parametric, match);

@ -53,26 +53,18 @@ public:
VTKM_CONT
CellLocatorUniformGrid(const vtkm::Bounds& bounds,
const vtkm::Vec<vtkm::FloatDefault, 3> rangeTransform,
const vtkm::Id planeSize,
const vtkm::Id rowSize,
const vtkm::Vec<vtkm::Id, 3> cellDims,
const vtkm::cont::CellSetStructured<3>& cellSet,
const vtkm::cont::ArrayHandleVirtualCoordinates& coords,
DeviceAdapter)
: Bounds(bounds)
, RangeTransform(rangeTransform)
, PlaneSize(planeSize)
, RowSize(rowSize)
, CellDims(cellDims)
, PlaneSize(cellDims[0] * cellDims[1])
, RowSize(cellDims[0])
, CellSet(cellSet.PrepareForInput(DeviceAdapter(), FromType(), ToType()))
, Coords(coords.PrepareForInput(DeviceAdapter()))
{
CellSet = cellSet.PrepareForInput(DeviceAdapter(), FromType(), ToType());
Coords = coords.PrepareForInput(DeviceAdapter());
vtkm::Id3 cellDims = cellSet.GetCellDimensions();
this->Scale[0] = static_cast<vtkm::FloatDefault>(cellDims[0] - 1) /
static_cast<vtkm::FloatDefault>(Bounds.X.Length());
this->Scale[1] = static_cast<vtkm::FloatDefault>(cellDims[1] - 1) /
static_cast<vtkm::FloatDefault>(Bounds.Y.Length());
this->Scale[2] = static_cast<vtkm::FloatDefault>(cellDims[2] - 1) /
static_cast<vtkm::FloatDefault>(Bounds.Z.Length());
}
VTKM_EXEC
@ -88,9 +80,15 @@ public:
}
// Get the Cell Id from the point.
vtkm::Vec<vtkm::Id, 3> logicalCell;
logicalCell[0] = static_cast<vtkm::Id>(vtkm::Floor((point[0] - Bounds.X.Min) * Scale[0]));
logicalCell[1] = static_cast<vtkm::Id>(vtkm::Floor((point[1] - Bounds.Y.Min) * Scale[1]));
logicalCell[2] = static_cast<vtkm::Id>(vtkm::Floor((point[2] - Bounds.Z.Min) * Scale[2]));
logicalCell[0] = (point[0] == Bounds.X.Max)
? CellDims[0] - 1
: static_cast<vtkm::Id>(vtkm::Floor((point[0] - Bounds.X.Min) * RangeTransform[0]));
logicalCell[1] = (point[1] == Bounds.Y.Max)
? CellDims[1] - 1
: static_cast<vtkm::Id>(vtkm::Floor((point[1] - Bounds.Y.Min) * RangeTransform[1]));
logicalCell[2] = (point[2] == Bounds.Z.Max)
? CellDims[2] - 1
: static_cast<vtkm::Id>(vtkm::Floor((point[2] - Bounds.Z.Min) * RangeTransform[2]));
// Get the actual cellId, from the logical cell index of the cell
cellId = logicalCell[2] * PlaneSize + logicalCell[1] * RowSize + logicalCell[0];
@ -108,11 +106,11 @@ public:
private:
vtkm::Bounds Bounds;
vtkm::Vec<vtkm::FloatDefault, 3> RangeTransform;
vtkm::Vec<vtkm::Id, 3> CellDims;
vtkm::Id PlaneSize;
vtkm::Id RowSize;
CellSetPortal CellSet;
CoordsPortal Coords;
vtkm::Vec<vtkm::FloatDefault, 3> Scale;
};
}
}