From a1e8d029c552543bda69c31091a6790bdbf237d1 Mon Sep 17 00:00:00 2001 From: Kenneth Moreland Date: Sat, 27 May 2023 15:22:49 -0600 Subject: [PATCH] Get the 3D index from a BoundaryState in WorkletPointNeighborhood There are occasions when you need a worklet to opeate on 2D or 3D indices. Most worklets operate on 1D indices, which requires recomputing the 3D index in each worklet instance. A workaround is to use a worklet that does a 3D scheduling and pull the working index from that. The problem was that there was no easy way to get this 3D index. To provide this option, a feature was added to the `BoundaryState` class that can be provided by `WorkletPointNeighborhood`. Thus, to get a 3D index in a worklet, use the `WorkletPointNeighborhood`, add `Boundary` as an argument to the `ExecutionSignature`, and then call `GetCenterIndex` on the `BoundaryState` object passed to the worklet operator. --- docs/changelog/index-from-boundary-state.md | 15 ++++ vtkm/exec/BoundaryState.h | 5 ++ .../UnitTestWorkletMapPointNeighborhood.cxx | 72 ++++++++++++++++++- 3 files changed, 91 insertions(+), 1 deletion(-) create mode 100644 docs/changelog/index-from-boundary-state.md diff --git a/docs/changelog/index-from-boundary-state.md b/docs/changelog/index-from-boundary-state.md new file mode 100644 index 000000000..ac3ee116e --- /dev/null +++ b/docs/changelog/index-from-boundary-state.md @@ -0,0 +1,15 @@ +# Get the 3D index from a BoundaryState in WorkletPointNeighborhood + +There are occasions when you need a worklet to opeate on 2D or 3D indices. +Most worklets operate on 1D indices, which requires recomputing the 3D +index in each worklet instance. A workaround is to use a worklet that does +a 3D scheduling and pull the working index from that. + +The problem was that there was no easy way to get this 3D index. To provide +this option, a feature was added to the `BoundaryState` class that can be +provided by `WorkletPointNeighborhood`. + +Thus, to get a 3D index in a worklet, use the `WorkletPointNeighborhood`, +add `Boundary` as an argument to the `ExecutionSignature`, and then call +`GetCenterIndex` on the `BoundaryState` object passed to the worklet +operator. diff --git a/vtkm/exec/BoundaryState.h b/vtkm/exec/BoundaryState.h index ebe2a50ce..27cbe037e 100644 --- a/vtkm/exec/BoundaryState.h +++ b/vtkm/exec/BoundaryState.h @@ -37,6 +37,11 @@ struct BoundaryState { } + /// Returns the center index of the neighborhood. This is typically the position of the + /// invocation of the worklet given this boundary condition. + /// + VTKM_EXEC const vtkm::Id3& GetCenterIndex() const { return this->IJK; } + ///@{ /// Returns true if a neighborhood of the given radius is contained within the bounds of the cell /// set in the X, Y, or Z direction. Returns false if the neighborhood extends outside of the diff --git a/vtkm/worklet/testing/UnitTestWorkletMapPointNeighborhood.cxx b/vtkm/worklet/testing/UnitTestWorkletMapPointNeighborhood.cxx index b74aa9f84..c7ad22d99 100644 --- a/vtkm/worklet/testing/UnitTestWorkletMapPointNeighborhood.cxx +++ b/vtkm/worklet/testing/UnitTestWorkletMapPointNeighborhood.cxx @@ -17,8 +17,11 @@ #include #include +#include +#include #include #include +#include #include #include @@ -178,7 +181,34 @@ struct ScatterUniformNeighbor : public vtkm::worklet::WorkletPointNeighborhood using ScatterType = vtkm::worklet::ScatterUniform<3>; }; -} + +// An example of using WorkletPointNeighborhood to iterate over a structured 3D cell +// domain rather than look at an actual neighborhood. It reduces a domain by subsampling +// every other item in the input field. +struct Subsample : public vtkm::worklet::WorkletPointNeighborhood +{ + using ControlSignature = + void(WholeCellSetIn inputTopology, + WholeArrayIn inputField, + CellSetIn outputTopology, + FieldOut sampledField); + using ExecutionSignature = void(_1, _2, Boundary, _4); + using InputDomain = _3; + + template + VTKM_EXEC void operator()(const vtkm::exec::ConnectivityStructured& inputTopology, + const InFieldPortal& inFieldPortal, + const vtkm::exec::BoundaryState& boundary, + T& sample) const + { + sample = + inFieldPortal.Get(inputTopology.LogicalToFlatVisitIndex(2 * boundary.GetCenterIndex())); + } +}; + +} // namespace test_pointneighborhood namespace { @@ -186,6 +216,7 @@ namespace static void TestMaxNeighborValue(); static void TestScatterIdentityNeighbor(); static void TestScatterUnfiormNeighbor(); +static void TestIndexing(); void TestWorkletPointNeighborhood(vtkm::cont::DeviceAdapterId id) { @@ -195,6 +226,7 @@ void TestWorkletPointNeighborhood(vtkm::cont::DeviceAdapterId id) TestMaxNeighborValue(); TestScatterIdentityNeighbor(); TestScatterUnfiormNeighbor(); + TestIndexing(); } static void TestMaxNeighborValue() @@ -276,6 +308,44 @@ static void TestScatterUnfiormNeighbor() dispatcher.Invoke(dataSet2D.GetCellSet(), dataSet2D.GetCoordinateSystem()); } +static void TestIndexing() +{ + std::cout << "Testing using PointNeighborhood for 3D indexing." << std::endl; + + constexpr vtkm::Id outDim = 4; + constexpr vtkm::Id inDim = outDim * 2; + + vtkm::cont::CellSetStructured<3> inCellSet; + inCellSet.SetPointDimensions({ inDim }); + vtkm::cont::CellSetStructured<3> outCellSet; + outCellSet.SetPointDimensions({ outDim }); + + vtkm::cont::ArrayHandleUniformPointCoordinates inField{ { inDim } }; + + vtkm::cont::ArrayHandle outField; + + vtkm::cont::Invoker invoke; + invoke(::test_pointneighborhood::Subsample{}, inCellSet, inField, outCellSet, outField); + + VTKM_TEST_ASSERT(outField.GetNumberOfValues() == (outDim * outDim * outDim)); + + vtkm::Id flatIndex = 0; + vtkm::Id3 IJK; + auto outFieldPortal = outField.WritePortal(); + for (IJK[2] = 0; IJK[2] < outDim; ++IJK[2]) + { + for (IJK[1] = 0; IJK[1] < outDim; ++IJK[1]) + { + for (IJK[0] = 0; IJK[0] < outDim; ++IJK[0]) + { + vtkm::Vec3f computed = outFieldPortal.Get(flatIndex); + VTKM_TEST_ASSERT(test_equal(computed, 2 * IJK)); + ++flatIndex; + } + } + } +} + } // anonymous namespace int UnitTestWorkletMapPointNeighborhood(int argc, char* argv[])