Commit 67198714 authored by ayenpure's avatar ayenpure

Adding working implementation of BIH with test

- Uses smart pointers
  - need to get rid of them
parent 61fdfac7
......@@ -35,6 +35,7 @@
#include <vtkm/cont/CellLocator.h>
#include <vtkm/cont/DeviceAdapterAlgorithm.h>
#include <vtkm/cont/Timer.h>
#include <vtkm/cont/internal/DeviceAdapterListHelpers.h>
#include <vtkm/worklet/DispatcherMapField.h>
#include <vtkm/worklet/DispatcherMapTopology.h>
#include <vtkm/worklet/WorkletMapField.h>
......@@ -869,6 +870,69 @@ private:
}
};
struct PrepareForExecutionFunctor
{
public:
template <typename DeviceAdapter>
VTKM_CONT void operator()(DeviceAdapter,
const vtkm::cont::BoundingIntervalHierarchy* bih,
std::shared_ptr<vtkm::exec::CellLocator>& bihExec) const
{
vtkm::cont::DynamicCellSet cellSet = bih->GetCellSet();
if (cellSet.IsType<vtkm::cont::CellSetExplicit<>>())
{
using CellSetType = vtkm::cont::CellSetExplicit<>;
using ExecutionType = vtkm::exec::BoundingIntervalHierarchyExec<DeviceAdapter, CellSetType>;
bihExec = make_shared<ExecutionType>(bih->Nodes,
bih->ProcessedCellIds,
cellSet.Cast<CellSetType>(),
bih->GetCoords().GetData(),
DeviceAdapter());
}
else if (cellSet.IsType<vtkm::cont::CellSetStructured<2>>())
{
using CellSetType = vtkm::cont::CellSetStructured<2>;
using ExecutionType = vtkm::exec::BoundingIntervalHierarchyExec<DeviceAdapter, CellSetType>;
bihExec = make_shared<ExecutionType>(bih->Nodes,
bih->ProcessedCellIds,
cellSet.Cast<CellSetType>(),
bih->GetCoords().GetData(),
DeviceAdapter());
}
else if (cellSet.IsType<vtkm::cont::CellSetStructured<3>>())
{
using CellSetType = vtkm::cont::CellSetStructured<3>;
using ExecutionType = vtkm::exec::BoundingIntervalHierarchyExec<DeviceAdapter, CellSetType>;
bihExec = make_shared<ExecutionType>(bih->Nodes,
bih->ProcessedCellIds,
cellSet.Cast<CellSetType>(),
bih->GetCoords().GetData(),
DeviceAdapter());
}
else if (cellSet.IsType<vtkm::cont::CellSetSingleType<>>())
{
using CellSetType = vtkm::cont::CellSetSingleType<>;
using ExecutionType = vtkm::exec::BoundingIntervalHierarchyExec<DeviceAdapter, CellSetType>;
bihExec = make_shared<ExecutionType>(bih->Nodes,
bih->ProcessedCellIds,
cellSet.Cast<CellSetType>(),
bih->GetCoords().GetData(),
DeviceAdapter());
}
else
{
throw vtkm::cont::ErrorBadType("Could not determine type to write out.");
}
}
template <typename T, typename... Ts>
std::shared_ptr<T> make_shared(Ts&&... args) const
{
return std::shared_ptr<T>(new T(std::forward<Ts>(args)...));
}
};
public:
VTKM_CONT
BoundingIntervalHierarchy(vtkm::IdComponent numPlanes = 4, vtkm::IdComponent maxLeafSize = 5)
......@@ -884,11 +948,20 @@ public:
vtkm::cont::TryExecute(functor);
}
protected:
VTKM_CONT
std::unique_ptr<vtkm::exec::CellLocator> PrepareForExecution(
vtkm::cont::DeviceAdapterId device) override
std::shared_ptr<vtkm::exec::CellLocator> PrepareForExecutionOnDevice(
vtkm::cont::DeviceAdapterId& device) const override
{
using DeviceList = vtkm::ListTagBase<vtkm::cont::DeviceAdapterTagCuda,
vtkm::cont::DeviceAdapterTagTBB,
vtkm::cont::DeviceAdapterTagSerial>;
//vtkm::exec::CellLocator* toReturn;
std::shared_ptr<vtkm::exec::CellLocator> toReturn;
vtkm::cont::internal::FindDeviceAdapterTagAndCall(
device, DeviceList(), PrepareForExecutionFunctor(), this, toReturn);
//return make_unique<vtkm::exec::CellLocator>(*toReturn);
return toReturn;
}
private:
......
......@@ -38,7 +38,7 @@ using CellIdArrayHandle = vtkm::cont::ArrayHandle<vtkm::Id>;
} // namespace
template <typename DeviceAdapter>
template <typename DeviceAdapter, typename CellSetType>
class BoundingIntervalHierarchyExec : public vtkm::exec::CellLocator
{
public:
......@@ -48,33 +48,35 @@ public:
VTKM_CONT
BoundingIntervalHierarchyExec(const NodeArrayHandle& nodes,
const CellIdArrayHandle& cellIds,
const CellSetType& cellSet,
const vtkm::cont::ArrayHandleVirtualCoordinates& coords,
DeviceAdapter)
: Nodes(nodes.PrepareForInput(DeviceAdapter()))
, CellIds(cellIds.PrepareForInput(DeviceAdapter()))
{
CellSet = cellSet.PrepareForInput(DeviceAdapter(), FromType(), ToType());
Coords = coords.PrepareForInput(DeviceAdapter());
}
template <typename CellSetType, typename PointPortal>
VTKM_EXEC vtkm::Id Find(const vtkm::Vec<vtkm::Float64, 3>& point,
const CellSetType& cellSet,
const PointPortal& points,
const vtkm::exec::FunctorBase& worklet) const
VTKM_EXEC
void FindCell(const vtkm::Vec<vtkm::FloatDefault, 3>& point,
vtkm::Id& cellId,
vtkm::Vec<vtkm::FloatDefault, 3>& parametric,
const vtkm::exec::FunctorBase& worklet) const override
{
return Find(0, point, cellSet, points, worklet);
cellId = Find(0, point, parametric, worklet);
}
private:
template <typename CellSetType, typename PointPortal>
VTKM_EXEC vtkm::Id Find(vtkm::Id index,
const vtkm::Vec<vtkm::Float64, 3>& point,
const CellSetType& cellSet,
const PointPortal& points,
const vtkm::Vec<vtkm::FloatDefault, 3>& point,
vtkm::Vec<vtkm::FloatDefault, 3>& parametric,
const vtkm::exec::FunctorBase& worklet) const
{
const vtkm::cont::BoundingIntervalHierarchyNode& node = Nodes.Get(index);
if (node.ChildIndex < 0)
{
return FindInLeaf(point, node, cellSet, points, worklet);
return FindInLeaf(point, parametric, node, worklet);
}
else
{
......@@ -83,11 +85,11 @@ private:
vtkm::Id id2 = -1;
if (c <= node.Node.LMax)
{
id1 = Find(node.ChildIndex, point, cellSet, points, worklet);
id1 = Find(node.ChildIndex, point, parametric, worklet);
}
if (id1 == -1 && c >= node.Node.RMin)
{
id2 = Find(node.ChildIndex + 1, point, cellSet, points, worklet);
id2 = Find(node.ChildIndex + 1, point, parametric, worklet);
}
if (id1 == -1 && id2 == -1)
{
......@@ -104,20 +106,18 @@ private:
}
}
template <typename CellSetType, typename PointPortal>
VTKM_EXEC vtkm::Id FindInLeaf(const vtkm::Vec<vtkm::Float64, 3>& point,
VTKM_EXEC vtkm::Id FindInLeaf(const vtkm::Vec<vtkm::FloatDefault, 3>& point,
vtkm::Vec<vtkm::FloatDefault, 3>& parametric,
const vtkm::cont::BoundingIntervalHierarchyNode& node,
const CellSetType& cellSet,
const PointPortal& points,
const vtkm::exec::FunctorBase& worklet) const
{
using IndicesType = typename CellSetType::IndicesType;
using IndicesType = typename CellSetPortal::IndicesType;
for (vtkm::Id i = node.Leaf.Start; i < node.Leaf.Start + node.Leaf.Size; ++i)
{
vtkm::Id cellId = CellIds.Get(i);
IndicesType cellPointIndices = cellSet.GetIndices(cellId);
vtkm::VecFromPortalPermute<IndicesType, PointPortal> cellPoints(&cellPointIndices, points);
if (IsPointInCell(point, cellSet.GetCellShape(cellId), cellPoints, worklet))
IndicesType cellPointIndices = CellSet.GetIndices(cellId);
vtkm::VecFromPortalPermute<IndicesType, CoordsPortal> cellPoints(&cellPointIndices, Coords);
if (IsPointInCell(point, parametric, CellSet.GetCellShape(cellId), cellPoints, worklet))
{
return cellId;
}
......@@ -126,24 +126,32 @@ private:
}
template <typename CoordsType, typename CellShapeTag>
VTKM_EXEC static bool IsPointInCell(const vtkm::Vec<vtkm::Float64, 3>& point,
VTKM_EXEC static bool IsPointInCell(const vtkm::Vec<vtkm::FloatDefault, 3>& point,
vtkm::Vec<vtkm::FloatDefault, 3>& parametric,
CellShapeTag cellShape,
const CoordsType& cellPoints,
const vtkm::exec::FunctorBase& worklet)
{
bool success = false;
vtkm::Vec<vtkm::Float64, 3> parametricCoords =
vtkm::exec::WorldCoordinatesToParametricCoordinates(
cellPoints, point, cellShape, success, worklet);
return success && vtkm::exec::CellInside(parametricCoords, cellShape);
parametric = vtkm::exec::WorldCoordinatesToParametricCoordinates(
cellPoints, point, cellShape, success, worklet);
return success && vtkm::exec::CellInside(parametric, cellShape);
}
using FromType = vtkm::TopologyElementTagPoint;
using ToType = vtkm::TopologyElementTagCell;
using NodePortal = typename NodeArrayHandle::template ExecutionTypes<DeviceAdapter>::PortalConst;
using CellIdPortal =
typename CellIdArrayHandle::template ExecutionTypes<DeviceAdapter>::PortalConst;
using CellSetPortal =
typename CellSetType::template ExecutionTypes<DeviceAdapter, FromType, ToType>::ExecObjectType;
using CoordsPortal = typename vtkm::cont::ArrayHandleVirtualCoordinates::template ExecutionTypes<
DeviceAdapter>::PortalConst;
NodePortal Nodes;
CellIdPortal CellIds;
CellSetPortal CellSet;
CoordsPortal Coords;
}; // class BoundingIntervalHierarchyExec
} // namespace exec
......
......@@ -20,11 +20,11 @@ namespace exec
class CellLocator
{
public:
VTKM_EXEC void FindCell(const vtkm::Vec<vtkm::FloatDefault, 3>& point,
vtkm::Id& cellId,
vtkm::Vec<vtkm::FloatDefault, 3>& parametric)
{
} //const = 0;
VTKM_EXEC
virtual void FindCell(const vtkm::Vec<vtkm::FloatDefault, 3>& point,
vtkm::Id& cellId,
vtkm::Vec<vtkm::FloatDefault, 3>& parametric,
const vtkm::exec::FunctorBase& worklet) const = 0;
};
} // namespace exec
......@@ -68,14 +68,15 @@ public:
}
template <typename DeviceAdapter>
VTKM_CONT std::unique_ptr<vtkm::exec::CellLocator> PrepareForExecution(DeviceAdapter device)
VTKM_CONT std::shared_ptr<vtkm::exec::CellLocator> PrepareForExecution(DeviceAdapter) const
{
vtkm::cont::DeviceAdapterId deviceId = vtkm::cont::DeviceAdapterTraits<DeviceAdapter>::GetId();
return PrepareForExecution(deviceId);
return PrepareForExecutionOnDevice(deviceId);
}
VTKM_CONT virtual std::unique_ptr<vtkm::exec::CellLocator> PrepareForExecution(
vtkm::cont::DeviceAdapterId device) = 0;
protected:
VTKM_CONT virtual std::shared_ptr<vtkm::exec::CellLocator> PrepareForExecutionOnDevice(
vtkm::cont::DeviceAdapterId& device) const = 0;
private:
vtkm::cont::DynamicCellSet CellSet;
......
......@@ -55,6 +55,7 @@ set(unit_tests
UnitTestArrayHandleUniformPointCoordinates.cxx
UnitTestArrayHandleConcatenate.cxx
UnitTestArrayPortalToIterators.cxx
UnitTestBoundingIntervalHierarchy.cxx
UnitTestCellLocator.cxx
UnitTestCellSetExplicit.cxx
UnitTestCellSetPermutation.cxx
......
......@@ -19,12 +19,12 @@
//============================================================================
#include <vtkm/cont/ArrayHandleConcatenate.h>
#include <vtkm/cont/BoundingIntervalHierarchy.h>
#include <vtkm/cont/BoundingIntervalHierarchyExec.h>
#include <vtkm/cont/DataSetBuilderUniform.h>
#include <vtkm/cont/Timer.h>
#include <vtkm/cont/testing/Testing.h>
#include <vtkm/io/reader/VTKDataSetReader.h>
#include <vtkm/worklet/spatialstructure/BoundingIntervalHierarchy.h>
#include <vtkm/worklet/spatialstructure/BoundingIntervalHierarchyBuilder.h>
namespace
{
......@@ -47,25 +47,17 @@ struct CellCentroidCalculator : public vtkm::worklet::WorkletMapPointToCell
struct BoundingIntervalHierarchyTester : public vtkm::worklet::WorkletMapField
{
typedef void ControlSignature(FieldIn<>,
ExecObject,
WholeCellSetIn<>,
WholeArrayIn<>,
FieldIn<>,
FieldOut<>);
typedef _6 ExecutionSignature(_1, _2, _3, _4, _5);
template <typename Point,
typename BoundingIntervalHierarchyExecObject,
typename CellSet,
typename CoordsPortal>
typedef void ControlSignature(FieldIn<>, ExecObject, FieldIn<>, FieldOut<>);
typedef _4 ExecutionSignature(_1, _2, _3);
template <typename Point, typename BoundingIntervalHierarchyExecObject>
VTKM_EXEC vtkm::IdComponent operator()(const Point& point,
const BoundingIntervalHierarchyExecObject& bih,
const CellSet& cellSet,
const CoordsPortal& coords,
const vtkm::Id expectedId) const
{
vtkm::Id cellId = bih.Find(point, cellSet, coords, *this);
vtkm::Vec<vtkm::FloatDefault, 3> parametric;
vtkm::Id cellId;
bih->FindCell(point, cellId, parametric, *this);
return (1 - static_cast<vtkm::IdComponent>(expectedId == cellId));
}
}; // struct BoundingIntervalHierarchyTester
......@@ -80,14 +72,15 @@ void TestBoundingIntervalHierarchy(vtkm::cont::DataSet dataSet, vtkm::IdComponen
using DeviceAdapter = VTKM_DEFAULT_DEVICE_ADAPTER_TAG;
using Algorithms = vtkm::cont::DeviceAdapterAlgorithm<DeviceAdapter>;
using Timer = vtkm::cont::Timer<DeviceAdapter>;
namespace spatial = vtkm::worklet::spatialstructure;
vtkm::cont::DynamicCellSet cellSet = dataSet.GetCellSet();
vtkm::cont::DynamicArrayHandleCoordinateSystem vertices = dataSet.GetCoordinateSystem().GetData();
vtkm::cont::ArrayHandleVirtualCoordinates vertices = dataSet.GetCoordinateSystem().GetData();
std::cout << "Using numPlanes: " << numPlanes << "\n";
spatial::BoundingIntervalHierarchy bih = spatial::BoundingIntervalHierarchyBuilder(numPlanes, 5)
.Build(cellSet, vertices, DeviceAdapter());
vtkm::cont::BoundingIntervalHierarchy bih = vtkm::cont::BoundingIntervalHierarchy(numPlanes, 5);
bih.SetCellSet(cellSet);
bih.SetCoords(dataSet.GetCoordinateSystem());
bih.Build();
Timer centroidsTimer;
vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::FloatDefault, 3>> centroids;
......@@ -95,9 +88,6 @@ void TestBoundingIntervalHierarchy(vtkm::cont::DataSet dataSet, vtkm::IdComponen
cellSet, vertices, centroids);
std::cout << "Centroids calculation time: " << centroidsTimer.GetElapsedTime() << "\n";
vtkm::worklet::spatialstructure::BoundingIntervalHierarchyExecutionObject<DeviceAdapter>
bihExecObject = bih.PrepareForInput<DeviceAdapter>();
vtkm::cont::ArrayHandleCounting<vtkm::Id> expectedCellIds(0, 1, cellSet.GetNumberOfCells());
Timer interpolationTimer;
......@@ -108,8 +98,15 @@ void TestBoundingIntervalHierarchy(vtkm::cont::DataSet dataSet, vtkm::IdComponen
cudaDeviceGetLimit(&stackSizeBackup, cudaLimitStackSize);
cudaDeviceSetLimit(cudaLimitStackSize, 1024 * 200);
#endif
/*vtkm::cont::DeviceAdapterId deviceId = vtkm::cont::DeviceAdapterTraits<DeviceAdapter>::GetId();
std::unique_ptr<vtkm::exec::CellLocator> temp = bih.PrepareForExecution(deviceId);
vtkm::exec::BoundingIntervalHierarchyExec<DeviceAdapter>& bihExec
= dynamic_cast<vtkm::exec::BoundingIntervalHierarchyExec<DeviceAdapter>&>(*temp);*/
vtkm::worklet::DispatcherMapField<BoundingIntervalHierarchyTester>().Invoke(
centroids, bihExecObject, cellSet, vertices, expectedCellIds, results);
centroids, bih, expectedCellIds, results);
#if VTKM_DEVICE_ADAPTER == VTKM_DEVICE_ADAPTER_CUDA
cudaDeviceSetLimit(cudaLimitStackSize, stackSizeBackup);
#endif
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment