Add CellNeighborhood

This commit is contained in:
Li-Ta Lo 2020-07-15 14:41:32 -06:00
parent 60339bffe6
commit e52b8fa88a
16 changed files with 733 additions and 111 deletions

@ -19,9 +19,10 @@ if(VTKm_ENABLE_EXAMPLES)
add_subdirectory(contour_tree_distributed)
add_subdirectory(cosmotools)
add_subdirectory(demo)
#add_subdirectory(game_of_life)
add_subdirectory(game_of_life)
add_subdirectory(hello_worklet)
add_subdirectory(histogram)
add_subdirectory(ising)
add_subdirectory(lagrangian)
add_subdirectory(mesh_quality)
add_subdirectory(multi_backend)

@ -0,0 +1,21 @@
##============================================================================
## 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.
##============================================================================
cmake_minimum_required(VERSION 3.12...3.15 FATAL_ERROR)
project(IsingModel CXX)
#Find the VTK-m package
find_package(VTKm REQUIRED QUIET)
add_executable(Ising Ising.cxx)
target_link_libraries(Ising PRIVATE vtkm_worklet vtkm_io vtkm_rendering)
vtkm_add_target_information(Ising
DROP_UNUSED_SYMBOLS MODIFY_CUDA_FLAGS
DEVICE_SOURCES Ising.cxx)

121
examples/ising/Ising.cxx Normal file

@ -0,0 +1,121 @@
//
// Created by ollie on 7/8/20.
//
//============================================================================
// 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.
//============================================================================
/// Simulation of ferromagnetism using the Ising Model
/// Reference: Computational Physics 2nd Edition, Nicholas Giordano & Hisao Nakanishi
#include <iomanip>
#include <vtkm/cont/ArrayHandleRandomUniformReal.h>
#include <vtkm/cont/DataSet.h>
#include <vtkm/cont/DataSetBuilderUniform.h>
#include <vtkm/cont/Initialize.h>
#include <vtkm/rendering/CanvasRayTracer.h>
#include <vtkm/rendering/MapperRayTracer.h>
#include <vtkm/rendering/Scene.h>
#include <vtkm/rendering/View2D.h>
struct UpDown
{
VTKM_EXEC_CONT vtkm::Float32 operator()(vtkm::Float32 p) const { return p > 0.5 ? 1.0f : -1.0f; }
};
vtkm::cont::DataSet SpinField(vtkm::Id2 dims)
{
auto result =
vtkm::cont::DataSetBuilderUniform::Create(dims, vtkm::Vec2f{ 0, 0 }, vtkm::Vec2f{ 1, 1 });
vtkm::cont::ArrayHandle<vtkm::Float32> spins;
vtkm::cont::ArrayCopy(
vtkm::cont::make_ArrayHandleTransform(
vtkm::cont::ArrayHandleRandomUniformReal<vtkm::Float32>(result.GetNumberOfCells()), UpDown{}),
spins);
result.AddCellField("spins", spins);
return result;
}
struct UpdateSpins : public vtkm::worklet::WorkletCellNeighborhood
{
using ControlSignature = void(CellSetIn,
FieldInNeighborhood prevspin,
FieldIn prob,
FieldOut spin);
using ExecutionSignature = void(_2, _3, _4);
template <typename NeighIn>
VTKM_EXEC_CONT void operator()(const NeighIn& prevspin,
vtkm::Float32 p,
vtkm::Float32& spin) const
{
// TODO: what is the real value and unit of the change constant J and Boltzmann constant kB?
const vtkm::Float32 J = 1.f;
const vtkm::Float32 kB = 1.f;
// TODO: temperature in Kelvin
const vtkm::Float32 T = 5.f;
const auto mySpin = prevspin.Get(0, 0, 0);
// 1. Calculate the energy of flipping, E_flip
vtkm::Float32 E_flip =
J * mySpin * (prevspin.Get(-1, -1, 0) + prevspin.Get(-1, 0, 0) + prevspin.Get(-1, 1, 0) +
prevspin.Get(0, -1, 0) + prevspin.Get(0, 1, 0) + prevspin.Get(1, -1, 0) +
prevspin.Get(1, 0, 0) + prevspin.Get(1, 1, 0));
if (E_flip <= 0)
{
// 2. If E_flip <= 0, just flip the spin
spin = -1.f * mySpin;
}
else
{
// 3. otherwise, flip the spin if the Boltzmann factor exp(-E_flip/kB*T) is larger than the
// uniform real random number p.
if (p <= vtkm::Exp(-E_flip / (kB * T)))
spin = -1.f * mySpin;
else
spin = mySpin;
}
}
};
int main(int argc, char** argv)
{
auto opts =
vtkm::cont::InitializeOptions::DefaultAnyDevice | vtkm::cont::InitializeOptions::Strict;
vtkm::cont::Initialize(argc, argv, opts);
auto dataSet = SpinField({ 5, 5 });
vtkm::cont::ArrayHandle<vtkm::Float32> spins;
dataSet.GetCellField("spins").GetData().CopyTo(spins);
vtkm::rendering::Scene scene;
vtkm::rendering::Actor actor(dataSet.GetCellSet(),
dataSet.GetCoordinateSystem(),
dataSet.GetCellField("spins"),
vtkm::cont::ColorTable("Cool To Warm"));
scene.AddActor(actor);
vtkm::rendering::CanvasRayTracer canvas(1024, 1024);
vtkm::rendering::MapperRayTracer mapper;
mapper.SetShadingOn(false);
vtkm::rendering::View2D view(scene, mapper, canvas);
view.Paint();
view.SaveAs("spin0.png");
vtkm::cont::Invoker invoker;
for (int i = 1; i < 10; ++i)
{
vtkm::cont::ArrayHandleRandomUniformReal<vtkm::Float32> prob(dataSet.GetNumberOfCells(), { i });
invoker(UpdateSpins{}, dataSet.GetCellSet(), spins, prob, spins);
view.Paint();
view.SaveAs("spin" + std::to_string(i) + ".png");
}
}

@ -10,6 +10,7 @@
#ifndef vtk_m_cont_Invoker_h
#define vtk_m_cont_Invoker_h
#include <vtkm/worklet/DispatcherCellNeighborhood.h>
#include <vtkm/worklet/DispatcherMapField.h>
#include <vtkm/worklet/DispatcherMapTopology.h>
#include <vtkm/worklet/DispatcherPointNeighborhood.h>

@ -115,6 +115,12 @@ public:
return this->Internals.GetPointDimensions();
}
VTKM_EXEC_CONT
vtkm::Vec<vtkm::Id, Dimension> GetCellDimensions() const
{
return this->Internals.GetCellDimensions();
}
VTKM_EXEC_CONT
SchedulingRangeType GetGlobalPointIndexStart() const
{

@ -27,7 +27,7 @@ namespace exec
/// be returned determined by the boundary behavior. A \c BoundaryState object can be used to
/// determine if the neighborhood extends beyond the boundary of the mesh.
///
/// This class is typically constructued using the \c FieldInNeighborhood tag in an
/// This class is typically constructed using the \c FieldInNeighborhood tag in an
/// \c ExecutionSignature. There is little reason to construct this in user code.
///
/// \c FieldNeighborhood is templated on the array portal from which field values are retrieved.

@ -33,6 +33,7 @@ set(headers
ThreadIndicesBasic.h
ThreadIndicesBasic3D.h
ThreadIndicesExtrude.h
ThreadIndicesCellNeighborhood.h
ThreadIndicesPointNeighborhood.h
ThreadIndicesReduceByKey.h
ThreadIndicesTopologyMap.h

@ -0,0 +1,89 @@
//============================================================================
// 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_arg_ThreadIndicesCellNeighborhood_h
#define vtk_m_exec_arg_ThreadIndicesCellNeighborhood_h
#include <vtkm/exec/BoundaryState.h>
#include <vtkm/exec/ConnectivityStructured.h>
#include <vtkm/exec/arg/ThreadIndicesBasic.h>
#include <vtkm/exec/arg/ThreadIndicesTopologyMap.h> //for Deflate and Inflate
#include <vtkm/Math.h>
namespace vtkm
{
namespace exec
{
namespace arg
{
/// \brief Container for thread information in a WorkletCellNeighborhood.
///
///
class ThreadIndicesCellNeighborhood : public vtkm::exec::arg::ThreadIndicesNeighborhood
{
using Superclass = vtkm::exec::arg::ThreadIndicesNeighborhood;
public:
template <vtkm::IdComponent Dimension>
VTKM_EXEC ThreadIndicesCellNeighborhood(
const vtkm::Id3& threadIndex3D,
vtkm::Id threadIndex1D,
const vtkm::exec::ConnectivityStructured<vtkm::TopologyElementTagPoint,
vtkm::TopologyElementTagCell,
Dimension>& connectivity)
: Superclass(
threadIndex1D,
vtkm::exec::BoundaryState{ threadIndex3D, detail::To3D(connectivity.GetCellDimensions()) })
{
}
template <vtkm::IdComponent Dimension>
VTKM_EXEC ThreadIndicesCellNeighborhood(
const vtkm::Id3& threadIndex3D,
vtkm::Id threadIndex1D,
vtkm::Id inputIndex,
vtkm::IdComponent visitIndex,
vtkm::Id outputIndex,
const vtkm::exec::ConnectivityStructured<vtkm::TopologyElementTagPoint,
vtkm::TopologyElementTagCell,
Dimension>& connectivity)
: Superclass(
threadIndex1D,
inputIndex,
visitIndex,
outputIndex,
vtkm::exec::BoundaryState{ threadIndex3D, detail::To3D(connectivity.GetCellDimensions()) })
{
}
template <vtkm::IdComponent Dimension>
VTKM_EXEC ThreadIndicesCellNeighborhood(
vtkm::Id threadIndex,
vtkm::Id inputIndex,
vtkm::IdComponent visitIndex,
vtkm::Id outputIndex,
const vtkm::exec::ConnectivityStructured<vtkm::TopologyElementTagPoint,
vtkm::TopologyElementTagCell,
Dimension>& connectivity)
: Superclass(
threadIndex,
inputIndex,
visitIndex,
outputIndex,
vtkm::exec::BoundaryState{ detail::To3D(connectivity.FlatToLogicalToIndex(inputIndex)),
detail::To3D(connectivity.GetCellDimensions()) })
{
}
};
}
}
} // namespace vtkm::exec::arg
#endif //vtk_m_exec_arg_ThreadIndicesCellNeighborhood_h

@ -0,0 +1,115 @@
//============================================================================
// 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_arg_ThreadIndicesNeighborhood_h
#define vtk_m_exec_arg_ThreadIndicesNeighborhood_h
#include <vtkm/exec/BoundaryState.h>
#include <vtkm/exec/ConnectivityStructured.h>
#include <vtkm/exec/arg/ThreadIndicesBasic.h>
#include <vtkm/exec/arg/ThreadIndicesTopologyMap.h> //for Deflate and Inflate
#include <vtkm/Math.h>
namespace vtkm
{
namespace exec
{
namespace arg
{
namespace detail
{
/// Given a \c Vec of (semi) arbitrary size, inflate it to a vtkm::Id3 by padding with zeros.
///
inline VTKM_EXEC vtkm::Id3 To3D(vtkm::Id3 index)
{
return index;
}
/// Given a \c Vec of (semi) arbitrary size, inflate it to a vtkm::Id3 by padding with zeros.
/// \overload
inline VTKM_EXEC vtkm::Id3 To3D(vtkm::Id2 index)
{
return vtkm::Id3(index[0], index[1], 1);
}
/// Given a \c Vec of (semi) arbitrary size, inflate it to a vtkm::Id3 by padding with zeros.
/// \overload
inline VTKM_EXEC vtkm::Id3 To3D(vtkm::Vec<vtkm::Id, 1> index)
{
return vtkm::Id3(index[0], 1, 1);
}
/// Given a \c Vec of (semi) arbitrary size, inflate it to a vtkm::Id3 by padding with zeros.
/// \overload
inline VTKM_EXEC vtkm::Id3 To3D(vtkm::Id index)
{
return vtkm::Id3(index, 1, 1);
}
}
class ThreadIndicesNeighborhood
{
public:
VTKM_EXEC ThreadIndicesNeighborhood(vtkm::Id threadIndex1D,
const vtkm::exec::BoundaryState& state)
: State(state)
, ThreadIndex(threadIndex1D)
, InputIndex(threadIndex1D)
, OutputIndex(threadIndex1D)
, VisitIndex(0)
{
}
VTKM_EXEC ThreadIndicesNeighborhood(vtkm::Id threadIndex1D,
vtkm::Id inputIndex,
vtkm::IdComponent visitIndex,
vtkm::Id outputIndex,
const vtkm::exec::BoundaryState& state)
: State(state)
, ThreadIndex(threadIndex1D)
, InputIndex(inputIndex)
, OutputIndex(outputIndex)
, VisitIndex(visitIndex)
{
}
VTKM_EXEC
const vtkm::exec::BoundaryState& GetBoundaryState() const { return this->State; }
VTKM_EXEC
vtkm::Id GetThreadIndex() const { return this->ThreadIndex; }
VTKM_EXEC
vtkm::Id GetInputIndex() const { return this->InputIndex; }
VTKM_EXEC
vtkm::Id3 GetInputIndex3D() const { return this->State.IJK; }
VTKM_EXEC
vtkm::Id GetOutputIndex() const { return this->OutputIndex; }
VTKM_EXEC
vtkm::IdComponent GetVisitIndex() const { return this->VisitIndex; }
private:
vtkm::exec::BoundaryState State;
vtkm::Id ThreadIndex;
vtkm::Id InputIndex;
vtkm::Id OutputIndex;
vtkm::IdComponent VisitIndex;
};
}
}
} // namespace vtkm::exec::arg
#endif //vtk_m_exec_arg_ThreadIndicesNeighborhood_h

@ -10,12 +10,7 @@
#ifndef vtk_m_exec_arg_ThreadIndicesPointNeighborhood_h
#define vtk_m_exec_arg_ThreadIndicesPointNeighborhood_h
#include <vtkm/exec/BoundaryState.h>
#include <vtkm/exec/ConnectivityStructured.h>
#include <vtkm/exec/arg/ThreadIndicesBasic.h>
#include <vtkm/exec/arg/ThreadIndicesTopologyMap.h> //for Deflate and Inflate
#include <vtkm/Math.h>
#include <vtkm/exec/arg/ThreadIndicesNeighborhood.h>
namespace vtkm
{
@ -23,43 +18,12 @@ namespace exec
{
namespace arg
{
namespace detail
{
/// Given a \c Vec of (semi) arbitrary size, inflate it to a vtkm::Id3 by padding with zeros.
///
inline VTKM_EXEC vtkm::Id3 To3D(vtkm::Id3 index)
{
return index;
}
/// Given a \c Vec of (semi) arbitrary size, inflate it to a vtkm::Id3 by padding with zeros.
/// \overload
inline VTKM_EXEC vtkm::Id3 To3D(vtkm::Id2 index)
{
return vtkm::Id3(index[0], index[1], 1);
}
/// Given a \c Vec of (semi) arbitrary size, inflate it to a vtkm::Id3 by padding with zeros.
/// \overload
inline VTKM_EXEC vtkm::Id3 To3D(vtkm::Vec<vtkm::Id, 1> index)
{
return vtkm::Id3(index[0], 1, 1);
}
/// Given a \c Vec of (semi) arbitrary size, inflate it to a vtkm::Id3 by padding with zeros.
/// \overload
inline VTKM_EXEC vtkm::Id3 To3D(vtkm::Id index)
{
return vtkm::Id3(index, 1, 1);
}
}
/// \brief Container for thread information in a WorkletPointNeighborhood.
///
///
class ThreadIndicesPointNeighborhood
class ThreadIndicesPointNeighborhood : public vtkm::exec::arg::ThreadIndicesNeighborhood
{
using Superclass = vtkm::exec::arg::ThreadIndicesNeighborhood;
public:
template <vtkm::IdComponent Dimension>
@ -69,11 +33,9 @@ public:
const vtkm::exec::ConnectivityStructured<vtkm::TopologyElementTagPoint,
vtkm::TopologyElementTagCell,
Dimension>& connectivity)
: State(threadIndex3D, detail::To3D(connectivity.GetPointDimensions()))
, ThreadIndex(threadIndex1D)
, InputIndex(threadIndex1D)
, OutputIndex(threadIndex1D)
, VisitIndex(0)
: Superclass(
threadIndex1D,
vtkm::exec::BoundaryState{ threadIndex3D, detail::To3D(connectivity.GetPointDimensions()) })
{
}
@ -87,11 +49,12 @@ public:
const vtkm::exec::ConnectivityStructured<vtkm::TopologyElementTagPoint,
vtkm::TopologyElementTagCell,
Dimension>& connectivity)
: State(threadIndex3D, detail::To3D(connectivity.GetPointDimensions()))
, ThreadIndex(threadIndex1D)
, InputIndex(inputIndex)
, OutputIndex(outputIndex)
, VisitIndex(visitIndex)
: Superclass(
threadIndex1D,
inputIndex,
visitIndex,
outputIndex,
vtkm::exec::BoundaryState{ threadIndex3D, detail::To3D(connectivity.GetPointDimensions()) })
{
}
@ -104,42 +67,17 @@ public:
const vtkm::exec::ConnectivityStructured<vtkm::TopologyElementTagPoint,
vtkm::TopologyElementTagCell,
Dimension>& connectivity)
: State(detail::To3D(connectivity.FlatToLogicalToIndex(inputIndex)),
detail::To3D(connectivity.GetPointDimensions()))
, ThreadIndex(threadIndex)
, InputIndex(inputIndex)
, OutputIndex(outputIndex)
, VisitIndex(visitIndex)
: Superclass(
threadIndex,
inputIndex,
visitIndex,
outputIndex,
vtkm::exec::BoundaryState{ detail::To3D(connectivity.FlatToLogicalToIndex(inputIndex)),
detail::To3D(connectivity.GetPointDimensions()) })
{
}
VTKM_EXEC
const vtkm::exec::BoundaryState& GetBoundaryState() const { return this->State; }
VTKM_EXEC
vtkm::Id GetThreadIndex() const { return this->ThreadIndex; }
VTKM_EXEC
vtkm::Id GetInputIndex() const { return this->InputIndex; }
VTKM_EXEC
vtkm::Id3 GetInputIndex3D() const { return this->State.IJK; }
VTKM_EXEC
vtkm::Id GetOutputIndex() const { return this->OutputIndex; }
VTKM_EXEC
vtkm::IdComponent GetVisitIndex() const { return this->VisitIndex; }
private:
vtkm::exec::BoundaryState State;
vtkm::Id ThreadIndex;
vtkm::Id InputIndex;
vtkm::Id OutputIndex;
vtkm::IdComponent VisitIndex;
};
}
}
} // namespace vtkm::exec::arg
} // arg
} // exec
} // vtkm
#endif //vtk_m_exec_arg_ThreadIndicesPointNeighborhood_h

@ -11,7 +11,7 @@
#define vtk_m_interop_internal_TransferToOpenGL_h
#include <vtkm/cont/ArrayHandle.h>
#include <vtkm/cont/StorageBasic.h>
#include <vtkm/cont/Storage.h>
#include <vtkm/cont/DeviceAdapterAlgorithm.h>
#include <vtkm/cont/serial/DeviceAdapterSerial.h>

@ -0,0 +1,41 @@
//============================================================================
// 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_worklet_BoundaryTypes_h
#define vtk_m_worklet_BoundaryTypes_h
namespace vtkm
{
namespace worklet
{
/// \brief Clamps boundary values to the nearest valid i,j,k value
///
/// BoundaryClamp always returns the nearest valid i,j,k value when at an
/// image boundary. This is a commonly used when solving differential equations.
///
/// For example, when used with WorkletCellNeighborhood3x3x3 when centered
/// on the point 1:
/// \code
/// * * *
/// * 1 2 (where * denotes points that lie outside of the image boundary)
/// * 3 5
/// \endcode
/// returns the following neighborhood of values:
/// \code
/// 1 1 2
/// 1 1 2
/// 3 3 5
/// \endcode
struct BoundaryClamp
{
};
}
}
#endif //vtk_m_worklet_BoundaryTypes_h

@ -10,6 +10,7 @@
set(headers
AverageByKey.h
BoundaryTypes.h
CellAverage.h
CellDeepCopy.h
CellMeasure.h
@ -21,6 +22,7 @@ set(headers
CrossProduct.h
DispatcherMapField.h
DispatcherMapTopology.h
DispatcherCellNeighborhood.h
DispatcherPointNeighborhood.h
DispatcherReduceByKey.h
DotProduct.h
@ -84,6 +86,7 @@ set(headers
WaveletCompressor.h
WorkletMapField.h
WorkletMapTopology.h
WorkletCellNeighborhood.h
WorkletPointNeighborhood.h
WorkletReduceByKey.h
ZFPCompressor.h

@ -0,0 +1,71 @@
//============================================================================
// 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_worklet_DispatcherCellNeighborhood_h
#define vtk_m_worklet_DispatcherCellNeighborhood_h
#include <vtkm/cont/DeviceAdapter.h>
#include <vtkm/worklet/WorkletCellNeighborhood.h>
#include <vtkm/worklet/internal/DispatcherBase.h>
namespace vtkm
{
namespace worklet
{
/// \brief Dispatcher for worklets that inherit from \c WorkletCellNeighborhood.
///
template <typename WorkletType>
class DispatcherCellNeighborhood
: public vtkm::worklet::internal::DispatcherBase<DispatcherCellNeighborhood<WorkletType>,
WorkletType,
vtkm::worklet::WorkletCellNeighborhoodBase>
{
using Superclass =
vtkm::worklet::internal::DispatcherBase<DispatcherCellNeighborhood<WorkletType>,
WorkletType,
vtkm::worklet::WorkletCellNeighborhoodBase>;
using ScatterType = typename Superclass::ScatterType;
public:
template <typename... T>
VTKM_CONT DispatcherCellNeighborhood(T&&... args)
: Superclass(std::forward<T>(args)...)
{
}
template <typename Invocation>
void DoInvoke(Invocation& invocation) const
{
using namespace vtkm::worklet::internal;
// This is the type for the input domain
using InputDomainType = typename Invocation::InputDomainType;
// If you get a compile error on this line, then you have tried to use
// something that is not a vtkm::cont::CellSet as the input domain to a
// topology operation (that operates on a cell set connection domain).
VTKM_IS_CELL_SET(InputDomainType);
// We can pull the input domain parameter (the data specifying the input
// domain) from the invocation object.
const InputDomainType& inputDomain = invocation.GetInputDomain();
auto inputRange = SchedulingRange(inputDomain, vtkm::TopologyElementTagCell{});
// This is pretty straightforward dispatch. Once we know the number
// of invocations, the superclass can take care of the rest.
this->BasicInvoke(invocation, inputRange);
}
};
}
} // namespace vtkm::worklet
#endif //vtk_m_worklet_DispatcherCellNeighborhood_h

@ -0,0 +1,237 @@
//============================================================================
// 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_worklet_WorkletCellNeighborhood_h
#define vtk_m_worklet_WorkletCellNeighborhood_h
/// \brief Worklet for volume algorithms that require a neighborhood
///
/// WorkletCellNeighborhood executes on every point inside a volume providing
/// access to the 3D neighborhood values. The neighborhood is always cubic in
/// nature and is fixed at compile time.
#include <vtkm/worklet/internal/WorkletBase.h>
#include <vtkm/TopologyElementTag.h>
#include <vtkm/cont/arg/ControlSignatureTagBase.h>
#include <vtkm/cont/arg/TransportTagArrayIn.h>
#include <vtkm/cont/arg/TransportTagArrayInOut.h>
#include <vtkm/cont/arg/TransportTagArrayOut.h>
#include <vtkm/cont/arg/TransportTagCellSetIn.h>
#include <vtkm/cont/arg/TypeCheckTagArrayIn.h>
#include <vtkm/cont/arg/TypeCheckTagArrayInOut.h>
#include <vtkm/cont/arg/TypeCheckTagArrayOut.h>
#include <vtkm/cont/arg/TypeCheckTagCellSetStructured.h>
#include <vtkm/exec/arg/Boundary.h>
#include <vtkm/exec/arg/FetchTagArrayDirectIn.h>
#include <vtkm/exec/arg/FetchTagArrayDirectInOut.h>
#include <vtkm/exec/arg/FetchTagArrayDirectOut.h>
#include <vtkm/exec/arg/FetchTagArrayNeighborhoodIn.h>
#include <vtkm/exec/arg/FetchTagCellSetIn.h>
#include <vtkm/exec/arg/ThreadIndicesCellNeighborhood.h>
#include <vtkm/worklet/BoundaryTypes.h>
#include <vtkm/worklet/ScatterIdentity.h>
namespace vtkm
{
namespace worklet
{
template <typename WorkletType>
class DispatcherCellNeighborhood;
class WorkletCellNeighborhoodBase : public vtkm::worklet::internal::WorkletBase
{
public:
template <typename Worklet>
using Dispatcher = vtkm::worklet::DispatcherCellNeighborhood<Worklet>;
/// \brief The \c ExecutionSignature tag to query if the current iteration is inside the boundary.
///
/// A \c WorkletCellNeighborhood operates by iterating over all points using a defined
/// neighborhood. This \c ExecutionSignature tag provides a \c BoundaryState object that allows
/// you to query whether the neighborhood of the current iteration is completely inside the
/// bounds of the mesh or if it extends beyond the mesh. This is important as when you are on a
/// boundary the neighboordhood will contain empty values for a certain subset of values, and in
/// this case the values returned will depend on the boundary behavior.
///
struct Boundary : vtkm::exec::arg::Boundary
{
};
/// All worklets must define their scatter operation.
using ScatterType = vtkm::worklet::ScatterIdentity;
/// All neighborhood worklets must define their boundary type operation.
/// The boundary type determines how loading on boundaries will work.
using BoundaryType = vtkm::worklet::BoundaryClamp;
/// In addition to defining the boundary type, the worklet must produce the
/// boundary condition. The default BoundaryClamp has no state, so just return an
/// instance.
/// Note: Currently only BoundaryClamp is implemented
VTKM_CONT
BoundaryType GetBoundaryCondition() const { return BoundaryType(); }
/// \brief A control signature tag for input point fields.
///
/// This tag takes a template argument that is a type list tag that limits
/// the possible value types in the array.
///
struct FieldIn : vtkm::cont::arg::ControlSignatureTagBase
{
using TypeCheckTag = vtkm::cont::arg::TypeCheckTagArrayIn;
using TransportTag = vtkm::cont::arg::TransportTagArrayIn;
using FetchTag = vtkm::exec::arg::FetchTagArrayDirectIn;
};
/// \brief A control signature tag for output point fields.
///
/// This tag takes a template argument that is a type list tag that limits
/// the possible value types in the array.
///
struct FieldOut : vtkm::cont::arg::ControlSignatureTagBase
{
using TypeCheckTag = vtkm::cont::arg::TypeCheckTagArrayOut;
using TransportTag = vtkm::cont::arg::TransportTagArrayOut;
using FetchTag = vtkm::exec::arg::FetchTagArrayDirectOut;
};
/// \brief A control signature tag for input-output (in-place) point fields.
///
/// This tag takes a template argument that is a type list tag that limits
/// the possible value types in the array.
///
struct FieldInOut : vtkm::cont::arg::ControlSignatureTagBase
{
using TypeCheckTag = vtkm::cont::arg::TypeCheckTagArrayInOut;
using TransportTag = vtkm::cont::arg::TransportTagArrayInOut;
using FetchTag = vtkm::exec::arg::FetchTagArrayDirectInOut;
};
/// \brief A control signature tag for input connectivity.
///
struct CellSetIn : vtkm::cont::arg::ControlSignatureTagBase
{
using TypeCheckTag = vtkm::cont::arg::TypeCheckTagCellSetStructured;
using TransportTag = vtkm::cont::arg::TransportTagCellSetIn<vtkm::TopologyElementTagPoint,
vtkm::TopologyElementTagCell>;
using FetchTag = vtkm::exec::arg::FetchTagCellSetIn;
};
};
class WorkletCellNeighborhood : public WorkletCellNeighborhoodBase
{
public:
/// \brief A control signature tag for neighborhood input values.
///
/// A \c WorkletCellNeighborhood operates allowing access to a adjacent point
/// values in a NxNxN patch called a neighborhood.
/// No matter the size of the neighborhood it is symmetric across its center
/// in each axis, and the current point value will be at the center
/// For example a 3x3x3 neighborhood would
///
/// This tag specifies an \c ArrayHandle object that holds the values. It is
/// an input array with entries for each point.
///
struct FieldInNeighborhood : vtkm::cont::arg::ControlSignatureTagBase
{
using TypeCheckTag = vtkm::cont::arg::TypeCheckTagArrayIn;
using TransportTag = vtkm::cont::arg::TransportTagArrayIn;
using FetchTag = vtkm::exec::arg::FetchTagArrayNeighborhoodIn;
};
/// Point neighborhood worklets use the related thread indices class.
///
VTKM_SUPPRESS_EXEC_WARNINGS
template <typename OutToInArrayType,
typename VisitArrayType,
typename ThreadToOutArrayType,
vtkm::IdComponent Dimension>
VTKM_EXEC vtkm::exec::arg::ThreadIndicesCellNeighborhood GetThreadIndices(
vtkm::Id threadIndex,
const OutToInArrayType& outToIn,
const VisitArrayType& visit,
const ThreadToOutArrayType& threadToOut,
const vtkm::exec::ConnectivityStructured<vtkm::TopologyElementTagPoint,
vtkm::TopologyElementTagCell,
Dimension>& inputDomain //this should be explicit
) const
{
const vtkm::Id outIndex = threadToOut.Get(threadIndex);
return vtkm::exec::arg::ThreadIndicesCellNeighborhood(
threadIndex, outToIn.Get(outIndex), visit.Get(outIndex), outIndex, inputDomain);
}
/// In the remaining methods and `constexpr` we determine at compilation time
/// which method definition will be actually used for GetThreadIndices.
///
/// We want to avoid further function calls when we use WorkletMapTopology in which
/// ScatterType is set as ScatterIdentity and MaskType as MaskNone.
/// Otherwise, we call the default method defined at the bottom of this class.
private:
static constexpr bool IsScatterIdentity =
std::is_same<ScatterType, vtkm::worklet::ScatterIdentity>::value;
static constexpr bool IsMaskNone = std::is_same<MaskType, vtkm::worklet::MaskNone>::value;
public:
template <bool Cond, typename ReturnType>
using EnableFnWhen = typename std::enable_if<Cond, ReturnType>::type;
VTKM_SUPPRESS_EXEC_WARNINGS
template <typename OutToInArrayType,
typename VisitArrayType,
typename ThreadToOutArrayType,
typename InputDomainType,
bool S = IsScatterIdentity,
bool M = IsMaskNone>
VTKM_EXEC EnableFnWhen<S && M, vtkm::exec::arg::ThreadIndicesCellNeighborhood> GetThreadIndices(
vtkm::Id threadIndex1D,
const vtkm::Id3& threadIndex3D,
const OutToInArrayType& vtkmNotUsed(outToIn),
const VisitArrayType& vtkmNotUsed(visit),
const ThreadToOutArrayType& vtkmNotUsed(threadToOut),
const InputDomainType& connectivity) const
{
return vtkm::exec::arg::ThreadIndicesCellNeighborhood(
threadIndex3D, threadIndex1D, connectivity);
}
VTKM_SUPPRESS_EXEC_WARNINGS
template <typename OutToInArrayType,
typename VisitArrayType,
typename ThreadToOutArrayType,
typename InputDomainType,
bool S = IsScatterIdentity,
bool M = IsMaskNone>
VTKM_EXEC EnableFnWhen<!(S && M), vtkm::exec::arg::ThreadIndicesCellNeighborhood>
GetThreadIndices(vtkm::Id threadIndex1D,
const vtkm::Id3& threadIndex3D,
const OutToInArrayType& outToIn,
const VisitArrayType& visit,
const ThreadToOutArrayType& threadToOut,
const InputDomainType& connectivity) const
{
const vtkm::Id outIndex = threadToOut.Get(threadIndex1D);
return vtkm::exec::arg::ThreadIndicesCellNeighborhood(threadIndex3D,
threadIndex1D,
outToIn.Get(outIndex),
visit.Get(outIndex),
outIndex,
connectivity);
}
};
}
}
#endif

@ -38,9 +38,9 @@
#include <vtkm/exec/arg/FetchTagCellSetIn.h>
#include <vtkm/exec/arg/ThreadIndicesPointNeighborhood.h>
#include <vtkm/worklet/BoundaryTypes.h>
#include <vtkm/worklet/ScatterIdentity.h>
namespace vtkm
{
namespace worklet
@ -49,29 +49,6 @@ namespace worklet
template <typename WorkletType>
class DispatcherPointNeighborhood;
/// \brief Clamps boundary values to the nearest valid i,j,k value
///
/// BoundaryClamp always returns the nearest valid i,j,k value when at an
/// image boundary. This is a commonly used when solving differential equations.
///
/// For example, when used with WorkletPointNeighborhood3x3x3 when centered
/// on the point 1:
/// \code
/// * * *
/// * 1 2 (where * denotes points that lie outside of the image boundary)
/// * 3 5
/// \endcode
/// returns the following neighborhood of values:
/// \code
/// 1 1 2
/// 1 1 2
/// 3 3 5
/// \endcode
struct BoundaryClamp
{
};
class WorkletPointNeighborhoodBase : public vtkm::worklet::internal::WorkletBase
{
public: