Merge topic 'particle_adv_virtuals_removal'
b316cf58f Fix compiler warnings 5f7287bc1 Removing virtuals from particle advection 2897b7a11 Fix missing include 4c781374c Removing virtuals v1 Acked-by: Kitware Robot <kwrobot@kitware.com> Acked-by: Kenneth Moreland <kmorel@acm.org> Merge-request: !2445
This commit is contained in:
commit
43400995dd
@ -25,9 +25,9 @@
|
|||||||
#include <vtkm/worklet/WorkletMapField.h>
|
#include <vtkm/worklet/WorkletMapField.h>
|
||||||
#include <vtkm/worklet/particleadvection/Field.h>
|
#include <vtkm/worklet/particleadvection/Field.h>
|
||||||
#include <vtkm/worklet/particleadvection/GridEvaluators.h>
|
#include <vtkm/worklet/particleadvection/GridEvaluators.h>
|
||||||
#include <vtkm/worklet/particleadvection/IntegratorBase.h>
|
|
||||||
#include <vtkm/worklet/particleadvection/Particles.h>
|
#include <vtkm/worklet/particleadvection/Particles.h>
|
||||||
#include <vtkm/worklet/particleadvection/RK4Integrator.h>
|
#include <vtkm/worklet/particleadvection/RK4Integrator.h>
|
||||||
|
#include <vtkm/worklet/particleadvection/Stepper.h>
|
||||||
|
|
||||||
#include <cstring>
|
#include <cstring>
|
||||||
#include <sstream>
|
#include <sstream>
|
||||||
@ -279,12 +279,14 @@ inline VTKM_CONT vtkm::cont::DataSet Lagrangian::DoExecute(
|
|||||||
using FieldType = vtkm::worklet::particleadvection::VelocityField<FieldHandle>;
|
using FieldType = vtkm::worklet::particleadvection::VelocityField<FieldHandle>;
|
||||||
using GridEvalType = vtkm::worklet::particleadvection::GridEvaluator<FieldType>;
|
using GridEvalType = vtkm::worklet::particleadvection::GridEvaluator<FieldType>;
|
||||||
using RK4Type = vtkm::worklet::particleadvection::RK4Integrator<GridEvalType>;
|
using RK4Type = vtkm::worklet::particleadvection::RK4Integrator<GridEvalType>;
|
||||||
|
using Stepper = vtkm::worklet::particleadvection::Stepper<RK4Type, GridEvalType>;
|
||||||
|
|
||||||
vtkm::worklet::ParticleAdvection particleadvection;
|
vtkm::worklet::ParticleAdvection particleadvection;
|
||||||
vtkm::worklet::ParticleAdvectionResult<vtkm::Particle> res;
|
vtkm::worklet::ParticleAdvectionResult<vtkm::Particle> res;
|
||||||
|
|
||||||
FieldType velocities(field);
|
FieldType velocities(field);
|
||||||
GridEvalType gridEval(coords, cells, velocities);
|
GridEvalType gridEval(coords, cells, velocities);
|
||||||
RK4Type rk4(gridEval, static_cast<vtkm::Float32>(this->stepSize));
|
Stepper rk4(gridEval, static_cast<vtkm::Float32>(this->stepSize));
|
||||||
|
|
||||||
res = particleadvection.Run(rk4, basisParticleArray, 1); // Taking a single step
|
res = particleadvection.Run(rk4, basisParticleArray, 1); // Taking a single step
|
||||||
auto particles = res.Particles;
|
auto particles = res.Particles;
|
||||||
|
@ -17,7 +17,7 @@
|
|||||||
|
|
||||||
#include <vtkm/worklet/ParticleAdvection.h>
|
#include <vtkm/worklet/ParticleAdvection.h>
|
||||||
#include <vtkm/worklet/particleadvection/GridEvaluators.h>
|
#include <vtkm/worklet/particleadvection/GridEvaluators.h>
|
||||||
#include <vtkm/worklet/particleadvection/IntegratorBase.h>
|
#include <vtkm/worklet/particleadvection/Stepper.h>
|
||||||
|
|
||||||
namespace vtkm
|
namespace vtkm
|
||||||
{
|
{
|
||||||
|
@ -16,9 +16,9 @@
|
|||||||
#include <vtkm/cont/Invoker.h>
|
#include <vtkm/cont/Invoker.h>
|
||||||
#include <vtkm/worklet/particleadvection/Field.h>
|
#include <vtkm/worklet/particleadvection/Field.h>
|
||||||
#include <vtkm/worklet/particleadvection/GridEvaluators.h>
|
#include <vtkm/worklet/particleadvection/GridEvaluators.h>
|
||||||
#include <vtkm/worklet/particleadvection/IntegratorBase.h>
|
|
||||||
#include <vtkm/worklet/particleadvection/Particles.h>
|
#include <vtkm/worklet/particleadvection/Particles.h>
|
||||||
#include <vtkm/worklet/particleadvection/RK4Integrator.h>
|
#include <vtkm/worklet/particleadvection/RK4Integrator.h>
|
||||||
|
#include <vtkm/worklet/particleadvection/Stepper.h>
|
||||||
|
|
||||||
#include <vtkm/worklet/LagrangianStructures.h>
|
#include <vtkm/worklet/LagrangianStructures.h>
|
||||||
|
|
||||||
@ -86,7 +86,8 @@ inline VTKM_CONT vtkm::cont::DataSet LagrangianStructures::DoExecute(
|
|||||||
using FieldHandle = vtkm::cont::ArrayHandle<vtkm::Vec<T, 3>, StorageType>;
|
using FieldHandle = vtkm::cont::ArrayHandle<vtkm::Vec<T, 3>, StorageType>;
|
||||||
using FieldType = vtkm::worklet::particleadvection::VelocityField<FieldHandle>;
|
using FieldType = vtkm::worklet::particleadvection::VelocityField<FieldHandle>;
|
||||||
using GridEvaluator = vtkm::worklet::particleadvection::GridEvaluator<FieldType>;
|
using GridEvaluator = vtkm::worklet::particleadvection::GridEvaluator<FieldType>;
|
||||||
using Integrator = vtkm::worklet::particleadvection::RK4Integrator<GridEvaluator>;
|
using IntegratorType = vtkm::worklet::particleadvection::RK4Integrator<GridEvaluator>;
|
||||||
|
using Stepper = vtkm::worklet::particleadvection::Stepper<IntegratorType, GridEvaluator>;
|
||||||
|
|
||||||
vtkm::FloatDefault stepSize = this->GetStepSize();
|
vtkm::FloatDefault stepSize = this->GetStepSize();
|
||||||
vtkm::Id numberOfSteps = this->GetNumberOfSteps();
|
vtkm::Id numberOfSteps = this->GetNumberOfSteps();
|
||||||
@ -138,7 +139,7 @@ inline VTKM_CONT vtkm::cont::DataSet LagrangianStructures::DoExecute(
|
|||||||
|
|
||||||
FieldType velocities(field);
|
FieldType velocities(field);
|
||||||
GridEvaluator evaluator(input.GetCoordinateSystem(), input.GetCellSet(), velocities);
|
GridEvaluator evaluator(input.GetCoordinateSystem(), input.GetCellSet(), velocities);
|
||||||
Integrator integrator(evaluator, stepSize);
|
Stepper integrator(evaluator, stepSize);
|
||||||
vtkm::worklet::ParticleAdvection particles;
|
vtkm::worklet::ParticleAdvection particles;
|
||||||
vtkm::worklet::ParticleAdvectionResult<vtkm::Particle> advectionResult;
|
vtkm::worklet::ParticleAdvectionResult<vtkm::Particle> advectionResult;
|
||||||
vtkm::cont::ArrayHandle<vtkm::Particle> advectionPoints;
|
vtkm::cont::ArrayHandle<vtkm::Particle> advectionPoints;
|
||||||
|
@ -14,7 +14,7 @@
|
|||||||
#include <vtkm/filter/FilterDataSetWithField.h>
|
#include <vtkm/filter/FilterDataSetWithField.h>
|
||||||
#include <vtkm/worklet/ParticleAdvection.h>
|
#include <vtkm/worklet/ParticleAdvection.h>
|
||||||
#include <vtkm/worklet/particleadvection/GridEvaluators.h>
|
#include <vtkm/worklet/particleadvection/GridEvaluators.h>
|
||||||
#include <vtkm/worklet/particleadvection/IntegratorBase.h>
|
#include <vtkm/worklet/particleadvection/Stepper.h>
|
||||||
|
|
||||||
namespace vtkm
|
namespace vtkm
|
||||||
{
|
{
|
||||||
|
@ -15,7 +15,7 @@
|
|||||||
#include <vtkm/worklet/ParticleAdvection.h>
|
#include <vtkm/worklet/ParticleAdvection.h>
|
||||||
#include <vtkm/worklet/StreamSurface.h>
|
#include <vtkm/worklet/StreamSurface.h>
|
||||||
#include <vtkm/worklet/particleadvection/GridEvaluators.h>
|
#include <vtkm/worklet/particleadvection/GridEvaluators.h>
|
||||||
#include <vtkm/worklet/particleadvection/IntegratorBase.h>
|
#include <vtkm/worklet/particleadvection/Stepper.h>
|
||||||
|
|
||||||
namespace vtkm
|
namespace vtkm
|
||||||
{
|
{
|
||||||
|
@ -21,6 +21,7 @@
|
|||||||
#include <vtkm/worklet/particleadvection/GridEvaluators.h>
|
#include <vtkm/worklet/particleadvection/GridEvaluators.h>
|
||||||
#include <vtkm/worklet/particleadvection/Particles.h>
|
#include <vtkm/worklet/particleadvection/Particles.h>
|
||||||
#include <vtkm/worklet/particleadvection/RK4Integrator.h>
|
#include <vtkm/worklet/particleadvection/RK4Integrator.h>
|
||||||
|
#include <vtkm/worklet/particleadvection/Stepper.h>
|
||||||
|
|
||||||
namespace vtkm
|
namespace vtkm
|
||||||
{
|
{
|
||||||
@ -57,11 +58,12 @@ inline VTKM_CONT vtkm::cont::DataSet StreamSurface::DoExecute(
|
|||||||
using FieldType = vtkm::worklet::particleadvection::VelocityField<FieldHandle>;
|
using FieldType = vtkm::worklet::particleadvection::VelocityField<FieldHandle>;
|
||||||
using GridEvalType = vtkm::worklet::particleadvection::GridEvaluator<FieldType>;
|
using GridEvalType = vtkm::worklet::particleadvection::GridEvaluator<FieldType>;
|
||||||
using RK4Type = vtkm::worklet::particleadvection::RK4Integrator<GridEvalType>;
|
using RK4Type = vtkm::worklet::particleadvection::RK4Integrator<GridEvalType>;
|
||||||
|
using Stepper = vtkm::worklet::particleadvection::Stepper<RK4Type, GridEvalType>;
|
||||||
|
|
||||||
//compute streamlines
|
//compute streamlines
|
||||||
FieldType velocities(field);
|
FieldType velocities(field);
|
||||||
GridEvalType eval(coords, cells, velocities);
|
GridEvalType eval(coords, cells, velocities);
|
||||||
RK4Type rk4(eval, this->StepSize);
|
Stepper rk4(eval, this->StepSize);
|
||||||
|
|
||||||
vtkm::worklet::Streamline streamline;
|
vtkm::worklet::Streamline streamline;
|
||||||
|
|
||||||
|
@ -15,6 +15,7 @@
|
|||||||
#include <vtkm/cont/DataSet.h>
|
#include <vtkm/cont/DataSet.h>
|
||||||
#include <vtkm/worklet/ParticleAdvection.h>
|
#include <vtkm/worklet/ParticleAdvection.h>
|
||||||
#include <vtkm/worklet/particleadvection/RK4Integrator.h>
|
#include <vtkm/worklet/particleadvection/RK4Integrator.h>
|
||||||
|
#include <vtkm/worklet/particleadvection/Stepper.h>
|
||||||
#include <vtkm/worklet/particleadvection/TemporalGridEvaluators.h>
|
#include <vtkm/worklet/particleadvection/TemporalGridEvaluators.h>
|
||||||
|
|
||||||
#include <memory>
|
#include <memory>
|
||||||
@ -51,17 +52,18 @@ public:
|
|||||||
{
|
{
|
||||||
auto copyFlag = (this->CopySeedArray ? vtkm::CopyFlag::On : vtkm::CopyFlag::Off);
|
auto copyFlag = (this->CopySeedArray ? vtkm::CopyFlag::On : vtkm::CopyFlag::Off);
|
||||||
auto seedArray = vtkm::cont::make_ArrayHandle(v, copyFlag);
|
auto seedArray = vtkm::cont::make_ArrayHandle(v, copyFlag);
|
||||||
RK4Type rk4(*this->Eval, stepSize);
|
Stepper rk4(*this->Eval, stepSize);
|
||||||
this->DoAdvect(seedArray, rk4, maxSteps, result);
|
this->DoAdvect(seedArray, rk4, maxSteps, result);
|
||||||
}
|
}
|
||||||
|
|
||||||
protected:
|
protected:
|
||||||
using RK4Type = vtkm::worklet::particleadvection::RK4Integrator<GridEvalType>;
|
using RK4Type = vtkm::worklet::particleadvection::RK4Integrator<GridEvalType>;
|
||||||
|
using Stepper = vtkm::worklet::particleadvection::Stepper<RK4Type, GridEvalType>;
|
||||||
using FieldHandleType = vtkm::cont::ArrayHandle<vtkm::Vec3f>;
|
using FieldHandleType = vtkm::cont::ArrayHandle<vtkm::Vec3f>;
|
||||||
|
|
||||||
template <typename ResultType>
|
template <typename ResultType>
|
||||||
inline void DoAdvect(vtkm::cont::ArrayHandle<vtkm::Particle>& seeds,
|
inline void DoAdvect(vtkm::cont::ArrayHandle<vtkm::Particle>& seeds,
|
||||||
const RK4Type& rk4,
|
const Stepper& rk4,
|
||||||
vtkm::Id maxSteps,
|
vtkm::Id maxSteps,
|
||||||
ResultType& result) const;
|
ResultType& result) const;
|
||||||
|
|
||||||
|
@ -23,7 +23,10 @@ using GridEvalType = vtkm::worklet::particleadvection::GridEvaluator<
|
|||||||
using TemporalGridEvalType = vtkm::worklet::particleadvection::TemporalGridEvaluator<
|
using TemporalGridEvalType = vtkm::worklet::particleadvection::TemporalGridEvaluator<
|
||||||
vtkm::worklet::particleadvection::VelocityField<vtkm::cont::ArrayHandle<vtkm::Vec3f>>>;
|
vtkm::worklet::particleadvection::VelocityField<vtkm::cont::ArrayHandle<vtkm::Vec3f>>>;
|
||||||
using RK4Type = vtkm::worklet::particleadvection::RK4Integrator<GridEvalType>;
|
using RK4Type = vtkm::worklet::particleadvection::RK4Integrator<GridEvalType>;
|
||||||
|
using Stepper = vtkm::worklet::particleadvection::Stepper<RK4Type, GridEvalType>;
|
||||||
using TemporalRK4Type = vtkm::worklet::particleadvection::RK4Integrator<TemporalGridEvalType>;
|
using TemporalRK4Type = vtkm::worklet::particleadvection::RK4Integrator<TemporalGridEvalType>;
|
||||||
|
using TemporalStepper =
|
||||||
|
vtkm::worklet::particleadvection::Stepper<TemporalRK4Type, TemporalGridEvalType>;
|
||||||
|
|
||||||
//-----
|
//-----
|
||||||
// Specialization for ParticleAdvection worklet
|
// Specialization for ParticleAdvection worklet
|
||||||
@ -31,7 +34,7 @@ template <>
|
|||||||
template <>
|
template <>
|
||||||
inline void DataSetIntegratorBase<GridEvalType>::DoAdvect(
|
inline void DataSetIntegratorBase<GridEvalType>::DoAdvect(
|
||||||
vtkm::cont::ArrayHandle<vtkm::Particle>& seeds,
|
vtkm::cont::ArrayHandle<vtkm::Particle>& seeds,
|
||||||
const RK4Type& rk4,
|
const Stepper& rk4,
|
||||||
vtkm::Id maxSteps,
|
vtkm::Id maxSteps,
|
||||||
vtkm::worklet::ParticleAdvectionResult<vtkm::Particle>& result) const
|
vtkm::worklet::ParticleAdvectionResult<vtkm::Particle>& result) const
|
||||||
{
|
{
|
||||||
@ -45,7 +48,7 @@ template <>
|
|||||||
template <>
|
template <>
|
||||||
inline void DataSetIntegratorBase<GridEvalType>::DoAdvect(
|
inline void DataSetIntegratorBase<GridEvalType>::DoAdvect(
|
||||||
vtkm::cont::ArrayHandle<vtkm::Particle>& seeds,
|
vtkm::cont::ArrayHandle<vtkm::Particle>& seeds,
|
||||||
const RK4Type& rk4,
|
const Stepper& rk4,
|
||||||
vtkm::Id maxSteps,
|
vtkm::Id maxSteps,
|
||||||
vtkm::worklet::StreamlineResult<vtkm::Particle>& result) const
|
vtkm::worklet::StreamlineResult<vtkm::Particle>& result) const
|
||||||
{
|
{
|
||||||
@ -59,7 +62,7 @@ template <>
|
|||||||
template <>
|
template <>
|
||||||
inline void DataSetIntegratorBase<TemporalGridEvalType>::DoAdvect(
|
inline void DataSetIntegratorBase<TemporalGridEvalType>::DoAdvect(
|
||||||
vtkm::cont::ArrayHandle<vtkm::Particle>& seeds,
|
vtkm::cont::ArrayHandle<vtkm::Particle>& seeds,
|
||||||
const TemporalRK4Type& rk4,
|
const TemporalStepper& rk4,
|
||||||
vtkm::Id maxSteps,
|
vtkm::Id maxSteps,
|
||||||
vtkm::worklet::StreamlineResult<vtkm::Particle>& result) const
|
vtkm::worklet::StreamlineResult<vtkm::Particle>& result) const
|
||||||
{
|
{
|
||||||
|
@ -14,7 +14,7 @@ set(headers
|
|||||||
Field.h
|
Field.h
|
||||||
GridEvaluators.h
|
GridEvaluators.h
|
||||||
GridEvaluatorStatus.h
|
GridEvaluatorStatus.h
|
||||||
IntegratorBase.h
|
Stepper.h
|
||||||
IntegratorStatus.h
|
IntegratorStatus.h
|
||||||
Particles.h
|
Particles.h
|
||||||
ParticleAdvectionWorklets.h
|
ParticleAdvectionWorklets.h
|
||||||
|
@ -31,158 +31,144 @@ namespace vtkm
|
|||||||
namespace exec
|
namespace exec
|
||||||
{
|
{
|
||||||
|
|
||||||
class CellInterpolationHelper : public vtkm::VirtualObjectBase
|
class CellInterpolationHelper
|
||||||
{
|
{
|
||||||
public:
|
|
||||||
VTKM_EXEC_CONT virtual ~CellInterpolationHelper() noexcept override
|
|
||||||
{
|
|
||||||
// This must not be defaulted, since defaulted virtual destructors are
|
|
||||||
// troublesome with CUDA __host__ __device__ markup.
|
|
||||||
}
|
|
||||||
|
|
||||||
VTKM_EXEC
|
|
||||||
virtual void GetCellInfo(const vtkm::Id& cellId,
|
|
||||||
vtkm::UInt8& cellShape,
|
|
||||||
vtkm::IdComponent& numVerts,
|
|
||||||
vtkm::VecVariable<vtkm::Id, 8>& indices) const = 0;
|
|
||||||
};
|
|
||||||
|
|
||||||
class StructuredCellInterpolationHelper : public vtkm::exec::CellInterpolationHelper
|
|
||||||
{
|
|
||||||
public:
|
|
||||||
StructuredCellInterpolationHelper() = default;
|
|
||||||
|
|
||||||
VTKM_CONT
|
|
||||||
StructuredCellInterpolationHelper(vtkm::Id3 cellDims, vtkm::Id3 pointDims, bool is3D)
|
|
||||||
: CellDims(cellDims)
|
|
||||||
, PointDims(pointDims)
|
|
||||||
, Is3D(is3D)
|
|
||||||
{
|
|
||||||
}
|
|
||||||
|
|
||||||
VTKM_EXEC
|
|
||||||
void GetCellInfo(const vtkm::Id& cellId,
|
|
||||||
vtkm::UInt8& cellShape,
|
|
||||||
vtkm::IdComponent& numVerts,
|
|
||||||
vtkm::VecVariable<vtkm::Id, 8>& indices) const override
|
|
||||||
{
|
|
||||||
vtkm::Id3 logicalCellId;
|
|
||||||
logicalCellId[0] = cellId % this->CellDims[0];
|
|
||||||
logicalCellId[1] = (cellId / this->CellDims[0]) % this->CellDims[1];
|
|
||||||
if (this->Is3D)
|
|
||||||
{
|
|
||||||
logicalCellId[2] = cellId / (this->CellDims[0] * this->CellDims[1]);
|
|
||||||
indices.Append((logicalCellId[2] * this->PointDims[1] + logicalCellId[1]) *
|
|
||||||
this->PointDims[0] +
|
|
||||||
logicalCellId[0]);
|
|
||||||
indices.Append(indices[0] + 1);
|
|
||||||
indices.Append(indices[1] + this->PointDims[0]);
|
|
||||||
indices.Append(indices[2] - 1);
|
|
||||||
indices.Append(indices[0] + this->PointDims[0] * this->PointDims[1]);
|
|
||||||
indices.Append(indices[4] + 1);
|
|
||||||
indices.Append(indices[5] + this->PointDims[0]);
|
|
||||||
indices.Append(indices[6] - 1);
|
|
||||||
cellShape = static_cast<vtkm::UInt8>(vtkm::CELL_SHAPE_HEXAHEDRON);
|
|
||||||
numVerts = 8;
|
|
||||||
}
|
|
||||||
else
|
|
||||||
{
|
|
||||||
indices.Append(logicalCellId[1] * this->PointDims[0] + logicalCellId[0]);
|
|
||||||
indices.Append(indices[0] + 1);
|
|
||||||
indices.Append(indices[1] + this->PointDims[0]);
|
|
||||||
indices.Append(indices[2] - 1);
|
|
||||||
cellShape = static_cast<vtkm::UInt8>(vtkm::CELL_SHAPE_QUAD);
|
|
||||||
numVerts = 4;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
private:
|
private:
|
||||||
vtkm::Id3 CellDims;
|
|
||||||
vtkm::Id3 PointDims;
|
|
||||||
bool Is3D = true;
|
|
||||||
};
|
|
||||||
|
|
||||||
class SingleCellTypeInterpolationHelper : public vtkm::exec::CellInterpolationHelper
|
|
||||||
{
|
|
||||||
using ConnType = vtkm::cont::ArrayHandle<vtkm::Id>;
|
|
||||||
using ConnPortalType = typename ConnType::ReadPortalType;
|
|
||||||
|
|
||||||
public:
|
|
||||||
SingleCellTypeInterpolationHelper() = default;
|
|
||||||
|
|
||||||
VTKM_CONT
|
|
||||||
SingleCellTypeInterpolationHelper(vtkm::UInt8 cellShape,
|
|
||||||
vtkm::IdComponent pointsPerCell,
|
|
||||||
const ConnType& connectivity,
|
|
||||||
vtkm::cont::DeviceAdapterId device,
|
|
||||||
vtkm::cont::Token& token)
|
|
||||||
: CellShape(cellShape)
|
|
||||||
, PointsPerCell(pointsPerCell)
|
|
||||||
, Connectivity(connectivity.PrepareForInput(device, token))
|
|
||||||
{
|
|
||||||
}
|
|
||||||
|
|
||||||
VTKM_EXEC
|
|
||||||
void GetCellInfo(const vtkm::Id& cellId,
|
|
||||||
vtkm::UInt8& cellShape,
|
|
||||||
vtkm::IdComponent& numVerts,
|
|
||||||
vtkm::VecVariable<vtkm::Id, 8>& indices) const override
|
|
||||||
{
|
|
||||||
cellShape = CellShape;
|
|
||||||
numVerts = PointsPerCell;
|
|
||||||
vtkm::Id n = static_cast<vtkm::Id>(PointsPerCell);
|
|
||||||
vtkm::Id offset = cellId * n;
|
|
||||||
|
|
||||||
for (vtkm::Id i = 0; i < n; i++)
|
|
||||||
indices.Append(Connectivity.Get(offset + i));
|
|
||||||
}
|
|
||||||
|
|
||||||
private:
|
|
||||||
vtkm::UInt8 CellShape;
|
|
||||||
vtkm::IdComponent PointsPerCell;
|
|
||||||
ConnPortalType Connectivity;
|
|
||||||
};
|
|
||||||
|
|
||||||
class ExplicitCellInterpolationHelper : public vtkm::exec::CellInterpolationHelper
|
|
||||||
{
|
|
||||||
using ShapeType = vtkm::cont::ArrayHandle<vtkm::UInt8>;
|
using ShapeType = vtkm::cont::ArrayHandle<vtkm::UInt8>;
|
||||||
using OffsetType = vtkm::cont::ArrayHandle<vtkm::Id>;
|
using OffsetType = vtkm::cont::ArrayHandle<vtkm::Id>;
|
||||||
using ConnType = vtkm::cont::ArrayHandle<vtkm::Id>;
|
using ConnType = vtkm::cont::ArrayHandle<vtkm::Id>;
|
||||||
|
|
||||||
using ShapePortalType = typename ShapeType::ReadPortalType;
|
using ShapePortalType = typename ShapeType::ReadPortalType;
|
||||||
using OffsetPortalType = typename OffsetType::ReadPortalType;
|
using OffsetPortalType = typename OffsetType::ReadPortalType;
|
||||||
using ConnPortalType = typename ConnType::ReadPortalType;
|
using ConnPortalType = typename ConnType::ReadPortalType;
|
||||||
|
|
||||||
public:
|
public:
|
||||||
ExplicitCellInterpolationHelper() = default;
|
enum class HelperType
|
||||||
|
{
|
||||||
|
STRUCTURED,
|
||||||
|
EXPSINGLE,
|
||||||
|
EXPLICIT
|
||||||
|
};
|
||||||
|
|
||||||
VTKM_CONT
|
VTKM_CONT
|
||||||
ExplicitCellInterpolationHelper(const ShapeType& shape,
|
CellInterpolationHelper() = default;
|
||||||
const OffsetType& offset,
|
|
||||||
const ConnType& connectivity,
|
VTKM_CONT
|
||||||
vtkm::cont::DeviceAdapterId device,
|
CellInterpolationHelper(const vtkm::Id3& cellDims, const vtkm::Id3& pointDims, bool is3D)
|
||||||
vtkm::cont::Token& token)
|
: CellDims(cellDims)
|
||||||
|
, PointDims(pointDims)
|
||||||
|
, Is3D(is3D)
|
||||||
|
{
|
||||||
|
this->Type = HelperType::STRUCTURED;
|
||||||
|
}
|
||||||
|
|
||||||
|
CellInterpolationHelper(const vtkm::UInt8 cellShape,
|
||||||
|
const vtkm::IdComponent pointsPerCell,
|
||||||
|
const ConnType& connectivity,
|
||||||
|
vtkm::cont::DeviceAdapterId device,
|
||||||
|
vtkm::cont::Token& token)
|
||||||
|
: CellShape(cellShape)
|
||||||
|
, PointsPerCell(pointsPerCell)
|
||||||
|
, Connectivity(connectivity.PrepareForInput(device, token))
|
||||||
|
{
|
||||||
|
this->Type = HelperType::EXPSINGLE;
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
VTKM_CONT
|
||||||
|
CellInterpolationHelper(const ShapeType& shape,
|
||||||
|
const OffsetType& offset,
|
||||||
|
const ConnType& connectivity,
|
||||||
|
vtkm::cont::DeviceAdapterId device,
|
||||||
|
vtkm::cont::Token& token)
|
||||||
: Shape(shape.PrepareForInput(device, token))
|
: Shape(shape.PrepareForInput(device, token))
|
||||||
, Offset(offset.PrepareForInput(device, token))
|
, Offset(offset.PrepareForInput(device, token))
|
||||||
, Connectivity(connectivity.PrepareForInput(device, token))
|
, Connectivity(connectivity.PrepareForInput(device, token))
|
||||||
{
|
{
|
||||||
|
this->Type = HelperType::EXPLICIT;
|
||||||
}
|
}
|
||||||
|
|
||||||
VTKM_EXEC
|
VTKM_EXEC
|
||||||
void GetCellInfo(const vtkm::Id& cellId,
|
void GetCellInfo(const vtkm::Id& cellId,
|
||||||
vtkm::UInt8& cellShape,
|
vtkm::UInt8& cellShape,
|
||||||
vtkm::IdComponent& numVerts,
|
vtkm::IdComponent& numVerts,
|
||||||
vtkm::VecVariable<vtkm::Id, 8>& indices) const override
|
vtkm::VecVariable<vtkm::Id, 8>& indices) const
|
||||||
{
|
{
|
||||||
cellShape = this->Shape.Get(cellId);
|
switch (this->Type)
|
||||||
const vtkm::Id offset = this->Offset.Get(cellId);
|
{
|
||||||
numVerts = static_cast<vtkm::IdComponent>(this->Offset.Get(cellId + 1) - offset);
|
case HelperType::STRUCTURED:
|
||||||
|
{
|
||||||
|
vtkm::Id3 logicalCellId;
|
||||||
|
logicalCellId[0] = cellId % this->CellDims[0];
|
||||||
|
logicalCellId[1] = (cellId / this->CellDims[0]) % this->CellDims[1];
|
||||||
|
if (this->Is3D)
|
||||||
|
{
|
||||||
|
logicalCellId[2] = cellId / (this->CellDims[0] * this->CellDims[1]);
|
||||||
|
indices.Append((logicalCellId[2] * this->PointDims[1] + logicalCellId[1]) *
|
||||||
|
this->PointDims[0] +
|
||||||
|
logicalCellId[0]);
|
||||||
|
indices.Append(indices[0] + 1);
|
||||||
|
indices.Append(indices[1] + this->PointDims[0]);
|
||||||
|
indices.Append(indices[2] - 1);
|
||||||
|
indices.Append(indices[0] + this->PointDims[0] * this->PointDims[1]);
|
||||||
|
indices.Append(indices[4] + 1);
|
||||||
|
indices.Append(indices[5] + this->PointDims[0]);
|
||||||
|
indices.Append(indices[6] - 1);
|
||||||
|
cellShape = static_cast<vtkm::UInt8>(vtkm::CELL_SHAPE_HEXAHEDRON);
|
||||||
|
numVerts = 8;
|
||||||
|
}
|
||||||
|
else
|
||||||
|
{
|
||||||
|
indices.Append(logicalCellId[1] * this->PointDims[0] + logicalCellId[0]);
|
||||||
|
indices.Append(indices[0] + 1);
|
||||||
|
indices.Append(indices[1] + this->PointDims[0]);
|
||||||
|
indices.Append(indices[2] - 1);
|
||||||
|
cellShape = static_cast<vtkm::UInt8>(vtkm::CELL_SHAPE_QUAD);
|
||||||
|
numVerts = 4;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
break;
|
||||||
|
|
||||||
for (vtkm::IdComponent i = 0; i < numVerts; i++)
|
case HelperType::EXPSINGLE:
|
||||||
indices.Append(this->Connectivity.Get(offset + i));
|
{
|
||||||
|
cellShape = this->CellShape;
|
||||||
|
numVerts = this->PointsPerCell;
|
||||||
|
vtkm::Id n = static_cast<vtkm::Id>(PointsPerCell);
|
||||||
|
vtkm::Id offset = cellId * n;
|
||||||
|
for (vtkm::Id i = 0; i < n; i++)
|
||||||
|
indices.Append(Connectivity.Get(offset + i));
|
||||||
|
}
|
||||||
|
break;
|
||||||
|
|
||||||
|
case HelperType::EXPLICIT:
|
||||||
|
{
|
||||||
|
cellShape = this->Shape.Get(cellId);
|
||||||
|
const vtkm::Id offset = this->Offset.Get(cellId);
|
||||||
|
numVerts = static_cast<vtkm::IdComponent>(this->Offset.Get(cellId + 1) - offset);
|
||||||
|
for (vtkm::IdComponent i = 0; i < numVerts; i++)
|
||||||
|
indices.Append(this->Connectivity.Get(offset + i));
|
||||||
|
}
|
||||||
|
break;
|
||||||
|
|
||||||
|
default:
|
||||||
|
{
|
||||||
|
// Code path not expected to execute in correct cases
|
||||||
|
// Supress unused variable warning
|
||||||
|
cellShape = vtkm::UInt8(0);
|
||||||
|
numVerts = vtkm::IdComponent(0);
|
||||||
|
}
|
||||||
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
private:
|
private:
|
||||||
|
HelperType Type;
|
||||||
|
// variables for structured type
|
||||||
|
vtkm::Id3 CellDims;
|
||||||
|
vtkm::Id3 PointDims;
|
||||||
|
bool Is3D = true;
|
||||||
|
// variables for single explicit type
|
||||||
|
vtkm::UInt8 CellShape;
|
||||||
|
vtkm::IdComponent PointsPerCell;
|
||||||
|
// variables for explicit type
|
||||||
ShapePortalType Shape;
|
ShapePortalType Shape;
|
||||||
OffsetPortalType Offset;
|
OffsetPortalType Offset;
|
||||||
ConnPortalType Connectivity;
|
ConnPortalType Connectivity;
|
||||||
@ -198,26 +184,19 @@ namespace cont
|
|||||||
|
|
||||||
class CellInterpolationHelper : public vtkm::cont::ExecutionObjectBase
|
class CellInterpolationHelper : public vtkm::cont::ExecutionObjectBase
|
||||||
{
|
{
|
||||||
public:
|
private:
|
||||||
using HandleType = vtkm::cont::VirtualObjectHandle<vtkm::exec::CellInterpolationHelper>;
|
using ExecutionType = vtkm::exec::CellInterpolationHelper;
|
||||||
|
|
||||||
virtual ~CellInterpolationHelper() = default;
|
|
||||||
|
|
||||||
VTKM_CONT virtual const vtkm::exec::CellInterpolationHelper* PrepareForExecution(
|
|
||||||
vtkm::cont::DeviceAdapterId device,
|
|
||||||
vtkm::cont::Token& token) const = 0;
|
|
||||||
};
|
|
||||||
|
|
||||||
class StructuredCellInterpolationHelper : public vtkm::cont::CellInterpolationHelper
|
|
||||||
{
|
|
||||||
public:
|
|
||||||
using Structured2DType = vtkm::cont::CellSetStructured<2>;
|
using Structured2DType = vtkm::cont::CellSetStructured<2>;
|
||||||
using Structured3DType = vtkm::cont::CellSetStructured<3>;
|
using Structured3DType = vtkm::cont::CellSetStructured<3>;
|
||||||
|
using SingleExplicitType = vtkm::cont::CellSetSingleType<>;
|
||||||
|
using ExplicitType = vtkm::cont::CellSetExplicit<>;
|
||||||
|
|
||||||
StructuredCellInterpolationHelper() = default;
|
public:
|
||||||
|
VTKM_CONT
|
||||||
|
CellInterpolationHelper() = default;
|
||||||
|
|
||||||
VTKM_CONT
|
VTKM_CONT
|
||||||
StructuredCellInterpolationHelper(const vtkm::cont::DynamicCellSet& cellSet)
|
CellInterpolationHelper(const vtkm::cont::DynamicCellSet& cellSet)
|
||||||
{
|
{
|
||||||
if (cellSet.IsSameType(Structured2DType()))
|
if (cellSet.IsSameType(Structured2DType()))
|
||||||
{
|
{
|
||||||
@ -228,6 +207,7 @@ public:
|
|||||||
cellSet.Cast<Structured2DType>().GetSchedulingRange(vtkm::TopologyElementTagPoint());
|
cellSet.Cast<Structured2DType>().GetSchedulingRange(vtkm::TopologyElementTagPoint());
|
||||||
this->CellDims = vtkm::Id3(cellDims[0], cellDims[1], 0);
|
this->CellDims = vtkm::Id3(cellDims[0], cellDims[1], 0);
|
||||||
this->PointDims = vtkm::Id3(pointDims[0], pointDims[1], 1);
|
this->PointDims = vtkm::Id3(pointDims[0], pointDims[1], 1);
|
||||||
|
this->Type = vtkm::exec::CellInterpolationHelper::HelperType::STRUCTURED;
|
||||||
}
|
}
|
||||||
else if (cellSet.IsSameType(Structured3DType()))
|
else if (cellSet.IsSameType(Structured3DType()))
|
||||||
{
|
{
|
||||||
@ -236,113 +216,22 @@ public:
|
|||||||
cellSet.Cast<Structured3DType>().GetSchedulingRange(vtkm::TopologyElementTagCell());
|
cellSet.Cast<Structured3DType>().GetSchedulingRange(vtkm::TopologyElementTagCell());
|
||||||
this->PointDims =
|
this->PointDims =
|
||||||
cellSet.Cast<Structured3DType>().GetSchedulingRange(vtkm::TopologyElementTagPoint());
|
cellSet.Cast<Structured3DType>().GetSchedulingRange(vtkm::TopologyElementTagPoint());
|
||||||
|
this->Type = vtkm::exec::CellInterpolationHelper::HelperType::STRUCTURED;
|
||||||
}
|
}
|
||||||
else
|
else if (cellSet.IsSameType(SingleExplicitType()))
|
||||||
throw vtkm::cont::ErrorBadType("Cell set is not of type CellSetStructured");
|
|
||||||
}
|
|
||||||
|
|
||||||
VTKM_CONT
|
|
||||||
const vtkm::exec::CellInterpolationHelper* PrepareForExecution(
|
|
||||||
vtkm::cont::DeviceAdapterId deviceId,
|
|
||||||
vtkm::cont::Token& token) const override
|
|
||||||
{
|
|
||||||
auto& tracker = vtkm::cont::GetRuntimeDeviceTracker();
|
|
||||||
const bool valid = tracker.CanRunOn(deviceId);
|
|
||||||
if (!valid)
|
|
||||||
{
|
|
||||||
throwFailedRuntimeDeviceTransfer("StructuredCellInterpolationHelper", deviceId);
|
|
||||||
}
|
|
||||||
|
|
||||||
using ExecutionType = vtkm::exec::StructuredCellInterpolationHelper;
|
|
||||||
ExecutionType* execObject = new ExecutionType(this->CellDims, this->PointDims, this->Is3D);
|
|
||||||
this->ExecHandle.Reset(execObject);
|
|
||||||
|
|
||||||
return this->ExecHandle.PrepareForExecution(deviceId, token);
|
|
||||||
}
|
|
||||||
|
|
||||||
private:
|
|
||||||
vtkm::Id3 CellDims;
|
|
||||||
vtkm::Id3 PointDims;
|
|
||||||
bool Is3D = true;
|
|
||||||
mutable HandleType ExecHandle;
|
|
||||||
};
|
|
||||||
|
|
||||||
class SingleCellTypeInterpolationHelper : public vtkm::cont::CellInterpolationHelper
|
|
||||||
{
|
|
||||||
public:
|
|
||||||
using SingleExplicitType = vtkm::cont::CellSetSingleType<>;
|
|
||||||
|
|
||||||
SingleCellTypeInterpolationHelper() = default;
|
|
||||||
|
|
||||||
VTKM_CONT
|
|
||||||
SingleCellTypeInterpolationHelper(const vtkm::cont::DynamicCellSet& cellSet)
|
|
||||||
{
|
|
||||||
if (cellSet.IsSameType(SingleExplicitType()))
|
|
||||||
{
|
{
|
||||||
SingleExplicitType CellSet = cellSet.Cast<SingleExplicitType>();
|
SingleExplicitType CellSet = cellSet.Cast<SingleExplicitType>();
|
||||||
|
|
||||||
const auto cellShapes =
|
const auto cellShapes =
|
||||||
CellSet.GetShapesArray(vtkm::TopologyElementTagCell(), vtkm::TopologyElementTagPoint());
|
CellSet.GetShapesArray(vtkm::TopologyElementTagCell(), vtkm::TopologyElementTagPoint());
|
||||||
const auto numIndices =
|
const auto numIndices =
|
||||||
CellSet.GetNumIndicesArray(vtkm::TopologyElementTagCell(), vtkm::TopologyElementTagPoint());
|
CellSet.GetNumIndicesArray(vtkm::TopologyElementTagCell(), vtkm::TopologyElementTagPoint());
|
||||||
|
|
||||||
CellShape = vtkm::cont::ArrayGetValue(0, cellShapes);
|
CellShape = vtkm::cont::ArrayGetValue(0, cellShapes);
|
||||||
PointsPerCell = vtkm::cont::ArrayGetValue(0, numIndices);
|
PointsPerCell = vtkm::cont::ArrayGetValue(0, numIndices);
|
||||||
Connectivity = CellSet.GetConnectivityArray(vtkm::TopologyElementTagCell(),
|
Connectivity = CellSet.GetConnectivityArray(vtkm::TopologyElementTagCell(),
|
||||||
vtkm::TopologyElementTagPoint());
|
vtkm::TopologyElementTagPoint());
|
||||||
|
this->Type = vtkm::exec::CellInterpolationHelper::HelperType::EXPSINGLE;
|
||||||
}
|
}
|
||||||
else
|
else if (cellSet.IsSameType(ExplicitType()))
|
||||||
throw vtkm::cont::ErrorBadType("Cell set is not of type CellSetSingleType");
|
|
||||||
}
|
|
||||||
|
|
||||||
struct SingleCellTypeFunctor
|
|
||||||
{
|
|
||||||
VTKM_CONT bool operator()(vtkm::cont::DeviceAdapterId device,
|
|
||||||
const vtkm::cont::SingleCellTypeInterpolationHelper& contInterpolator,
|
|
||||||
HandleType& execInterpolator,
|
|
||||||
vtkm::cont::Token& token) const
|
|
||||||
{
|
|
||||||
using ExecutionType = vtkm::exec::SingleCellTypeInterpolationHelper;
|
|
||||||
ExecutionType* execObject = new ExecutionType(contInterpolator.CellShape,
|
|
||||||
contInterpolator.PointsPerCell,
|
|
||||||
contInterpolator.Connectivity,
|
|
||||||
device,
|
|
||||||
token);
|
|
||||||
execInterpolator.Reset(execObject);
|
|
||||||
return true;
|
|
||||||
}
|
|
||||||
};
|
|
||||||
|
|
||||||
VTKM_CONT
|
|
||||||
const vtkm::exec::CellInterpolationHelper* PrepareForExecution(
|
|
||||||
vtkm::cont::DeviceAdapterId deviceId,
|
|
||||||
vtkm::cont::Token& token) const override
|
|
||||||
{
|
|
||||||
const bool success = vtkm::cont::TryExecuteOnDevice(
|
|
||||||
deviceId, SingleCellTypeFunctor(), *this, this->ExecHandle, token);
|
|
||||||
if (!success)
|
|
||||||
{
|
|
||||||
throwFailedRuntimeDeviceTransfer("SingleCellTypeInterpolationHelper", deviceId);
|
|
||||||
}
|
|
||||||
return this->ExecHandle.PrepareForExecution(deviceId, token);
|
|
||||||
}
|
|
||||||
|
|
||||||
private:
|
|
||||||
vtkm::UInt8 CellShape;
|
|
||||||
vtkm::IdComponent PointsPerCell;
|
|
||||||
vtkm::cont::ArrayHandle<vtkm::Id> Connectivity;
|
|
||||||
mutable HandleType ExecHandle;
|
|
||||||
};
|
|
||||||
|
|
||||||
class ExplicitCellInterpolationHelper : public vtkm::cont::CellInterpolationHelper
|
|
||||||
{
|
|
||||||
public:
|
|
||||||
ExplicitCellInterpolationHelper() = default;
|
|
||||||
|
|
||||||
VTKM_CONT
|
|
||||||
ExplicitCellInterpolationHelper(const vtkm::cont::DynamicCellSet& cellSet)
|
|
||||||
{
|
|
||||||
if (cellSet.IsSameType(vtkm::cont::CellSetExplicit<>()))
|
|
||||||
{
|
{
|
||||||
vtkm::cont::CellSetExplicit<> CellSet = cellSet.Cast<vtkm::cont::CellSetExplicit<>>();
|
vtkm::cont::CellSetExplicit<> CellSet = cellSet.Cast<vtkm::cont::CellSetExplicit<>>();
|
||||||
Shape =
|
Shape =
|
||||||
@ -351,28 +240,44 @@ public:
|
|||||||
CellSet.GetOffsetsArray(vtkm::TopologyElementTagCell(), vtkm::TopologyElementTagPoint());
|
CellSet.GetOffsetsArray(vtkm::TopologyElementTagCell(), vtkm::TopologyElementTagPoint());
|
||||||
Connectivity = CellSet.GetConnectivityArray(vtkm::TopologyElementTagCell(),
|
Connectivity = CellSet.GetConnectivityArray(vtkm::TopologyElementTagCell(),
|
||||||
vtkm::TopologyElementTagPoint());
|
vtkm::TopologyElementTagPoint());
|
||||||
|
this->Type = vtkm::exec::CellInterpolationHelper::HelperType::EXPLICIT;
|
||||||
}
|
}
|
||||||
else
|
else
|
||||||
throw vtkm::cont::ErrorBadType("Cell set is not of type CellSetExplicit");
|
throw vtkm::cont::ErrorInternal("Unsupported cellset type");
|
||||||
}
|
}
|
||||||
|
|
||||||
VTKM_CONT
|
VTKM_CONT
|
||||||
const vtkm::exec::CellInterpolationHelper* PrepareForExecution(
|
const ExecutionType PrepareForExecution(vtkm::cont::DeviceAdapterId device,
|
||||||
vtkm::cont::DeviceAdapterId device,
|
vtkm::cont::Token& token) const
|
||||||
vtkm::cont::Token& token) const override
|
|
||||||
{
|
{
|
||||||
using ExecutionType = vtkm::exec::ExplicitCellInterpolationHelper;
|
switch (this->Type)
|
||||||
ExecutionType* execObject =
|
{
|
||||||
new ExecutionType(this->Shape, this->Offset, this->Connectivity, device, token);
|
case ExecutionType::HelperType::STRUCTURED:
|
||||||
this->ExecHandle.Reset(execObject);
|
return ExecutionType(this->CellDims, this->PointDims, this->Is3D);
|
||||||
return this->ExecHandle.PrepareForExecution(device, token);
|
|
||||||
|
case ExecutionType::HelperType::EXPSINGLE:
|
||||||
|
return ExecutionType(
|
||||||
|
this->CellShape, this->PointsPerCell, this->Connectivity, device, token);
|
||||||
|
|
||||||
|
case ExecutionType::HelperType::EXPLICIT:
|
||||||
|
return ExecutionType(this->Shape, this->Offset, this->Connectivity, device, token);
|
||||||
|
}
|
||||||
|
throw vtkm::cont::ErrorInternal("Undefined case for building cell interpolation helper");
|
||||||
}
|
}
|
||||||
|
|
||||||
private:
|
private:
|
||||||
|
// Variables required for strucutred grids
|
||||||
|
vtkm::Id3 CellDims;
|
||||||
|
vtkm::Id3 PointDims;
|
||||||
|
bool Is3D = true;
|
||||||
|
// variables for single explicit type
|
||||||
|
vtkm::UInt8 CellShape;
|
||||||
|
vtkm::IdComponent PointsPerCell;
|
||||||
|
// Variables required for unstructured grids
|
||||||
vtkm::cont::ArrayHandle<vtkm::UInt8> Shape;
|
vtkm::cont::ArrayHandle<vtkm::UInt8> Shape;
|
||||||
vtkm::cont::ArrayHandle<vtkm::Id> Offset;
|
vtkm::cont::ArrayHandle<vtkm::Id> Offset;
|
||||||
vtkm::cont::ArrayHandle<vtkm::Id> Connectivity;
|
vtkm::cont::ArrayHandle<vtkm::Id> Connectivity;
|
||||||
mutable HandleType ExecHandle;
|
ExecutionType::HelperType Type;
|
||||||
};
|
};
|
||||||
|
|
||||||
} //namespace cont
|
} //namespace cont
|
||||||
|
@ -13,8 +13,6 @@
|
|||||||
#ifndef vtk_m_worklet_particleadvection_EulerIntegrator_h
|
#ifndef vtk_m_worklet_particleadvection_EulerIntegrator_h
|
||||||
#define vtk_m_worklet_particleadvection_EulerIntegrator_h
|
#define vtk_m_worklet_particleadvection_EulerIntegrator_h
|
||||||
|
|
||||||
#include <vtkm/worklet/particleadvection/IntegratorBase.h>
|
|
||||||
|
|
||||||
namespace vtkm
|
namespace vtkm
|
||||||
{
|
{
|
||||||
namespace worklet
|
namespace worklet
|
||||||
@ -22,73 +20,56 @@ namespace worklet
|
|||||||
namespace particleadvection
|
namespace particleadvection
|
||||||
{
|
{
|
||||||
|
|
||||||
template <typename FieldEvaluateType>
|
template <typename EvaluatorType>
|
||||||
class EulerIntegrator : public IntegratorBase
|
class ExecEulerIntegrator
|
||||||
{
|
{
|
||||||
|
public:
|
||||||
|
VTKM_EXEC_CONT
|
||||||
|
ExecEulerIntegrator(const EvaluatorType& evaluator)
|
||||||
|
: Evaluator(evaluator)
|
||||||
|
{
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename Particle>
|
||||||
|
VTKM_EXEC IntegratorStatus CheckStep(Particle& particle,
|
||||||
|
vtkm::FloatDefault stepLength,
|
||||||
|
vtkm::Vec3f& velocity) const
|
||||||
|
{
|
||||||
|
auto time = particle.Time;
|
||||||
|
auto inpos = particle.Pos;
|
||||||
|
vtkm::VecVariable<vtkm::Vec3f, 2> vectors;
|
||||||
|
GridEvaluatorStatus status = this->Evaluator.Evaluate(inpos, time, vectors);
|
||||||
|
if (status.CheckOk())
|
||||||
|
velocity = particle.Velocity(vectors, stepLength);
|
||||||
|
return IntegratorStatus(status);
|
||||||
|
}
|
||||||
|
|
||||||
|
private:
|
||||||
|
EvaluatorType Evaluator;
|
||||||
|
};
|
||||||
|
|
||||||
|
template <typename EvaluatorType>
|
||||||
|
class EulerIntegrator
|
||||||
|
{
|
||||||
|
private:
|
||||||
|
EvaluatorType Evaluator;
|
||||||
|
|
||||||
public:
|
public:
|
||||||
EulerIntegrator() = default;
|
EulerIntegrator() = default;
|
||||||
|
|
||||||
VTKM_CONT
|
VTKM_CONT
|
||||||
EulerIntegrator(const FieldEvaluateType& evaluator, const vtkm::FloatDefault stepLength)
|
EulerIntegrator(const EvaluatorType& evaluator)
|
||||||
: IntegratorBase(stepLength)
|
: Evaluator(evaluator)
|
||||||
, Evaluator(evaluator)
|
|
||||||
{
|
{
|
||||||
}
|
}
|
||||||
|
|
||||||
template <typename Device>
|
VTKM_CONT auto PrepareForExecution(vtkm::cont::DeviceAdapterId device,
|
||||||
class ExecObject
|
vtkm::cont::Token& token) const
|
||||||
: public IntegratorBase::ExecObjectBaseImpl<
|
-> ExecEulerIntegrator<decltype(this->Evaluator.PrepareForExecution(device, token))>
|
||||||
vtkm::cont::internal::ExecutionObjectType<FieldEvaluateType, Device>,
|
|
||||||
typename EulerIntegrator::template ExecObject<Device>>
|
|
||||||
{
|
{
|
||||||
VTKM_IS_DEVICE_ADAPTER_TAG(Device);
|
auto evaluator = this->Evaluator.PrepareForExecution(device, token);
|
||||||
|
using ExecEvaluatorType = decltype(evaluator);
|
||||||
using FieldEvaluateExecType =
|
return ExecEulerIntegrator<ExecEvaluatorType>(evaluator);
|
||||||
vtkm::cont::internal::ExecutionObjectType<FieldEvaluateType, Device>;
|
|
||||||
using Superclass =
|
|
||||||
IntegratorBase::ExecObjectBaseImpl<FieldEvaluateExecType,
|
|
||||||
typename EulerIntegrator::template ExecObject<Device>>;
|
|
||||||
|
|
||||||
public:
|
|
||||||
VTKM_EXEC_CONT
|
|
||||||
ExecObject(const FieldEvaluateExecType& evaluator,
|
|
||||||
vtkm::FloatDefault stepLength,
|
|
||||||
vtkm::FloatDefault tolerance)
|
|
||||||
: Superclass(evaluator, stepLength, tolerance)
|
|
||||||
{
|
|
||||||
}
|
|
||||||
|
|
||||||
VTKM_EXEC
|
|
||||||
IntegratorStatus CheckStep(vtkm::Particle* particle,
|
|
||||||
vtkm::FloatDefault stepLength,
|
|
||||||
vtkm::Vec3f& velocity) const
|
|
||||||
{
|
|
||||||
auto time = particle->Time;
|
|
||||||
auto inpos = particle->Pos;
|
|
||||||
vtkm::VecVariable<vtkm::Vec3f, 2> vectors;
|
|
||||||
GridEvaluatorStatus status = this->Evaluator.Evaluate(inpos, time, vectors);
|
|
||||||
if (status.CheckOk())
|
|
||||||
velocity = particle->Velocity(vectors, stepLength);
|
|
||||||
return IntegratorStatus(status);
|
|
||||||
}
|
|
||||||
};
|
|
||||||
|
|
||||||
private:
|
|
||||||
FieldEvaluateType Evaluator;
|
|
||||||
|
|
||||||
protected:
|
|
||||||
VTKM_CONT virtual void PrepareForExecutionImpl(
|
|
||||||
vtkm::cont::DeviceAdapterId device,
|
|
||||||
vtkm::cont::VirtualObjectHandle<IntegratorBase::ExecObject>& execObjectHandle,
|
|
||||||
vtkm::cont::Token& token) const override
|
|
||||||
{
|
|
||||||
vtkm::cont::TryExecuteOnDevice(device,
|
|
||||||
detail::IntegratorPrepareForExecutionFunctor<ExecObject>(),
|
|
||||||
execObjectHandle,
|
|
||||||
this->Evaluator,
|
|
||||||
this->StepLength,
|
|
||||||
this->Tolerance,
|
|
||||||
token);
|
|
||||||
}
|
}
|
||||||
}; //EulerIntegrator
|
}; //EulerIntegrator
|
||||||
|
|
||||||
|
@ -26,22 +26,8 @@ namespace worklet
|
|||||||
namespace particleadvection
|
namespace particleadvection
|
||||||
{
|
{
|
||||||
|
|
||||||
class ExecutionField : public vtkm::VirtualObjectBase
|
|
||||||
{
|
|
||||||
public:
|
|
||||||
VTKM_EXEC_CONT
|
|
||||||
virtual ~ExecutionField() noexcept override {}
|
|
||||||
|
|
||||||
VTKM_EXEC
|
|
||||||
virtual void GetValue(const vtkm::VecVariable<vtkm::Id, 8>& indices,
|
|
||||||
const vtkm::Id vertices,
|
|
||||||
const vtkm::Vec3f& parametric,
|
|
||||||
const vtkm::UInt8 cellShape,
|
|
||||||
vtkm::VecVariable<vtkm::Vec3f, 2>& value) const = 0;
|
|
||||||
};
|
|
||||||
|
|
||||||
template <typename FieldArrayType>
|
template <typename FieldArrayType>
|
||||||
class ExecutionVelocityField : public vtkm::worklet::particleadvection::ExecutionField
|
class ExecutionVelocityField
|
||||||
{
|
{
|
||||||
public:
|
public:
|
||||||
using FieldPortalType = typename FieldArrayType::ReadPortalType;
|
using FieldPortalType = typename FieldArrayType::ReadPortalType;
|
||||||
@ -73,7 +59,7 @@ private:
|
|||||||
};
|
};
|
||||||
|
|
||||||
template <typename FieldArrayType>
|
template <typename FieldArrayType>
|
||||||
class ExecutionElectroMagneticField : public vtkm::worklet::particleadvection::ExecutionField
|
class ExecutionElectroMagneticField
|
||||||
{
|
{
|
||||||
public:
|
public:
|
||||||
using FieldPortalType = typename FieldArrayType::ReadPortalType;
|
using FieldPortalType = typename FieldArrayType::ReadPortalType;
|
||||||
@ -112,22 +98,15 @@ private:
|
|||||||
FieldPortalType MagneticValues;
|
FieldPortalType MagneticValues;
|
||||||
};
|
};
|
||||||
|
|
||||||
class Field : public vtkm::cont::ExecutionObjectBase
|
template <typename FieldArrayType>
|
||||||
|
class VelocityField : public vtkm::cont::ExecutionObjectBase
|
||||||
{
|
{
|
||||||
public:
|
public:
|
||||||
using HandleType = vtkm::cont::VirtualObjectHandle<ExecutionField>;
|
using ExecutionType = ExecutionVelocityField<FieldArrayType>;
|
||||||
|
|
||||||
virtual ~Field() = default;
|
|
||||||
|
|
||||||
VTKM_CONT
|
VTKM_CONT
|
||||||
virtual const ExecutionField* PrepareForExecution(vtkm::cont::DeviceAdapterId deviceId,
|
VelocityField() = default;
|
||||||
vtkm::cont::Token& token) const = 0;
|
|
||||||
};
|
|
||||||
|
|
||||||
template <typename FieldArrayType>
|
|
||||||
class VelocityField : public vtkm::worklet::particleadvection::Field
|
|
||||||
{
|
|
||||||
public:
|
|
||||||
VTKM_CONT
|
VTKM_CONT
|
||||||
VelocityField(const FieldArrayType& fieldValues)
|
VelocityField(const FieldArrayType& fieldValues)
|
||||||
: FieldValues(fieldValues)
|
: FieldValues(fieldValues)
|
||||||
@ -135,24 +114,25 @@ public:
|
|||||||
}
|
}
|
||||||
|
|
||||||
VTKM_CONT
|
VTKM_CONT
|
||||||
const ExecutionField* PrepareForExecution(vtkm::cont::DeviceAdapterId device,
|
const ExecutionType PrepareForExecution(vtkm::cont::DeviceAdapterId device,
|
||||||
vtkm::cont::Token& token) const override
|
vtkm::cont::Token& token) const
|
||||||
{
|
{
|
||||||
using ExecutionType = ExecutionVelocityField<FieldArrayType>;
|
return ExecutionType(this->FieldValues, device, token);
|
||||||
ExecutionType* execObject = new ExecutionType(this->FieldValues, device, token);
|
|
||||||
this->ExecHandle.Reset(execObject);
|
|
||||||
return this->ExecHandle.PrepareForExecution(device, token);
|
|
||||||
}
|
}
|
||||||
|
|
||||||
private:
|
private:
|
||||||
FieldArrayType FieldValues;
|
FieldArrayType FieldValues;
|
||||||
mutable HandleType ExecHandle;
|
|
||||||
};
|
};
|
||||||
|
|
||||||
template <typename FieldArrayType>
|
template <typename FieldArrayType>
|
||||||
class ElectroMagneticField : public vtkm::worklet::particleadvection::Field
|
class ElectroMagneticField : public vtkm::cont::ExecutionObjectBase
|
||||||
{
|
{
|
||||||
public:
|
public:
|
||||||
|
using ExecutionType = ExecutionElectroMagneticField<FieldArrayType>;
|
||||||
|
|
||||||
|
VTKM_CONT
|
||||||
|
ElectroMagneticField() = default;
|
||||||
|
|
||||||
VTKM_CONT
|
VTKM_CONT
|
||||||
ElectroMagneticField(const FieldArrayType& electricField, const FieldArrayType& magneticField)
|
ElectroMagneticField(const FieldArrayType& electricField, const FieldArrayType& magneticField)
|
||||||
: ElectricField(electricField)
|
: ElectricField(electricField)
|
||||||
@ -161,20 +141,15 @@ public:
|
|||||||
}
|
}
|
||||||
|
|
||||||
VTKM_CONT
|
VTKM_CONT
|
||||||
const ExecutionField* PrepareForExecution(vtkm::cont::DeviceAdapterId device,
|
const ExecutionType PrepareForExecution(vtkm::cont::DeviceAdapterId device,
|
||||||
vtkm::cont::Token& token) const override
|
vtkm::cont::Token& token) const
|
||||||
{
|
{
|
||||||
using ExecutionType = ExecutionElectroMagneticField<FieldArrayType>;
|
return ExecutionType(this->ElectricField, this->MagneticField, device, token);
|
||||||
ExecutionType* execObject =
|
|
||||||
new ExecutionType(this->ElectricField, this->MagneticField, device, token);
|
|
||||||
this->ExecHandle.Reset(execObject);
|
|
||||||
return this->ExecHandle.PrepareForExecution(device, token);
|
|
||||||
}
|
}
|
||||||
|
|
||||||
private:
|
private:
|
||||||
FieldArrayType ElectricField;
|
FieldArrayType ElectricField;
|
||||||
FieldArrayType MagneticField;
|
FieldArrayType MagneticField;
|
||||||
mutable HandleType ExecHandle;
|
|
||||||
};
|
};
|
||||||
|
|
||||||
} // namespace particleadvection
|
} // namespace particleadvection
|
||||||
|
@ -22,9 +22,6 @@
|
|||||||
#include <vtkm/cont/DataSet.h>
|
#include <vtkm/cont/DataSet.h>
|
||||||
#include <vtkm/cont/DeviceAdapter.h>
|
#include <vtkm/cont/DeviceAdapter.h>
|
||||||
|
|
||||||
#include <vtkm/worklet/particleadvection/CellInterpolationHelper.h>
|
|
||||||
#include <vtkm/worklet/particleadvection/IntegratorBase.h>
|
|
||||||
|
|
||||||
namespace vtkm
|
namespace vtkm
|
||||||
{
|
{
|
||||||
namespace worklet
|
namespace worklet
|
||||||
|
@ -26,7 +26,6 @@
|
|||||||
#include <vtkm/worklet/particleadvection/CellInterpolationHelper.h>
|
#include <vtkm/worklet/particleadvection/CellInterpolationHelper.h>
|
||||||
#include <vtkm/worklet/particleadvection/Field.h>
|
#include <vtkm/worklet/particleadvection/Field.h>
|
||||||
#include <vtkm/worklet/particleadvection/GridEvaluatorStatus.h>
|
#include <vtkm/worklet/particleadvection/GridEvaluatorStatus.h>
|
||||||
#include <vtkm/worklet/particleadvection/IntegratorBase.h>
|
|
||||||
|
|
||||||
namespace vtkm
|
namespace vtkm
|
||||||
{
|
{
|
||||||
@ -46,7 +45,7 @@ public:
|
|||||||
|
|
||||||
VTKM_CONT
|
VTKM_CONT
|
||||||
ExecutionGridEvaluator(const vtkm::cont::CellLocatorGeneral& locator,
|
ExecutionGridEvaluator(const vtkm::cont::CellLocatorGeneral& locator,
|
||||||
std::shared_ptr<vtkm::cont::CellInterpolationHelper> interpolationHelper,
|
const vtkm::cont::CellInterpolationHelper interpolationHelper,
|
||||||
const vtkm::Bounds& bounds,
|
const vtkm::Bounds& bounds,
|
||||||
const FieldType& field,
|
const FieldType& field,
|
||||||
const GhostCellArrayType& ghostCells,
|
const GhostCellArrayType& ghostCells,
|
||||||
@ -56,7 +55,7 @@ public:
|
|||||||
, Field(field.PrepareForExecution(device, token))
|
, Field(field.PrepareForExecution(device, token))
|
||||||
, GhostCells(ghostCells.PrepareForInput(device, token))
|
, GhostCells(ghostCells.PrepareForInput(device, token))
|
||||||
, HaveGhostCells(ghostCells.GetNumberOfValues() > 0)
|
, HaveGhostCells(ghostCells.GetNumberOfValues() > 0)
|
||||||
, InterpolationHelper(interpolationHelper->PrepareForExecution(device, token))
|
, InterpolationHelper(interpolationHelper.PrepareForExecution(device, token))
|
||||||
, Locator(locator.PrepareForExecution(device, token))
|
, Locator(locator.PrepareForExecution(device, token))
|
||||||
{
|
{
|
||||||
}
|
}
|
||||||
@ -125,9 +124,8 @@ public:
|
|||||||
vtkm::IdComponent nVerts;
|
vtkm::IdComponent nVerts;
|
||||||
vtkm::VecVariable<vtkm::Id, 8> ptIndices;
|
vtkm::VecVariable<vtkm::Id, 8> ptIndices;
|
||||||
vtkm::VecVariable<vtkm::Vec3f, 8> fieldValues;
|
vtkm::VecVariable<vtkm::Vec3f, 8> fieldValues;
|
||||||
InterpolationHelper->GetCellInfo(cellId, cellShape, nVerts, ptIndices);
|
this->InterpolationHelper.GetCellInfo(cellId, cellShape, nVerts, ptIndices);
|
||||||
|
this->Field.GetValue(ptIndices, nVerts, parametric, cellShape, out);
|
||||||
this->Field->GetValue(ptIndices, nVerts, parametric, cellShape, out);
|
|
||||||
status.SetOk();
|
status.SetOk();
|
||||||
}
|
}
|
||||||
|
|
||||||
@ -146,10 +144,10 @@ private:
|
|||||||
using GhostCellPortal = typename vtkm::cont::ArrayHandle<vtkm::UInt8>::ReadPortalType;
|
using GhostCellPortal = typename vtkm::cont::ArrayHandle<vtkm::UInt8>::ReadPortalType;
|
||||||
|
|
||||||
vtkm::Bounds Bounds;
|
vtkm::Bounds Bounds;
|
||||||
const vtkm::worklet::particleadvection::ExecutionField* Field;
|
typename FieldType::ExecutionType Field;
|
||||||
GhostCellPortal GhostCells;
|
GhostCellPortal GhostCells;
|
||||||
bool HaveGhostCells;
|
bool HaveGhostCells;
|
||||||
const vtkm::exec::CellInterpolationHelper* InterpolationHelper;
|
vtkm::exec::CellInterpolationHelper InterpolationHelper;
|
||||||
typename vtkm::cont::CellLocatorGeneral::ExecObjType Locator;
|
typename vtkm::cont::CellLocatorGeneral::ExecObjType Locator;
|
||||||
};
|
};
|
||||||
|
|
||||||
@ -217,32 +215,13 @@ private:
|
|||||||
this->Locator.SetCoordinates(coordinates);
|
this->Locator.SetCoordinates(coordinates);
|
||||||
this->Locator.SetCellSet(cellset);
|
this->Locator.SetCellSet(cellset);
|
||||||
this->Locator.Update();
|
this->Locator.Update();
|
||||||
if (cellset.IsSameType(Structured2DType()) || cellset.IsSameType(Structured3DType()))
|
this->InterpolationHelper = vtkm::cont::CellInterpolationHelper(cellset);
|
||||||
{
|
|
||||||
vtkm::cont::StructuredCellInterpolationHelper interpolationHelper(cellset);
|
|
||||||
this->InterpolationHelper =
|
|
||||||
std::make_shared<vtkm::cont::StructuredCellInterpolationHelper>(interpolationHelper);
|
|
||||||
}
|
|
||||||
else if (cellset.IsSameType(vtkm::cont::CellSetSingleType<>()))
|
|
||||||
{
|
|
||||||
vtkm::cont::SingleCellTypeInterpolationHelper interpolationHelper(cellset);
|
|
||||||
this->InterpolationHelper =
|
|
||||||
std::make_shared<vtkm::cont::SingleCellTypeInterpolationHelper>(interpolationHelper);
|
|
||||||
}
|
|
||||||
else if (cellset.IsSameType(vtkm::cont::CellSetExplicit<>()))
|
|
||||||
{
|
|
||||||
vtkm::cont::ExplicitCellInterpolationHelper interpolationHelper(cellset);
|
|
||||||
this->InterpolationHelper =
|
|
||||||
std::make_shared<vtkm::cont::ExplicitCellInterpolationHelper>(interpolationHelper);
|
|
||||||
}
|
|
||||||
else
|
|
||||||
throw vtkm::cont::ErrorInternal("Unsupported cellset type.");
|
|
||||||
}
|
}
|
||||||
|
|
||||||
vtkm::Bounds Bounds;
|
vtkm::Bounds Bounds;
|
||||||
FieldType Field;
|
FieldType Field;
|
||||||
GhostCellArrayType GhostCellArray;
|
GhostCellArrayType GhostCellArray;
|
||||||
std::shared_ptr<vtkm::cont::CellInterpolationHelper> InterpolationHelper;
|
vtkm::cont::CellInterpolationHelper InterpolationHelper;
|
||||||
vtkm::cont::CellLocatorGeneral Locator;
|
vtkm::cont::CellLocatorGeneral Locator;
|
||||||
};
|
};
|
||||||
|
|
||||||
|
@ -1,252 +0,0 @@
|
|||||||
//=============================================================================
|
|
||||||
//
|
|
||||||
// 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_particleadvection_IntegratorBase_h
|
|
||||||
#define vtk_m_worklet_particleadvection_IntegratorBase_h
|
|
||||||
|
|
||||||
#include <limits>
|
|
||||||
|
|
||||||
#include <vtkm/Bitset.h>
|
|
||||||
#include <vtkm/TypeTraits.h>
|
|
||||||
#include <vtkm/Types.h>
|
|
||||||
#include <vtkm/VectorAnalysis.h>
|
|
||||||
|
|
||||||
#include <vtkm/cont/DataSet.h>
|
|
||||||
#include <vtkm/cont/VirtualObjectHandle.h>
|
|
||||||
|
|
||||||
#include <vtkm/worklet/particleadvection/GridEvaluators.h>
|
|
||||||
#include <vtkm/worklet/particleadvection/IntegratorStatus.h>
|
|
||||||
#include <vtkm/worklet/particleadvection/Particles.h>
|
|
||||||
|
|
||||||
namespace vtkm
|
|
||||||
{
|
|
||||||
namespace worklet
|
|
||||||
{
|
|
||||||
namespace particleadvection
|
|
||||||
{
|
|
||||||
|
|
||||||
class IntegratorBase : public vtkm::cont::ExecutionObjectBase
|
|
||||||
{
|
|
||||||
protected:
|
|
||||||
VTKM_CONT
|
|
||||||
IntegratorBase() = default;
|
|
||||||
|
|
||||||
VTKM_CONT
|
|
||||||
IntegratorBase(vtkm::FloatDefault stepLength)
|
|
||||||
: StepLength(stepLength)
|
|
||||||
{
|
|
||||||
}
|
|
||||||
|
|
||||||
public:
|
|
||||||
class ExecObject : public vtkm::VirtualObjectBase
|
|
||||||
{
|
|
||||||
protected:
|
|
||||||
VTKM_EXEC_CONT
|
|
||||||
ExecObject(const vtkm::FloatDefault stepLength, vtkm::FloatDefault tolerance)
|
|
||||||
: StepLength(stepLength)
|
|
||||||
, Tolerance(tolerance)
|
|
||||||
{
|
|
||||||
}
|
|
||||||
|
|
||||||
public:
|
|
||||||
VTKM_EXEC
|
|
||||||
virtual IntegratorStatus Step(vtkm::Particle* inpos,
|
|
||||||
vtkm::FloatDefault& time,
|
|
||||||
vtkm::Vec3f& outpos) const = 0;
|
|
||||||
|
|
||||||
VTKM_EXEC
|
|
||||||
virtual IntegratorStatus SmallStep(vtkm::Particle* inpos,
|
|
||||||
vtkm::FloatDefault& time,
|
|
||||||
vtkm::Vec3f& outpos) const = 0;
|
|
||||||
|
|
||||||
protected:
|
|
||||||
vtkm::FloatDefault StepLength = 1.0f;
|
|
||||||
vtkm::FloatDefault Tolerance = 0.001f;
|
|
||||||
};
|
|
||||||
|
|
||||||
template <typename Device>
|
|
||||||
VTKM_CONT const ExecObject* PrepareForExecution(Device, vtkm::cont::Token& token) const
|
|
||||||
{
|
|
||||||
this->PrepareForExecutionImpl(
|
|
||||||
Device(),
|
|
||||||
const_cast<vtkm::cont::VirtualObjectHandle<ExecObject>&>(this->ExecObjectHandle),
|
|
||||||
token);
|
|
||||||
return this->ExecObjectHandle.PrepareForExecution(Device(), token);
|
|
||||||
}
|
|
||||||
|
|
||||||
private:
|
|
||||||
vtkm::cont::VirtualObjectHandle<ExecObject> ExecObjectHandle;
|
|
||||||
|
|
||||||
protected:
|
|
||||||
vtkm::FloatDefault StepLength;
|
|
||||||
vtkm::FloatDefault Tolerance =
|
|
||||||
std::numeric_limits<vtkm::FloatDefault>::epsilon() * static_cast<vtkm::FloatDefault>(100.0f);
|
|
||||||
|
|
||||||
VTKM_CONT virtual void PrepareForExecutionImpl(
|
|
||||||
vtkm::cont::DeviceAdapterId device,
|
|
||||||
vtkm::cont::VirtualObjectHandle<ExecObject>& execObjectHandle,
|
|
||||||
vtkm::cont::Token& token) const = 0;
|
|
||||||
|
|
||||||
template <typename FieldEvaluateType, typename DerivedType>
|
|
||||||
class ExecObjectBaseImpl : public ExecObject
|
|
||||||
{
|
|
||||||
protected:
|
|
||||||
VTKM_EXEC_CONT
|
|
||||||
ExecObjectBaseImpl(const FieldEvaluateType& evaluator,
|
|
||||||
vtkm::FloatDefault stepLength,
|
|
||||||
vtkm::FloatDefault tolerance)
|
|
||||||
: ExecObject(stepLength, tolerance)
|
|
||||||
, Evaluator(evaluator)
|
|
||||||
{
|
|
||||||
}
|
|
||||||
|
|
||||||
public:
|
|
||||||
VTKM_EXEC
|
|
||||||
IntegratorStatus Step(vtkm::Particle* particle,
|
|
||||||
vtkm::FloatDefault& time,
|
|
||||||
vtkm::Vec3f& outpos) const override
|
|
||||||
{
|
|
||||||
vtkm::Vec3f velocity(0, 0, 0);
|
|
||||||
auto status = this->CheckStep(particle, this->StepLength, velocity);
|
|
||||||
if (status.CheckOk())
|
|
||||||
{
|
|
||||||
outpos = particle->Pos + this->StepLength * velocity;
|
|
||||||
time += this->StepLength;
|
|
||||||
}
|
|
||||||
else
|
|
||||||
outpos = particle->Pos;
|
|
||||||
|
|
||||||
return status;
|
|
||||||
}
|
|
||||||
|
|
||||||
VTKM_EXEC
|
|
||||||
IntegratorStatus SmallStep(vtkm::Particle* particle,
|
|
||||||
vtkm::FloatDefault& time,
|
|
||||||
vtkm::Vec3f& outpos) const override
|
|
||||||
{
|
|
||||||
//Stepping by this->StepLength goes beyond the bounds of the dataset.
|
|
||||||
//We need to take an Euler step that goes outside of the dataset.
|
|
||||||
//Use a binary search to find the largest step INSIDE the dataset.
|
|
||||||
//Binary search uses a shrinking bracket of inside / outside, so when
|
|
||||||
//we terminate, the outside value is the stepsize that will nudge
|
|
||||||
//the particle outside the dataset.
|
|
||||||
|
|
||||||
//The binary search will be between {0, this->StepLength}
|
|
||||||
vtkm::FloatDefault stepRange[2] = { 0, this->StepLength };
|
|
||||||
|
|
||||||
vtkm::Vec3f currPos(particle->Pos);
|
|
||||||
vtkm::Vec3f currVelocity(0, 0, 0);
|
|
||||||
vtkm::VecVariable<vtkm::Vec3f, 2> currValue, tmp;
|
|
||||||
auto evalStatus = this->Evaluator.Evaluate(currPos, particle->Time, currValue);
|
|
||||||
if (evalStatus.CheckFail())
|
|
||||||
return IntegratorStatus(evalStatus);
|
|
||||||
|
|
||||||
const vtkm::FloatDefault eps = vtkm::Epsilon<vtkm::FloatDefault>() * 10;
|
|
||||||
vtkm::FloatDefault div = 1;
|
|
||||||
while ((stepRange[1] - stepRange[0]) > eps)
|
|
||||||
{
|
|
||||||
//Try a step midway between stepRange[0] and stepRange[1]
|
|
||||||
div *= 2;
|
|
||||||
vtkm::FloatDefault currStep = stepRange[0] + (this->StepLength / div);
|
|
||||||
|
|
||||||
//See if we can step by currStep
|
|
||||||
IntegratorStatus status = this->CheckStep(particle, currStep, currVelocity);
|
|
||||||
|
|
||||||
if (status.CheckOk()) //Integration step succedded.
|
|
||||||
{
|
|
||||||
//See if this point is in/out.
|
|
||||||
auto newPos = particle->Pos + currStep * currVelocity;
|
|
||||||
evalStatus = this->Evaluator.Evaluate(newPos, particle->Time + currStep, tmp);
|
|
||||||
if (evalStatus.CheckOk())
|
|
||||||
{
|
|
||||||
//Point still in. Update currPos and set range to {currStep, stepRange[1]}
|
|
||||||
currPos = newPos;
|
|
||||||
stepRange[0] = currStep;
|
|
||||||
}
|
|
||||||
else
|
|
||||||
{
|
|
||||||
//The step succedded, but the next point is outside.
|
|
||||||
//Step too long. Set range to: {stepRange[0], currStep} and continue.
|
|
||||||
stepRange[1] = currStep;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
else
|
|
||||||
{
|
|
||||||
//Step too long. Set range to: {stepRange[0], stepCurr} and continue.
|
|
||||||
stepRange[1] = currStep;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
evalStatus = this->Evaluator.Evaluate(currPos, particle->Time + stepRange[0], currValue);
|
|
||||||
//The eval at Time + stepRange[0] better be *inside*
|
|
||||||
VTKM_ASSERT(evalStatus.CheckOk() && !evalStatus.CheckSpatialBounds());
|
|
||||||
if (evalStatus.CheckFail() || evalStatus.CheckSpatialBounds())
|
|
||||||
return IntegratorStatus(evalStatus);
|
|
||||||
|
|
||||||
//Update the position and time.
|
|
||||||
outpos = currPos + stepRange[1] * particle->Velocity(currValue, stepRange[1]);
|
|
||||||
time += stepRange[1];
|
|
||||||
|
|
||||||
//Get the evaluation status for the point that is *just* outside of the data.
|
|
||||||
evalStatus = this->Evaluator.Evaluate(outpos, time, currValue);
|
|
||||||
|
|
||||||
//The eval should fail, and the point should be outside either spatially or temporally.
|
|
||||||
VTKM_ASSERT(evalStatus.CheckFail() &&
|
|
||||||
(evalStatus.CheckSpatialBounds() || evalStatus.CheckTemporalBounds()));
|
|
||||||
|
|
||||||
IntegratorStatus status(evalStatus);
|
|
||||||
status.SetOk(); //status is ok.
|
|
||||||
|
|
||||||
return status;
|
|
||||||
}
|
|
||||||
|
|
||||||
VTKM_EXEC
|
|
||||||
IntegratorStatus CheckStep(vtkm::Particle* particle,
|
|
||||||
vtkm::FloatDefault stepLength,
|
|
||||||
vtkm::Vec3f& velocity) const
|
|
||||||
{
|
|
||||||
return static_cast<const DerivedType*>(this)->CheckStep(particle, stepLength, velocity);
|
|
||||||
}
|
|
||||||
|
|
||||||
protected:
|
|
||||||
FieldEvaluateType Evaluator;
|
|
||||||
};
|
|
||||||
};
|
|
||||||
|
|
||||||
namespace detail
|
|
||||||
{
|
|
||||||
|
|
||||||
template <template <typename> class IntegratorType>
|
|
||||||
struct IntegratorPrepareForExecutionFunctor
|
|
||||||
{
|
|
||||||
template <typename Device, typename EvaluatorType>
|
|
||||||
VTKM_CONT bool operator()(
|
|
||||||
Device,
|
|
||||||
vtkm::cont::VirtualObjectHandle<IntegratorBase::ExecObject>& execObjectHandle,
|
|
||||||
const EvaluatorType& evaluator,
|
|
||||||
vtkm::FloatDefault stepLength,
|
|
||||||
vtkm::FloatDefault tolerance,
|
|
||||||
vtkm::cont::Token& token) const
|
|
||||||
{
|
|
||||||
IntegratorType<Device>* integrator = new IntegratorType<Device>(
|
|
||||||
evaluator.PrepareForExecution(Device(), token), stepLength, tolerance);
|
|
||||||
execObjectHandle.Reset(integrator);
|
|
||||||
return true;
|
|
||||||
}
|
|
||||||
};
|
|
||||||
|
|
||||||
} //namespace detail
|
|
||||||
} //namespace particleadvection
|
|
||||||
} //namespace worklet
|
|
||||||
} //namespace vtkm
|
|
||||||
|
|
||||||
#endif // vtk_m_worklet_particleadvection_IntegratorBase_h
|
|
@ -21,6 +21,8 @@
|
|||||||
#include <vtkm/Types.h>
|
#include <vtkm/Types.h>
|
||||||
#include <vtkm/VectorAnalysis.h>
|
#include <vtkm/VectorAnalysis.h>
|
||||||
|
|
||||||
|
#include <vtkm/worklet/particleadvection/GridEvaluatorStatus.h>
|
||||||
|
|
||||||
namespace vtkm
|
namespace vtkm
|
||||||
{
|
{
|
||||||
namespace worklet
|
namespace worklet
|
||||||
|
@ -22,8 +22,8 @@
|
|||||||
|
|
||||||
#include <vtkm/Particle.h>
|
#include <vtkm/Particle.h>
|
||||||
#include <vtkm/worklet/WorkletMapField.h>
|
#include <vtkm/worklet/WorkletMapField.h>
|
||||||
#include <vtkm/worklet/particleadvection/IntegratorBase.h>
|
|
||||||
#include <vtkm/worklet/particleadvection/Particles.h>
|
#include <vtkm/worklet/particleadvection/Particles.h>
|
||||||
|
#include <vtkm/worklet/particleadvection/Stepper.h>
|
||||||
|
|
||||||
#ifdef VTKM_CUDA
|
#ifdef VTKM_CUDA
|
||||||
#include <vtkm/cont/cuda/internal/ScopedCudaStackSize.h>
|
#include <vtkm/cont/cuda/internal/ScopedCudaStackSize.h>
|
||||||
@ -48,7 +48,7 @@ public:
|
|||||||
|
|
||||||
template <typename IntegratorType, typename IntegralCurveType>
|
template <typename IntegratorType, typename IntegralCurveType>
|
||||||
VTKM_EXEC void operator()(const vtkm::Id& idx,
|
VTKM_EXEC void operator()(const vtkm::Id& idx,
|
||||||
const IntegratorType* integrator,
|
const IntegratorType& integrator,
|
||||||
IntegralCurveType& integralCurve,
|
IntegralCurveType& integralCurve,
|
||||||
const vtkm::Id& maxSteps) const
|
const vtkm::Id& maxSteps) const
|
||||||
{
|
{
|
||||||
@ -66,7 +66,7 @@ public:
|
|||||||
do
|
do
|
||||||
{
|
{
|
||||||
vtkm::Vec3f outpos;
|
vtkm::Vec3f outpos;
|
||||||
auto status = integrator->Step(&particle, time, outpos);
|
auto status = integrator.Step(particle, time, outpos);
|
||||||
if (status.CheckOk())
|
if (status.CheckOk())
|
||||||
{
|
{
|
||||||
integralCurve.StepUpdate(idx, time, outpos);
|
integralCurve.StepUpdate(idx, time, outpos);
|
||||||
@ -78,7 +78,7 @@ public:
|
|||||||
//Try and take a step just past the boundary.
|
//Try and take a step just past the boundary.
|
||||||
else if (status.CheckSpatialBounds())
|
else if (status.CheckSpatialBounds())
|
||||||
{
|
{
|
||||||
status = integrator->SmallStep(&particle, time, outpos);
|
status = integrator.SmallStep(particle, time, outpos);
|
||||||
if (status.CheckOk())
|
if (status.CheckOk())
|
||||||
{
|
{
|
||||||
integralCurve.StepUpdate(idx, time, outpos);
|
integralCurve.StepUpdate(idx, time, outpos);
|
||||||
|
@ -13,8 +13,6 @@
|
|||||||
#ifndef vtk_m_worklet_particleadvection_RK4Integrator_h
|
#ifndef vtk_m_worklet_particleadvection_RK4Integrator_h
|
||||||
#define vtk_m_worklet_particleadvection_RK4Integrator_h
|
#define vtk_m_worklet_particleadvection_RK4Integrator_h
|
||||||
|
|
||||||
#include <vtkm/worklet/particleadvection/IntegratorBase.h>
|
|
||||||
|
|
||||||
namespace vtkm
|
namespace vtkm
|
||||||
{
|
{
|
||||||
namespace worklet
|
namespace worklet
|
||||||
@ -22,116 +20,98 @@ namespace worklet
|
|||||||
namespace particleadvection
|
namespace particleadvection
|
||||||
{
|
{
|
||||||
|
|
||||||
template <typename FieldEvaluateType>
|
template <typename ExecEvaluatorType>
|
||||||
class RK4Integrator : public IntegratorBase
|
class ExecRK4Integrator
|
||||||
{
|
{
|
||||||
|
public:
|
||||||
|
VTKM_EXEC_CONT
|
||||||
|
ExecRK4Integrator(const ExecEvaluatorType& evaluator)
|
||||||
|
: Evaluator(evaluator)
|
||||||
|
{
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename Particle>
|
||||||
|
VTKM_EXEC IntegratorStatus CheckStep(Particle& particle,
|
||||||
|
vtkm::FloatDefault stepLength,
|
||||||
|
vtkm::Vec3f& velocity) const
|
||||||
|
{
|
||||||
|
auto time = particle.Time;
|
||||||
|
auto inpos = particle.Pos;
|
||||||
|
vtkm::FloatDefault boundary = this->Evaluator.GetTemporalBoundary(static_cast<vtkm::Id>(1));
|
||||||
|
if ((time + stepLength + vtkm::Epsilon<vtkm::FloatDefault>() - boundary) > 0.0)
|
||||||
|
stepLength = boundary - time;
|
||||||
|
|
||||||
|
//k1 = F(p,t)
|
||||||
|
//k2 = F(p+hk1/2, t+h/2
|
||||||
|
//k3 = F(p+hk2/2, t+h/2
|
||||||
|
//k4 = F(p+hk3, t+h)
|
||||||
|
//Yn+1 = Yn + 1/6 h (k1+2k2+2k3+k4)
|
||||||
|
|
||||||
|
vtkm::FloatDefault var1 = (stepLength / static_cast<vtkm::FloatDefault>(2));
|
||||||
|
vtkm::FloatDefault var2 = time + var1;
|
||||||
|
vtkm::FloatDefault var3 = time + stepLength;
|
||||||
|
|
||||||
|
vtkm::Vec3f v1 = vtkm::TypeTraits<vtkm::Vec3f>::ZeroInitialization();
|
||||||
|
vtkm::Vec3f v2 = v1, v3 = v1, v4 = v1;
|
||||||
|
vtkm::VecVariable<vtkm::Vec3f, 2> k1, k2, k3, k4;
|
||||||
|
|
||||||
|
GridEvaluatorStatus evalStatus;
|
||||||
|
|
||||||
|
evalStatus = this->Evaluator.Evaluate(inpos, time, k1);
|
||||||
|
if (evalStatus.CheckFail())
|
||||||
|
return IntegratorStatus(evalStatus);
|
||||||
|
v1 = particle.Velocity(k1, stepLength);
|
||||||
|
|
||||||
|
evalStatus = this->Evaluator.Evaluate(inpos + var1 * v1, var2, k2);
|
||||||
|
if (evalStatus.CheckFail())
|
||||||
|
return IntegratorStatus(evalStatus);
|
||||||
|
v2 = particle.Velocity(k2, stepLength);
|
||||||
|
|
||||||
|
evalStatus = this->Evaluator.Evaluate(inpos + var1 * v2, var2, k3);
|
||||||
|
if (evalStatus.CheckFail())
|
||||||
|
return IntegratorStatus(evalStatus);
|
||||||
|
v3 = particle.Velocity(k3, stepLength);
|
||||||
|
|
||||||
|
evalStatus = this->Evaluator.Evaluate(inpos + stepLength * v3, var3, k4);
|
||||||
|
if (evalStatus.CheckFail())
|
||||||
|
return IntegratorStatus(evalStatus);
|
||||||
|
v4 = particle.Velocity(k4, stepLength);
|
||||||
|
|
||||||
|
velocity = (v1 + 2 * v2 + 2 * v3 + v4) / static_cast<vtkm::FloatDefault>(6);
|
||||||
|
|
||||||
|
return IntegratorStatus(evalStatus);
|
||||||
|
}
|
||||||
|
|
||||||
|
private:
|
||||||
|
ExecEvaluatorType Evaluator;
|
||||||
|
};
|
||||||
|
|
||||||
|
template <typename EvaluatorType>
|
||||||
|
class RK4Integrator
|
||||||
|
{
|
||||||
|
private:
|
||||||
|
EvaluatorType Evaluator;
|
||||||
|
|
||||||
public:
|
public:
|
||||||
VTKM_CONT
|
VTKM_CONT
|
||||||
RK4Integrator() = default;
|
RK4Integrator() = default;
|
||||||
|
|
||||||
VTKM_CONT
|
VTKM_CONT
|
||||||
RK4Integrator(const FieldEvaluateType& evaluator, vtkm::FloatDefault stepLength)
|
RK4Integrator(const EvaluatorType& evaluator)
|
||||||
: IntegratorBase(stepLength)
|
: Evaluator(evaluator)
|
||||||
, Evaluator(evaluator)
|
|
||||||
{
|
{
|
||||||
}
|
}
|
||||||
|
|
||||||
template <typename Device>
|
VTKM_CONT auto PrepareForExecution(vtkm::cont::DeviceAdapterId device,
|
||||||
class ExecObject
|
vtkm::cont::Token& token) const
|
||||||
: public IntegratorBase::ExecObjectBaseImpl<
|
-> ExecRK4Integrator<decltype(this->Evaluator.PrepareForExecution(device, token))>
|
||||||
vtkm::cont::internal::ExecutionObjectType<FieldEvaluateType, Device>,
|
|
||||||
typename RK4Integrator::template ExecObject<Device>>
|
|
||||||
{
|
{
|
||||||
VTKM_IS_DEVICE_ADAPTER_TAG(Device);
|
auto evaluator = this->Evaluator.PrepareForExecution(device, token);
|
||||||
|
using ExecEvaluatorType = decltype(evaluator);
|
||||||
using FieldEvaluateExecType =
|
return ExecRK4Integrator<ExecEvaluatorType>(evaluator);
|
||||||
vtkm::cont::internal::ExecutionObjectType<FieldEvaluateType, Device>;
|
|
||||||
using Superclass =
|
|
||||||
IntegratorBase::ExecObjectBaseImpl<FieldEvaluateExecType,
|
|
||||||
typename RK4Integrator::template ExecObject<Device>>;
|
|
||||||
|
|
||||||
public:
|
|
||||||
VTKM_EXEC_CONT
|
|
||||||
ExecObject(const FieldEvaluateExecType& evaluator,
|
|
||||||
vtkm::FloatDefault stepLength,
|
|
||||||
vtkm::FloatDefault tolerance)
|
|
||||||
: Superclass(evaluator, stepLength, tolerance)
|
|
||||||
{
|
|
||||||
}
|
|
||||||
|
|
||||||
VTKM_EXEC
|
|
||||||
IntegratorStatus CheckStep(vtkm::Particle* particle,
|
|
||||||
vtkm::FloatDefault stepLength,
|
|
||||||
vtkm::Vec3f& velocity) const
|
|
||||||
{
|
|
||||||
auto time = particle->Time;
|
|
||||||
auto inpos = particle->Pos;
|
|
||||||
vtkm::FloatDefault boundary = this->Evaluator.GetTemporalBoundary(static_cast<vtkm::Id>(1));
|
|
||||||
if ((time + stepLength + vtkm::Epsilon<vtkm::FloatDefault>() - boundary) > 0.0)
|
|
||||||
stepLength = boundary - time;
|
|
||||||
|
|
||||||
//k1 = F(p,t)
|
|
||||||
//k2 = F(p+hk1/2, t+h/2
|
|
||||||
//k3 = F(p+hk2/2, t+h/2
|
|
||||||
//k4 = F(p+hk3, t+h)
|
|
||||||
//Yn+1 = Yn + 1/6 h (k1+2k2+2k3+k4)
|
|
||||||
|
|
||||||
vtkm::FloatDefault var1 = (stepLength / static_cast<vtkm::FloatDefault>(2));
|
|
||||||
vtkm::FloatDefault var2 = time + var1;
|
|
||||||
vtkm::FloatDefault var3 = time + stepLength;
|
|
||||||
|
|
||||||
vtkm::Vec3f v1 = vtkm::TypeTraits<vtkm::Vec3f>::ZeroInitialization();
|
|
||||||
vtkm::Vec3f v2 = v1, v3 = v1, v4 = v1;
|
|
||||||
vtkm::VecVariable<vtkm::Vec3f, 2> k1, k2, k3, k4;
|
|
||||||
|
|
||||||
GridEvaluatorStatus evalStatus;
|
|
||||||
|
|
||||||
evalStatus = this->Evaluator.Evaluate(inpos, time, k1);
|
|
||||||
if (evalStatus.CheckFail())
|
|
||||||
return IntegratorStatus(evalStatus);
|
|
||||||
v1 = particle->Velocity(k1, stepLength);
|
|
||||||
|
|
||||||
evalStatus = this->Evaluator.Evaluate(inpos + var1 * v1, var2, k2);
|
|
||||||
if (evalStatus.CheckFail())
|
|
||||||
return IntegratorStatus(evalStatus);
|
|
||||||
v2 = particle->Velocity(k2, stepLength);
|
|
||||||
|
|
||||||
evalStatus = this->Evaluator.Evaluate(inpos + var1 * v2, var2, k3);
|
|
||||||
if (evalStatus.CheckFail())
|
|
||||||
return IntegratorStatus(evalStatus);
|
|
||||||
v3 = particle->Velocity(k3, stepLength);
|
|
||||||
|
|
||||||
evalStatus = this->Evaluator.Evaluate(inpos + stepLength * v3, var3, k4);
|
|
||||||
if (evalStatus.CheckFail())
|
|
||||||
return IntegratorStatus(evalStatus);
|
|
||||||
v4 = particle->Velocity(k4, stepLength);
|
|
||||||
|
|
||||||
velocity = (v1 + 2 * v2 + 2 * v3 + v4) / static_cast<vtkm::FloatDefault>(6);
|
|
||||||
|
|
||||||
return IntegratorStatus(evalStatus);
|
|
||||||
}
|
|
||||||
};
|
|
||||||
|
|
||||||
private:
|
|
||||||
FieldEvaluateType Evaluator;
|
|
||||||
|
|
||||||
protected:
|
|
||||||
VTKM_CONT virtual void PrepareForExecutionImpl(
|
|
||||||
vtkm::cont::DeviceAdapterId device,
|
|
||||||
vtkm::cont::VirtualObjectHandle<IntegratorBase::ExecObject>& execObjectHandle,
|
|
||||||
vtkm::cont::Token& token) const override
|
|
||||||
{
|
|
||||||
vtkm::cont::TryExecuteOnDevice(device,
|
|
||||||
detail::IntegratorPrepareForExecutionFunctor<ExecObject>(),
|
|
||||||
execObjectHandle,
|
|
||||||
this->Evaluator,
|
|
||||||
this->StepLength,
|
|
||||||
this->Tolerance,
|
|
||||||
token);
|
|
||||||
}
|
}
|
||||||
};
|
};
|
||||||
|
|
||||||
|
|
||||||
} //namespace particleadvection
|
} //namespace particleadvection
|
||||||
} //namespace worklet
|
} //namespace worklet
|
||||||
} //namespace vtkm
|
} //namespace vtkm
|
||||||
|
206
vtkm/worklet/particleadvection/Stepper.h
Normal file
206
vtkm/worklet/particleadvection/Stepper.h
Normal file
@ -0,0 +1,206 @@
|
|||||||
|
//=============================================================================
|
||||||
|
//
|
||||||
|
// 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_particleadvection_Stepper_h
|
||||||
|
#define vtk_m_worklet_particleadvection_Stepper_h
|
||||||
|
|
||||||
|
#include <limits>
|
||||||
|
|
||||||
|
#include <vtkm/Bitset.h>
|
||||||
|
#include <vtkm/TypeTraits.h>
|
||||||
|
#include <vtkm/Types.h>
|
||||||
|
#include <vtkm/VectorAnalysis.h>
|
||||||
|
|
||||||
|
#include <vtkm/cont/DataSet.h>
|
||||||
|
#include <vtkm/cont/VirtualObjectHandle.h>
|
||||||
|
|
||||||
|
#include <vtkm/worklet/particleadvection/GridEvaluators.h>
|
||||||
|
#include <vtkm/worklet/particleadvection/IntegratorStatus.h>
|
||||||
|
#include <vtkm/worklet/particleadvection/Particles.h>
|
||||||
|
|
||||||
|
namespace vtkm
|
||||||
|
{
|
||||||
|
namespace worklet
|
||||||
|
{
|
||||||
|
namespace particleadvection
|
||||||
|
{
|
||||||
|
|
||||||
|
template <typename ExecIntegratorType, typename ExecEvaluatorType>
|
||||||
|
class StepperImpl
|
||||||
|
{
|
||||||
|
private:
|
||||||
|
ExecIntegratorType Integrator;
|
||||||
|
ExecEvaluatorType Evaluator;
|
||||||
|
vtkm::FloatDefault DeltaT;
|
||||||
|
vtkm::FloatDefault Tolerance;
|
||||||
|
|
||||||
|
public:
|
||||||
|
VTKM_EXEC_CONT
|
||||||
|
StepperImpl(const ExecIntegratorType& integrator,
|
||||||
|
const ExecEvaluatorType& evaluator,
|
||||||
|
const vtkm::FloatDefault deltaT,
|
||||||
|
const vtkm::FloatDefault tolerance)
|
||||||
|
: Integrator(integrator)
|
||||||
|
, Evaluator(evaluator)
|
||||||
|
, DeltaT(deltaT)
|
||||||
|
, Tolerance(tolerance)
|
||||||
|
{
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename Particle>
|
||||||
|
VTKM_EXEC IntegratorStatus Step(Particle& particle,
|
||||||
|
vtkm::FloatDefault& time,
|
||||||
|
vtkm::Vec3f& outpos) const
|
||||||
|
{
|
||||||
|
vtkm::Vec3f velocity(0, 0, 0);
|
||||||
|
auto status = this->Integrator.CheckStep(particle, this->DeltaT, velocity);
|
||||||
|
if (status.CheckOk())
|
||||||
|
{
|
||||||
|
outpos = particle.Pos + this->DeltaT * velocity;
|
||||||
|
time += this->DeltaT;
|
||||||
|
}
|
||||||
|
else
|
||||||
|
outpos = particle.Pos;
|
||||||
|
|
||||||
|
return status;
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename Particle>
|
||||||
|
VTKM_EXEC IntegratorStatus SmallStep(Particle particle,
|
||||||
|
vtkm::FloatDefault& time,
|
||||||
|
vtkm::Vec3f& outpos) const
|
||||||
|
{
|
||||||
|
//Stepping by this->DeltaT goes beyond the bounds of the dataset.
|
||||||
|
//We need to take an Euler step that goes outside of the dataset.
|
||||||
|
//Use a binary search to find the largest step INSIDE the dataset.
|
||||||
|
//Binary search uses a shrinking bracket of inside / outside, so when
|
||||||
|
//we terminate, the outside value is the stepsize that will nudge
|
||||||
|
//the particle outside the dataset.
|
||||||
|
|
||||||
|
//The binary search will be between {0, this->DeltaT}
|
||||||
|
vtkm::FloatDefault stepRange[2] = { 0, this->DeltaT };
|
||||||
|
|
||||||
|
vtkm::Vec3f currPos(particle.Pos);
|
||||||
|
vtkm::Vec3f currVelocity(0, 0, 0);
|
||||||
|
vtkm::VecVariable<vtkm::Vec3f, 2> currValue, tmp;
|
||||||
|
auto evalStatus = this->Evaluator.Evaluate(currPos, particle.Time, currValue);
|
||||||
|
if (evalStatus.CheckFail())
|
||||||
|
return IntegratorStatus(evalStatus);
|
||||||
|
|
||||||
|
const vtkm::FloatDefault eps = vtkm::Epsilon<vtkm::FloatDefault>() * 10;
|
||||||
|
vtkm::FloatDefault div = 1;
|
||||||
|
while ((stepRange[1] - stepRange[0]) > eps)
|
||||||
|
{
|
||||||
|
//Try a step midway between stepRange[0] and stepRange[1]
|
||||||
|
div *= 2;
|
||||||
|
vtkm::FloatDefault currStep = stepRange[0] + (this->DeltaT / div);
|
||||||
|
|
||||||
|
//See if we can step by currStep
|
||||||
|
IntegratorStatus status = this->Integrator.CheckStep(particle, currStep, currVelocity);
|
||||||
|
|
||||||
|
if (status.CheckOk()) //Integration step succedded.
|
||||||
|
{
|
||||||
|
//See if this point is in/out.
|
||||||
|
auto newPos = particle.Pos + currStep * currVelocity;
|
||||||
|
evalStatus = this->Evaluator.Evaluate(newPos, particle.Time + currStep, tmp);
|
||||||
|
if (evalStatus.CheckOk())
|
||||||
|
{
|
||||||
|
//Point still in. Update currPos and set range to {currStep, stepRange[1]}
|
||||||
|
currPos = newPos;
|
||||||
|
stepRange[0] = currStep;
|
||||||
|
}
|
||||||
|
else
|
||||||
|
{
|
||||||
|
//The step succedded, but the next point is outside.
|
||||||
|
//Step too long. Set range to: {stepRange[0], currStep} and continue.
|
||||||
|
stepRange[1] = currStep;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
else
|
||||||
|
{
|
||||||
|
//Step too long. Set range to: {stepRange[0], stepCurr} and continue.
|
||||||
|
stepRange[1] = currStep;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
evalStatus = this->Evaluator.Evaluate(currPos, particle.Time + stepRange[0], currValue);
|
||||||
|
// The eval at Time + stepRange[0] better be *inside*
|
||||||
|
VTKM_ASSERT(evalStatus.CheckOk() && !evalStatus.CheckSpatialBounds());
|
||||||
|
if (evalStatus.CheckFail() || evalStatus.CheckSpatialBounds())
|
||||||
|
return IntegratorStatus(evalStatus);
|
||||||
|
|
||||||
|
// Update the position and time.
|
||||||
|
outpos = currPos + stepRange[1] * particle.Velocity(currValue, stepRange[1]);
|
||||||
|
time += stepRange[1];
|
||||||
|
|
||||||
|
// Get the evaluation status for the point that is *just* outside of the data.
|
||||||
|
evalStatus = this->Evaluator.Evaluate(outpos, time, currValue);
|
||||||
|
|
||||||
|
// The eval should fail, and the point should be outside either spatially or temporally.
|
||||||
|
VTKM_ASSERT(evalStatus.CheckFail() &&
|
||||||
|
(evalStatus.CheckSpatialBounds() || evalStatus.CheckTemporalBounds()));
|
||||||
|
|
||||||
|
IntegratorStatus status(evalStatus);
|
||||||
|
status.SetOk(); //status is ok.
|
||||||
|
|
||||||
|
return status;
|
||||||
|
}
|
||||||
|
};
|
||||||
|
|
||||||
|
|
||||||
|
template <typename IntegratorType, typename EvaluatorType>
|
||||||
|
class Stepper : public vtkm::cont::ExecutionObjectBase
|
||||||
|
{
|
||||||
|
private:
|
||||||
|
IntegratorType Integrator;
|
||||||
|
EvaluatorType Evaluator;
|
||||||
|
vtkm::FloatDefault DeltaT;
|
||||||
|
vtkm::FloatDefault Tolerance =
|
||||||
|
std::numeric_limits<vtkm::FloatDefault>::epsilon() * static_cast<vtkm::FloatDefault>(100.0f);
|
||||||
|
|
||||||
|
public:
|
||||||
|
VTKM_CONT
|
||||||
|
Stepper() = default;
|
||||||
|
|
||||||
|
VTKM_CONT
|
||||||
|
Stepper(const EvaluatorType& evaluator, const vtkm::FloatDefault deltaT)
|
||||||
|
: Integrator(IntegratorType(evaluator))
|
||||||
|
, Evaluator(evaluator)
|
||||||
|
, DeltaT(deltaT)
|
||||||
|
{
|
||||||
|
}
|
||||||
|
|
||||||
|
VTKM_CONT
|
||||||
|
void SetTolerance(vtkm::FloatDefault tolerance) { this->Tolerance = tolerance; }
|
||||||
|
|
||||||
|
public:
|
||||||
|
/// Return the StepperImpl object
|
||||||
|
/// Prepares the execution object of Stepper
|
||||||
|
VTKM_CONT auto PrepareForExecution(vtkm::cont::DeviceAdapterId device,
|
||||||
|
vtkm::cont::Token& token) const
|
||||||
|
-> StepperImpl<decltype(this->Integrator.PrepareForExecution(device, token)),
|
||||||
|
decltype(this->Evaluator.PrepareForExecution(device, token))>
|
||||||
|
{
|
||||||
|
auto integrator = this->Integrator.PrepareForExecution(device, token);
|
||||||
|
auto evaluator = this->Evaluator.PrepareForExecution(device, token);
|
||||||
|
using ExecIntegratorType = decltype(integrator);
|
||||||
|
using ExecEvaluatorType = decltype(evaluator);
|
||||||
|
return StepperImpl<ExecIntegratorType, ExecEvaluatorType>(
|
||||||
|
integrator, evaluator, this->DeltaT, this->Tolerance);
|
||||||
|
}
|
||||||
|
};
|
||||||
|
|
||||||
|
} //namespace particleadvection
|
||||||
|
} //namespace worklet
|
||||||
|
} //namespace vtkm
|
||||||
|
|
||||||
|
#endif // vtk_m_worklet_particleadvection_Stepper_h
|
@ -19,9 +19,9 @@
|
|||||||
#include <vtkm/worklet/particleadvection/EulerIntegrator.h>
|
#include <vtkm/worklet/particleadvection/EulerIntegrator.h>
|
||||||
#include <vtkm/worklet/particleadvection/Field.h>
|
#include <vtkm/worklet/particleadvection/Field.h>
|
||||||
#include <vtkm/worklet/particleadvection/GridEvaluators.h>
|
#include <vtkm/worklet/particleadvection/GridEvaluators.h>
|
||||||
#include <vtkm/worklet/particleadvection/IntegratorBase.h>
|
|
||||||
#include <vtkm/worklet/particleadvection/Particles.h>
|
#include <vtkm/worklet/particleadvection/Particles.h>
|
||||||
#include <vtkm/worklet/particleadvection/RK4Integrator.h>
|
#include <vtkm/worklet/particleadvection/RK4Integrator.h>
|
||||||
|
#include <vtkm/worklet/particleadvection/Stepper.h>
|
||||||
#include <vtkm/worklet/testing/GenerateTestDataSets.h>
|
#include <vtkm/worklet/testing/GenerateTestDataSets.h>
|
||||||
|
|
||||||
#include <random>
|
#include <random>
|
||||||
@ -189,16 +189,16 @@ public:
|
|||||||
|
|
||||||
using ExecutionSignature = void(_1, _2, _3, _4);
|
using ExecutionSignature = void(_1, _2, _3, _4);
|
||||||
|
|
||||||
template <typename IntegratorType>
|
template <typename Particle, typename IntegratorType>
|
||||||
VTKM_EXEC void operator()(vtkm::Particle& pointIn,
|
VTKM_EXEC void operator()(Particle& pointIn,
|
||||||
const IntegratorType* integrator,
|
const IntegratorType integrator,
|
||||||
vtkm::worklet::particleadvection::IntegratorStatus& status,
|
vtkm::worklet::particleadvection::IntegratorStatus& status,
|
||||||
vtkm::Vec3f& pointOut) const
|
vtkm::Vec3f& pointOut) const
|
||||||
{
|
{
|
||||||
vtkm::FloatDefault time = 0;
|
vtkm::FloatDefault time = 0;
|
||||||
status = integrator->Step(&pointIn, time, pointOut);
|
status = integrator.Step(pointIn, time, pointOut);
|
||||||
if (status.CheckSpatialBounds())
|
if (status.CheckSpatialBounds())
|
||||||
status = integrator->SmallStep(&pointIn, time, pointOut);
|
status = integrator.SmallStep(pointIn, time, pointOut);
|
||||||
}
|
}
|
||||||
};
|
};
|
||||||
|
|
||||||
@ -272,6 +272,7 @@ void TestEvaluators()
|
|||||||
using FieldType = vtkm::worklet::particleadvection::VelocityField<FieldHandle>;
|
using FieldType = vtkm::worklet::particleadvection::VelocityField<FieldHandle>;
|
||||||
using GridEvalType = vtkm::worklet::particleadvection::GridEvaluator<FieldType>;
|
using GridEvalType = vtkm::worklet::particleadvection::GridEvaluator<FieldType>;
|
||||||
using RK4Type = vtkm::worklet::particleadvection::RK4Integrator<GridEvalType>;
|
using RK4Type = vtkm::worklet::particleadvection::RK4Integrator<GridEvalType>;
|
||||||
|
using Stepper = vtkm::worklet::particleadvection::Stepper<RK4Type, GridEvalType>;
|
||||||
|
|
||||||
std::vector<vtkm::Vec3f> vecs;
|
std::vector<vtkm::Vec3f> vecs;
|
||||||
vtkm::FloatDefault vals[3] = { -1., 0., 1. };
|
vtkm::FloatDefault vals[3] = { -1., 0., 1. };
|
||||||
@ -351,7 +352,7 @@ void TestEvaluators()
|
|||||||
GridEvalType gridEval(ds.GetCoordinateSystem(), ds.GetCellSet(), velocities);
|
GridEvalType gridEval(ds.GetCoordinateSystem(), ds.GetCellSet(), velocities);
|
||||||
ValidateEvaluator(gridEval, pointIns, vec, "grid evaluator");
|
ValidateEvaluator(gridEval, pointIns, vec, "grid evaluator");
|
||||||
|
|
||||||
RK4Type rk4(gridEval, stepSize);
|
Stepper rk4(gridEval, stepSize);
|
||||||
ValidateIntegrator(rk4, pointIns, stepResult, "constant vector RK4");
|
ValidateIntegrator(rk4, pointIns, stepResult, "constant vector RK4");
|
||||||
ValidateIntegratorForBoundary(bound, rk4, boundaryPoints, "constant vector RK4");
|
ValidateIntegratorForBoundary(bound, rk4, boundaryPoints, "constant vector RK4");
|
||||||
}
|
}
|
||||||
@ -366,6 +367,7 @@ void TestGhostCellEvaluators()
|
|||||||
using FieldType = vtkm::worklet::particleadvection::VelocityField<FieldHandle>;
|
using FieldType = vtkm::worklet::particleadvection::VelocityField<FieldHandle>;
|
||||||
using GridEvalType = vtkm::worklet::particleadvection::GridEvaluator<FieldType>;
|
using GridEvalType = vtkm::worklet::particleadvection::GridEvaluator<FieldType>;
|
||||||
using RK4Type = vtkm::worklet::particleadvection::RK4Integrator<GridEvalType>;
|
using RK4Type = vtkm::worklet::particleadvection::RK4Integrator<GridEvalType>;
|
||||||
|
using Stepper = vtkm::worklet::particleadvection::Stepper<RK4Type, GridEvalType>;
|
||||||
|
|
||||||
constexpr vtkm::Id nX = 6;
|
constexpr vtkm::Id nX = 6;
|
||||||
constexpr vtkm::Id nY = 6;
|
constexpr vtkm::Id nY = 6;
|
||||||
@ -391,7 +393,7 @@ void TestGhostCellEvaluators()
|
|||||||
GridEvalType gridEval(ds, velocities);
|
GridEvalType gridEval(ds, velocities);
|
||||||
|
|
||||||
vtkm::FloatDefault stepSize = static_cast<vtkm::FloatDefault>(0.1);
|
vtkm::FloatDefault stepSize = static_cast<vtkm::FloatDefault>(0.1);
|
||||||
RK4Type rk4(gridEval, stepSize);
|
Stepper rk4(gridEval, stepSize);
|
||||||
|
|
||||||
vtkm::worklet::ParticleAdvection pa;
|
vtkm::worklet::ParticleAdvection pa;
|
||||||
std::vector<vtkm::Particle> seeds;
|
std::vector<vtkm::Particle> seeds;
|
||||||
@ -494,14 +496,16 @@ void TestIntegrators()
|
|||||||
{
|
{
|
||||||
auto seeds = vtkm::cont::make_ArrayHandle(points, vtkm::CopyFlag::On);
|
auto seeds = vtkm::cont::make_ArrayHandle(points, vtkm::CopyFlag::On);
|
||||||
using IntegratorType = vtkm::worklet::particleadvection::RK4Integrator<GridEvalType>;
|
using IntegratorType = vtkm::worklet::particleadvection::RK4Integrator<GridEvalType>;
|
||||||
IntegratorType rk4(eval, stepSize);
|
using Stepper = vtkm::worklet::particleadvection::Stepper<IntegratorType, GridEvalType>;
|
||||||
|
Stepper rk4(eval, stepSize);
|
||||||
res = pa.Run(rk4, seeds, maxSteps);
|
res = pa.Run(rk4, seeds, maxSteps);
|
||||||
ValidateParticleAdvectionResult(res, nSeeds, maxSteps);
|
ValidateParticleAdvectionResult(res, nSeeds, maxSteps);
|
||||||
}
|
}
|
||||||
{
|
{
|
||||||
auto seeds = vtkm::cont::make_ArrayHandle(points, vtkm::CopyFlag::On);
|
auto seeds = vtkm::cont::make_ArrayHandle(points, vtkm::CopyFlag::On);
|
||||||
using IntegratorType = vtkm::worklet::particleadvection::EulerIntegrator<GridEvalType>;
|
using IntegratorType = vtkm::worklet::particleadvection::EulerIntegrator<GridEvalType>;
|
||||||
IntegratorType euler(eval, stepSize);
|
using Stepper = vtkm::worklet::particleadvection::Stepper<IntegratorType, GridEvalType>;
|
||||||
|
Stepper euler(eval, stepSize);
|
||||||
res = pa.Run(euler, seeds, maxSteps);
|
res = pa.Run(euler, seeds, maxSteps);
|
||||||
ValidateParticleAdvectionResult(res, nSeeds, maxSteps);
|
ValidateParticleAdvectionResult(res, nSeeds, maxSteps);
|
||||||
}
|
}
|
||||||
@ -514,6 +518,7 @@ void TestParticleWorkletsWithDataSetTypes()
|
|||||||
using FieldType = vtkm::worklet::particleadvection::VelocityField<FieldHandle>;
|
using FieldType = vtkm::worklet::particleadvection::VelocityField<FieldHandle>;
|
||||||
using GridEvalType = vtkm::worklet::particleadvection::GridEvaluator<FieldType>;
|
using GridEvalType = vtkm::worklet::particleadvection::GridEvaluator<FieldType>;
|
||||||
using RK4Type = vtkm::worklet::particleadvection::RK4Integrator<GridEvalType>;
|
using RK4Type = vtkm::worklet::particleadvection::RK4Integrator<GridEvalType>;
|
||||||
|
using Stepper = vtkm::worklet::particleadvection::Stepper<RK4Type, GridEvalType>;
|
||||||
vtkm::FloatDefault stepSize = 0.01f;
|
vtkm::FloatDefault stepSize = 0.01f;
|
||||||
|
|
||||||
const vtkm::Id3 dims(5, 5, 5);
|
const vtkm::Id3 dims(5, 5, 5);
|
||||||
@ -555,7 +560,7 @@ void TestParticleWorkletsWithDataSetTypes()
|
|||||||
for (auto& ds : dataSets)
|
for (auto& ds : dataSets)
|
||||||
{
|
{
|
||||||
GridEvalType eval(ds.GetCoordinateSystem(), ds.GetCellSet(), velocities);
|
GridEvalType eval(ds.GetCoordinateSystem(), ds.GetCellSet(), velocities);
|
||||||
RK4Type rk4(eval, stepSize);
|
Stepper rk4(eval, stepSize);
|
||||||
|
|
||||||
//Do 4 tests on each dataset.
|
//Do 4 tests on each dataset.
|
||||||
//Particle advection worklet with and without steps taken.
|
//Particle advection worklet with and without steps taken.
|
||||||
@ -616,13 +621,15 @@ void TestParticleStatus()
|
|||||||
using FieldType = vtkm::worklet::particleadvection::VelocityField<FieldHandle>;
|
using FieldType = vtkm::worklet::particleadvection::VelocityField<FieldHandle>;
|
||||||
using GridEvalType = vtkm::worklet::particleadvection::GridEvaluator<FieldType>;
|
using GridEvalType = vtkm::worklet::particleadvection::GridEvaluator<FieldType>;
|
||||||
using RK4Type = vtkm::worklet::particleadvection::RK4Integrator<GridEvalType>;
|
using RK4Type = vtkm::worklet::particleadvection::RK4Integrator<GridEvalType>;
|
||||||
|
using Stepper = vtkm::worklet::particleadvection::Stepper<RK4Type, GridEvalType>;
|
||||||
|
|
||||||
vtkm::Id maxSteps = 1000;
|
vtkm::Id maxSteps = 1000;
|
||||||
vtkm::FloatDefault stepSize = 0.01f;
|
vtkm::FloatDefault stepSize = 0.01f;
|
||||||
|
|
||||||
FieldType velocities(fieldArray);
|
FieldType velocities(fieldArray);
|
||||||
|
|
||||||
GridEvalType eval(ds, velocities);
|
GridEvalType eval(ds, velocities);
|
||||||
RK4Type rk4(eval, stepSize);
|
Stepper rk4(eval, stepSize);
|
||||||
|
|
||||||
vtkm::worklet::ParticleAdvection pa;
|
vtkm::worklet::ParticleAdvection pa;
|
||||||
std::vector<vtkm::Particle> pts;
|
std::vector<vtkm::Particle> pts;
|
||||||
@ -645,6 +652,7 @@ void TestWorkletsBasic()
|
|||||||
using FieldType = vtkm::worklet::particleadvection::VelocityField<FieldHandle>;
|
using FieldType = vtkm::worklet::particleadvection::VelocityField<FieldHandle>;
|
||||||
using GridEvalType = vtkm::worklet::particleadvection::GridEvaluator<FieldType>;
|
using GridEvalType = vtkm::worklet::particleadvection::GridEvaluator<FieldType>;
|
||||||
using RK4Type = vtkm::worklet::particleadvection::RK4Integrator<GridEvalType>;
|
using RK4Type = vtkm::worklet::particleadvection::RK4Integrator<GridEvalType>;
|
||||||
|
using Stepper = vtkm::worklet::particleadvection::Stepper<RK4Type, GridEvalType>;
|
||||||
vtkm::FloatDefault stepSize = 0.01f;
|
vtkm::FloatDefault stepSize = 0.01f;
|
||||||
|
|
||||||
const vtkm::Id3 dims(5, 5, 5);
|
const vtkm::Id3 dims(5, 5, 5);
|
||||||
@ -665,7 +673,7 @@ void TestWorkletsBasic()
|
|||||||
for (auto& ds : dataSets)
|
for (auto& ds : dataSets)
|
||||||
{
|
{
|
||||||
GridEvalType eval(ds, velocities);
|
GridEvalType eval(ds, velocities);
|
||||||
RK4Type rk4(eval, stepSize);
|
Stepper rk4(eval, stepSize);
|
||||||
|
|
||||||
vtkm::Id maxSteps = 83;
|
vtkm::Id maxSteps = 83;
|
||||||
std::vector<std::string> workletTypes = { "particleAdvection", "streamline" };
|
std::vector<std::string> workletTypes = { "particleAdvection", "streamline" };
|
||||||
@ -825,6 +833,7 @@ void TestParticleAdvectionFile(const std::string& fname,
|
|||||||
using FieldType = vtkm::worklet::particleadvection::VelocityField<FieldHandle>;
|
using FieldType = vtkm::worklet::particleadvection::VelocityField<FieldHandle>;
|
||||||
using GridEvalType = vtkm::worklet::particleadvection::GridEvaluator<FieldType>;
|
using GridEvalType = vtkm::worklet::particleadvection::GridEvaluator<FieldType>;
|
||||||
using RK4Type = vtkm::worklet::particleadvection::RK4Integrator<GridEvalType>;
|
using RK4Type = vtkm::worklet::particleadvection::RK4Integrator<GridEvalType>;
|
||||||
|
using Stepper = vtkm::worklet::particleadvection::Stepper<RK4Type, GridEvalType>;
|
||||||
|
|
||||||
VTKM_TEST_ASSERT(ds.HasField("vec"), "Data set missing a field named 'vec'");
|
VTKM_TEST_ASSERT(ds.HasField("vec"), "Data set missing a field named 'vec'");
|
||||||
vtkm::cont::Field& field = ds.GetField("vec");
|
vtkm::cont::Field& field = ds.GetField("vec");
|
||||||
@ -842,7 +851,7 @@ void TestParticleAdvectionFile(const std::string& fname,
|
|||||||
|
|
||||||
FieldType velocities(fieldArray);
|
FieldType velocities(fieldArray);
|
||||||
GridEvalType eval(ds.GetCoordinateSystem(), ds.GetCellSet(), velocities);
|
GridEvalType eval(ds.GetCoordinateSystem(), ds.GetCellSet(), velocities);
|
||||||
RK4Type rk4(eval, stepSize);
|
Stepper rk4(eval, stepSize);
|
||||||
|
|
||||||
for (int i = 0; i < 2; i++)
|
for (int i = 0; i < 2; i++)
|
||||||
{
|
{
|
||||||
|
@ -19,8 +19,8 @@
|
|||||||
#include <vtkm/cont/testing/Testing.h>
|
#include <vtkm/cont/testing/Testing.h>
|
||||||
#include <vtkm/worklet/ParticleAdvection.h>
|
#include <vtkm/worklet/ParticleAdvection.h>
|
||||||
#include <vtkm/worklet/particleadvection/Field.h>
|
#include <vtkm/worklet/particleadvection/Field.h>
|
||||||
#include <vtkm/worklet/particleadvection/IntegratorBase.h>
|
|
||||||
#include <vtkm/worklet/particleadvection/Particles.h>
|
#include <vtkm/worklet/particleadvection/Particles.h>
|
||||||
|
#include <vtkm/worklet/particleadvection/Stepper.h>
|
||||||
#include <vtkm/worklet/particleadvection/TemporalGridEvaluators.h>
|
#include <vtkm/worklet/particleadvection/TemporalGridEvaluators.h>
|
||||||
|
|
||||||
template <typename ScalarType>
|
template <typename ScalarType>
|
||||||
|
Loading…
Reference in New Issue
Block a user