Commit c597bab5 authored by Kenneth Moreland's avatar Kenneth Moreland Committed by Kitware Robot
Browse files

Merge topic 'virtual-integrator'

f3b923a5 Temporarily remove UnitTestParticleAdvection.cxx from CUDA tests
723792a4 Use new integrators and evaluators for advection
d49c29ce Make particle advection integrators and evaluators ExecObject
a2602183

 Make integrators have a virtual superclass
Acked-by: Kitware Robot's avatarKitware Robot <kwrobot@kitware.com>
Acked-by: default avatarRobert Maynard <robert.maynard@kitware.com>
Merge-request: !1315
parents 56121b10 f3b923a5
......@@ -70,16 +70,12 @@ void RunTest(const std::string& fname,
using FieldType = vtkm::Float32;
using FieldHandle = vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>>;
using FieldPortalConstType =
typename FieldHandle::template ExecutionTypes<DeviceAdapter>::PortalConst;
vtkm::io::reader::BOVDataSetReader rdr(fname);
vtkm::cont::DataSet ds = rdr.ReadDataSet();
using RGEvalType = vtkm::worklet::particleadvection::UniformGridEvaluate<FieldPortalConstType,
FieldType,
DeviceAdapter>;
using RK4RGType = vtkm::worklet::particleadvection::RK4Integrator<RGEvalType, FieldType>;
using RGEvalType = vtkm::worklet::particleadvection::UniformGridEvaluate<FieldHandle>;
using RK4RGType = vtkm::worklet::particleadvection::RK4Integrator<RGEvalType>;
vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>> fieldArray;
ds.GetField(0).GetData().CopyTo(fieldArray);
......
......@@ -52,8 +52,6 @@ void RunTest(vtkm::Id numSteps, vtkm::Float32 stepSize, vtkm::Id advectType)
using DeviceAdapter = VTKM_DEFAULT_DEVICE_ADAPTER_TAG;
using FieldType = vtkm::Float32;
using FieldHandle = vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>>;
using FieldPortalConstType =
typename FieldHandle::template ExecutionTypes<DeviceAdapter>::PortalConst;
// These lines read two datasets, which are BOVs.
// Currently VTKm does not support providing time series datasets
......@@ -73,11 +71,8 @@ void RunTest(vtkm::Id numSteps, vtkm::Float32 stepSize, vtkm::Id advectType)
// The only change in this example and the vanilla particle advection example is
// this example makes use of the TemporalGridEvaluator.
using GridEvaluator =
vtkm::worklet::particleadvection::TemporalGridEvaluator<FieldPortalConstType,
FieldType,
DeviceAdapter>;
using Integrator = vtkm::worklet::particleadvection::EulerIntegrator<GridEvaluator, FieldType>;
using GridEvaluator = vtkm::worklet::particleadvection::TemporalGridEvaluator<FieldHandle>;
using Integrator = vtkm::worklet::particleadvection::EulerIntegrator<GridEvaluator>;
GridEvaluator eval(ds1.GetCoordinateSystem(),
ds1.GetCellSet(0),
......@@ -114,7 +109,7 @@ void RunTest(vtkm::Id numSteps, vtkm::Float32 stepSize, vtkm::Id advectType)
else
{
vtkm::worklet::Streamline streamline;
vtkm::worklet::StreamlineResult<FieldType> res =
vtkm::worklet::StreamlineResult res =
streamline.Run(integrator, seedArray, numSteps, DeviceAdapter());
vtkm::cont::DataSet outData;
vtkm::cont::CoordinateSystem outputCoords("coordinates", res.positions);
......
......@@ -242,21 +242,16 @@ inline VTKM_CONT vtkm::cont::DataSet Lagrangian::DoExecute(
vtkm::cont::ArrayHandleCartesianProduct<AxisHandle, AxisHandle, AxisHandle>;
using UniformType = vtkm::cont::ArrayHandleUniformPointCoordinates;
using FieldHandle = vtkm::cont::ArrayHandle<vtkm::Vec<T, 3>, StorageType>;
using FieldPortalConstType =
typename FieldHandle::template ExecutionTypes<DeviceAdapter>::PortalConst;
using PortalType_Position = typename vtkm::cont::ArrayHandle<vtkm::Vec<T, 3>>::PortalControl;
using PortalType_DoublePosition =
typename vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Float64, 3>>::PortalControl;
vtkm::worklet::ParticleAdvection particleadvection;
vtkm::worklet::ParticleAdvectionResult<T> res;
vtkm::worklet::ParticleAdvectionResult res;
if (coords.GetData().IsType<RectilinearType>())
{
using RectilinearGridEvalType = vtkm::worklet::particleadvection::
RectilinearGridEvaluate<FieldPortalConstType, T, DeviceAdapter, StorageType>;
using RectilinearGridEvalType =
vtkm::worklet::particleadvection::RectilinearGridEvaluate<FieldHandle>;
using RK4IntegratorType =
vtkm::worklet::particleadvection::RK4Integrator<RectilinearGridEvalType, T>;
vtkm::worklet::particleadvection::RK4Integrator<RectilinearGridEvalType>;
/*
* If Euler step is preferred.
using EulerIntegratorType = vtkm::worklet::particleadvection::EulerIntegrator<RectilinearGridEvalType, T>;
......@@ -271,10 +266,8 @@ inline VTKM_CONT vtkm::cont::DataSet Lagrangian::DoExecute(
}
else if (coords.GetData().IsType<UniformType>())
{
using UniformGridEvalType = vtkm::worklet::particleadvection::
UniformGridEvaluate<FieldPortalConstType, T, DeviceAdapter, StorageType>;
using RK4IntegratorType =
vtkm::worklet::particleadvection::RK4Integrator<UniformGridEvalType, T>;
using UniformGridEvalType = vtkm::worklet::particleadvection::UniformGridEvaluate<FieldHandle>;
using RK4IntegratorType = vtkm::worklet::particleadvection::RK4Integrator<UniformGridEvalType>;
/*
* If Euler step is preferred.
using EulerIntegratorType = vtkm::worklet::particleadvection::EulerIntegrator<UniformGridEvalType, T>;
......@@ -292,17 +285,14 @@ inline VTKM_CONT vtkm::cont::DataSet Lagrangian::DoExecute(
std::cout << "Data set type is not rectilinear or uniform." << std::endl;
}
vtkm::cont::ArrayHandle<vtkm::Vec<T, 3>> particle_positions = res.positions;
vtkm::cont::ArrayHandle<vtkm::Id> particle_stepstaken = res.stepsTaken;
using PortalType_Position = typename vtkm::cont::ArrayHandle<vtkm::Vec<T, 3>>::PortalControl;
using PortalType_ID = typename vtkm::cont::ArrayHandle<vtkm::Id>::PortalControl;
auto particle_positions = res.positions;
auto particle_stepstaken = res.stepsTaken;
PortalType_DoublePosition start_position = BasisParticlesOriginal.GetPortalControl();
PortalType_Position end_position = particle_positions.GetPortalControl();
auto start_position = BasisParticlesOriginal.GetPortalControl();
auto end_position = particle_positions.GetPortalControl();
PortalType_ID portal_stepstaken = particle_stepstaken.GetPortalControl();
PortalType_ID portal_validity = BasisParticlesValidity.GetPortalControl();
auto portal_stepstaken = particle_stepstaken.GetPortalControl();
auto portal_validity = BasisParticlesValidity.GetPortalControl();
vtkm::cont::DataSet outputData;
vtkm::cont::DataSetBuilderExplicit dataSetBuilder;
......
......@@ -23,6 +23,7 @@
#include <vtkm/filter/FilterDataSetWithField.h>
#include <vtkm/worklet/ParticleAdvection.h>
#include <vtkm/worklet/particleadvection/Integrators.h>
namespace vtkm
{
......@@ -39,7 +40,7 @@ public:
Streamline();
VTKM_CONT
void SetStepSize(vtkm::Float64 s) { this->StepSize = s; }
void SetStepSize(vtkm::worklet::particleadvection::ScalarType s) { this->StepSize = s; }
VTKM_CONT
void SetNumberOfSteps(vtkm::Id n) { this->NumberOfSteps = n; }
......@@ -66,7 +67,7 @@ public:
private:
vtkm::worklet::Streamline Worklet;
vtkm::Float64 StepSize;
vtkm::worklet::particleadvection::ScalarType StepSize;
vtkm::Id NumberOfSteps;
vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::FloatDefault, 3>> Seeds;
};
......
......@@ -87,18 +87,15 @@ inline VTKM_CONT vtkm::cont::DataSet Streamline::DoExecute(
//todo: add check for rectilinear.
using FieldHandle = vtkm::cont::ArrayHandle<vtkm::Vec<T, 3>, StorageType>;
using FieldPortalConstType =
typename FieldHandle::template ExecutionTypes<DeviceAdapter>::PortalConst;
using RGEvalType = vtkm::worklet::particleadvection::
UniformGridEvaluate<FieldPortalConstType, T, DeviceAdapter, StorageType>;
using RK4RGType = vtkm::worklet::particleadvection::RK4Integrator<RGEvalType, T>;
using RGEvalType = vtkm::worklet::particleadvection::UniformGridEvaluate<FieldHandle>;
using RK4RGType = vtkm::worklet::particleadvection::RK4Integrator<RGEvalType>;
//RGEvalType eval(input.GetCoordinateSystem(), input.GetCellSet(0), field);
RGEvalType eval(coords, cells, field);
RK4RGType rk4(eval, static_cast<T>(this->StepSize));
vtkm::worklet::Streamline streamline;
vtkm::worklet::StreamlineResult<T> res;
vtkm::worklet::StreamlineResult res;
vtkm::cont::ArrayHandle<vtkm::Vec<T, 3>> seedArray;
vtkm::cont::DeviceAdapterAlgorithm<DeviceAdapter>::Copy(this->Seeds, seedArray);
......
......@@ -31,9 +31,10 @@ namespace vtkm
namespace worklet
{
template <typename FieldType>
struct ParticleAdvectionResult
{
using ScalarType = vtkm::worklet::particleadvection::ScalarType;
ParticleAdvectionResult()
: positions()
, status()
......@@ -42,7 +43,7 @@ struct ParticleAdvectionResult
{
}
ParticleAdvectionResult(const vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>>& pos,
ParticleAdvectionResult(const vtkm::cont::ArrayHandle<vtkm::Vec<ScalarType, 3>>& pos,
const vtkm::cont::ArrayHandle<vtkm::Id>& stat,
const vtkm::cont::ArrayHandle<vtkm::Id>& steps)
: positions(pos)
......@@ -51,10 +52,10 @@ struct ParticleAdvectionResult
{
}
ParticleAdvectionResult(const vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>>& pos,
ParticleAdvectionResult(const vtkm::cont::ArrayHandle<vtkm::Vec<ScalarType, 3>>& pos,
const vtkm::cont::ArrayHandle<vtkm::Id>& stat,
const vtkm::cont::ArrayHandle<vtkm::Id>& steps,
const vtkm::cont::ArrayHandle<FieldType>& timeArray)
const vtkm::cont::ArrayHandle<ScalarType>& timeArray)
: positions(pos)
, status(stat)
, stepsTaken(steps)
......@@ -62,40 +63,39 @@ struct ParticleAdvectionResult
{
}
vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>> positions;
vtkm::cont::ArrayHandle<vtkm::Vec<ScalarType, 3>> positions;
vtkm::cont::ArrayHandle<vtkm::Id> status;
vtkm::cont::ArrayHandle<vtkm::Id> stepsTaken;
vtkm::cont::ArrayHandle<FieldType> times;
vtkm::cont::ArrayHandle<ScalarType> times;
};
class ParticleAdvection
{
public:
using ScalarType = vtkm::worklet::particleadvection::ScalarType;
ParticleAdvection() {}
template <typename IntegratorType,
typename FieldType,
typename PointStorage,
typename DeviceAdapter>
ParticleAdvectionResult<FieldType> Run(
const IntegratorType& it,
const vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>, PointStorage>& pts,
const vtkm::Id& nSteps,
DeviceAdapter tag)
ParticleAdvectionResult Run(const IntegratorType& it,
vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>, PointStorage>& pts,
const vtkm::Id& nSteps,
DeviceAdapter tag)
{
vtkm::Id numSeeds = static_cast<vtkm::Id>(pts.GetNumberOfValues());
vtkm::cont::ArrayHandle<vtkm::Id> stepsTaken;
vtkm::cont::ArrayHandle<FieldType> timeArray;
vtkm::cont::ArrayHandle<ScalarType> timeArray;
//Allocate status and steps arrays.
vtkm::cont::ArrayHandleConstant<vtkm::Id> init(0, numSeeds);
stepsTaken.Allocate(numSeeds);
vtkm::cont::ArrayCopy(init, stepsTaken, tag);
//Allocate memory to store the time for temporal integration.
vtkm::cont::ArrayHandleConstant<FieldType> time(0, numSeeds);
timeArray.Allocate(numSeeds);
vtkm::cont::ArrayHandleConstant<ScalarType> time(0, numSeeds);
vtkm::cont::ArrayCopy(time, timeArray, tag);
return Run(it, pts, stepsTaken, timeArray, nSteps, tag);
......@@ -105,37 +105,32 @@ public:
typename FieldType,
typename PointStorage,
typename DeviceAdapter>
ParticleAdvectionResult<FieldType> Run(
const IntegratorType& it,
const vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>, PointStorage>& pts,
vtkm::cont::ArrayHandle<vtkm::Id>& inputSteps,
const vtkm::Id& nSteps,
DeviceAdapter tag)
ParticleAdvectionResult Run(const IntegratorType& it,
vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>, PointStorage>& pts,
vtkm::cont::ArrayHandle<vtkm::Id>& inputSteps,
const vtkm::Id& nSteps,
DeviceAdapter tag)
{
vtkm::Id numSeeds = static_cast<vtkm::Id>(pts.GetNumberOfValues());
vtkm::cont::ArrayHandle<FieldType> timeArray;
vtkm::cont::ArrayHandle<ScalarType> timeArray;
//Allocate memory to store the time for temporal integration.
vtkm::cont::ArrayHandleConstant<FieldType> time(0, numSeeds);
vtkm::cont::ArrayHandleConstant<ScalarType> time(0, numSeeds);
timeArray.Allocate(numSeeds);
vtkm::cont::ArrayCopy(time, timeArray, tag);
return Run(it, pts, inputSteps, timeArray, nSteps, tag);
}
template <typename IntegratorType,
typename FieldType,
typename PointStorage,
typename DeviceAdapter>
ParticleAdvectionResult<FieldType> Run(
const IntegratorType& it,
const vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>, PointStorage>& pts,
vtkm::cont::ArrayHandle<vtkm::Id>& inputSteps,
vtkm::cont::ArrayHandle<FieldType>& inputTime,
const vtkm::Id& nSteps,
DeviceAdapter tag)
template <typename IntegratorType, typename DeviceAdapter>
ParticleAdvectionResult Run(const IntegratorType& it,
vtkm::cont::ArrayHandle<vtkm::Vec<ScalarType, 3>>& pts,
vtkm::cont::ArrayHandle<vtkm::Id>& inputSteps,
vtkm::cont::ArrayHandle<ScalarType>& inputTime,
const vtkm::Id& nSteps,
DeviceAdapter tag)
{
vtkm::worklet::particleadvection::ParticleAdvectionWorklet<IntegratorType, FieldType> worklet;
vtkm::worklet::particleadvection::ParticleAdvectionWorklet<IntegratorType> worklet;
vtkm::Id numSeeds = static_cast<vtkm::Id>(pts.GetNumberOfValues());
......@@ -147,14 +142,31 @@ public:
worklet.Run(it, pts, nSteps, status, inputSteps, inputTime, tag);
//Create output.
ParticleAdvectionResult<FieldType> res(pts, status, inputSteps, inputTime);
return res;
return ParticleAdvectionResult(pts, status, inputSteps, inputTime);
}
template <typename IntegratorType,
typename FieldType,
typename PointStorage,
typename DeviceAdapter>
ParticleAdvectionResult Run(const IntegratorType& it,
vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>, PointStorage>& pts,
vtkm::cont::ArrayHandle<vtkm::Id>& inputSteps,
vtkm::cont::ArrayHandle<ScalarType>& inputTime,
const vtkm::Id& nSteps,
DeviceAdapter tag)
{
vtkm::cont::ArrayHandle<vtkm::Vec<ScalarType, 3>> ptsCopy;
vtkm::cont::ArrayCopy(pts, ptsCopy, tag);
return Run(it, ptsCopy, inputSteps, inputTime, nSteps, tag);
}
};
template <typename FieldType>
struct StreamlineResult
{
using ScalarType = vtkm::worklet::particleadvection::ScalarType;
using VectorType = vtkm::Vec<ScalarType, 3>;
StreamlineResult()
: positions()
, polyLines()
......@@ -164,7 +176,7 @@ struct StreamlineResult
{
}
StreamlineResult(const vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>>& pos,
StreamlineResult(const vtkm::cont::ArrayHandle<vtkm::Vec<ScalarType, 3>>& pos,
const vtkm::cont::CellSetExplicit<>& lines,
const vtkm::cont::ArrayHandle<vtkm::Id>& stat,
const vtkm::cont::ArrayHandle<vtkm::Id>& steps)
......@@ -175,11 +187,11 @@ struct StreamlineResult
{
}
StreamlineResult(const vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>>& pos,
StreamlineResult(const vtkm::cont::ArrayHandle<vtkm::Vec<ScalarType, 3>>& pos,
const vtkm::cont::CellSetExplicit<>& lines,
const vtkm::cont::ArrayHandle<vtkm::Id>& stat,
const vtkm::cont::ArrayHandle<vtkm::Id>& steps,
const vtkm::cont::ArrayHandle<FieldType>& timeArray)
const vtkm::cont::ArrayHandle<ScalarType>& timeArray)
: positions(pos)
, polyLines(lines)
......@@ -189,27 +201,28 @@ struct StreamlineResult
{
}
vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>> positions;
vtkm::cont::ArrayHandle<vtkm::Vec<ScalarType, 3>> positions;
vtkm::cont::CellSetExplicit<> polyLines;
vtkm::cont::ArrayHandle<vtkm::Id> status;
vtkm::cont::ArrayHandle<vtkm::Id> stepsTaken;
vtkm::cont::ArrayHandle<FieldType> times;
vtkm::cont::ArrayHandle<ScalarType> times;
};
class Streamline
{
public:
using ScalarType = vtkm::worklet::particleadvection::ScalarType;
Streamline() {}
template <typename IntegratorType,
typename FieldType,
typename PointStorage,
typename DeviceAdapter>
StreamlineResult<FieldType> Run(
const IntegratorType& it,
const vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>, PointStorage>& seedArray,
const vtkm::Id& nSteps,
DeviceAdapter tag)
StreamlineResult Run(const IntegratorType& it,
vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>, PointStorage>& seedArray,
const vtkm::Id& nSteps,
DeviceAdapter tag)
{
using DeviceAlgorithm = typename vtkm::cont::DeviceAdapterAlgorithm<DeviceAdapter>;
......@@ -227,8 +240,8 @@ public:
DeviceAlgorithm::Copy(zero, steps);
//Allocate memory to store the time for temporal integration.
vtkm::cont::ArrayHandle<FieldType> timeArray;
vtkm::cont::ArrayHandleConstant<FieldType> time(0, numSeeds);
vtkm::cont::ArrayHandle<ScalarType> timeArray;
vtkm::cont::ArrayHandleConstant<ScalarType> time(0, numSeeds);
timeArray.Allocate(numSeeds);
DeviceAlgorithm::Copy(time, timeArray);
......@@ -239,12 +252,11 @@ public:
typename FieldType,
typename PointStorage,
typename DeviceAdapter>
StreamlineResult<FieldType> Run(
const IntegratorType& it,
const vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>, PointStorage>& seedArray,
vtkm::cont::ArrayHandle<vtkm::Id>& inputSteps,
const vtkm::Id& nSteps,
DeviceAdapter tag)
StreamlineResult Run(const IntegratorType& it,
vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>, PointStorage>& seedArray,
vtkm::cont::ArrayHandle<vtkm::Id>& inputSteps,
const vtkm::Id& nSteps,
DeviceAdapter tag)
{
using DeviceAlgorithm = typename vtkm::cont::DeviceAdapterAlgorithm<DeviceAdapter>;
......@@ -257,30 +269,26 @@ public:
DeviceAlgorithm::Copy(statusOK, status);
//Allocate memory to store the time for temporal integration.
vtkm::cont::ArrayHandle<FieldType> timeArray;
vtkm::cont::ArrayHandleConstant<FieldType> time(0, numSeeds);
vtkm::cont::ArrayHandle<ScalarType> timeArray;
vtkm::cont::ArrayHandleConstant<ScalarType> time(0, numSeeds);
timeArray.Allocate(numSeeds);
DeviceAlgorithm::Copy(time, timeArray);
return Run(it, seedArray, inputSteps, timeArray, nSteps, tag);
}
template <typename IntegratorType,
typename FieldType,
typename PointStorage,
typename DeviceAdapter>
StreamlineResult<FieldType> Run(
const IntegratorType& it,
const vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>, PointStorage>& seedArray,
vtkm::cont::ArrayHandle<vtkm::Id>& inputSteps,
vtkm::cont::ArrayHandle<FieldType>& inputTime,
const vtkm::Id& nSteps,
DeviceAdapter tag)
template <typename IntegratorType, typename DeviceAdapter>
StreamlineResult Run(const IntegratorType& it,
vtkm::cont::ArrayHandle<vtkm::Vec<ScalarType, 3>>& seedArray,
vtkm::cont::ArrayHandle<vtkm::Id>& inputSteps,
vtkm::cont::ArrayHandle<ScalarType>& inputTime,
const vtkm::Id& nSteps,
DeviceAdapter tag)
{
using DeviceAlgorithm = typename vtkm::cont::DeviceAdapterAlgorithm<DeviceAdapter>;
vtkm::worklet::particleadvection::StreamlineWorklet<IntegratorType, FieldType> worklet;
vtkm::worklet::particleadvection::StreamlineWorklet<IntegratorType> worklet;
vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>, PointStorage> positions;
vtkm::cont::ArrayHandle<vtkm::Vec<ScalarType, 3>> positions;
vtkm::cont::CellSetExplicit<> polyLines;
//Allocate and initialize status array.
......@@ -292,8 +300,23 @@ public:
worklet.Run(it, seedArray, nSteps, positions, polyLines, status, inputSteps, inputTime, tag);
StreamlineResult<FieldType> res(positions, polyLines, status, inputSteps, inputTime);
return res;
return StreamlineResult(positions, polyLines, status, inputSteps, inputTime);
}
template <typename IntegratorType,
typename FieldType,
typename PointStorage,
typename DeviceAdapter>
StreamlineResult Run(const IntegratorType& it,
vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>, PointStorage>& seedArray,
vtkm::cont::ArrayHandle<vtkm::Id>& inputSteps,
vtkm::cont::ArrayHandle<ScalarType>& inputTime,
const vtkm::Id& nSteps,
DeviceAdapter tag)
{
vtkm::cont::ArrayHandle<vtkm::Vec<ScalarType, 3>> seedCopy;
vtkm::cont::ArrayCopy(seedArray, seedCopy, tag);
return Run(it, seedCopy, inputSteps, inputTime, nSteps, tag);
}
};
}
......
......@@ -22,6 +22,7 @@
#define vtk_m_worklet_particleadvection_ParticleAdvectionWorklets_h
#include <vtkm/Types.h>
#include <vtkm/cont/ArrayCopy.h>
#include <vtkm/cont/ArrayHandle.h>
#include <vtkm/cont/ArrayHandleCast.h>
#include <vtkm/cont/ArrayHandleCounting.h>
......@@ -38,33 +39,36 @@ namespace worklet
namespace particleadvection
{
template <typename IntegratorType, typename FieldType>
class ParticleAdvectWorklet : public vtkm::worklet::WorkletMapField
{
public:
using ControlSignature = void(FieldIn<IdType> idx, ExecObject ic);
using ExecutionSignature = void(_1, _2);
using ControlSignature = void(FieldIn<IdType> idx,
ExecObject integrator,
ExecObject integralCurve);
using ExecutionSignature = void(_1, _2, _3);
using InputDomain = _1;
template <typename IntegralCurveType>
VTKM_EXEC void operator()(const vtkm::Id& idx, IntegralCurveType& ic) const
template <typename IntegratorType, typename IntegralCurveType>
VTKM_EXEC void operator()(const vtkm::Id& idx,
const IntegratorType* integrator,
IntegralCurveType& integralCurve) const
{
vtkm::Vec<FieldType, 3> inpos = ic.GetPos(idx);
vtkm::Vec<FieldType, 3> outpos;
FieldType time = ic.GetTime(idx);
vtkm::Vec<ScalarType, 3> inpos = integralCurve.GetPos(idx);
vtkm::Vec<ScalarType, 3> outpos;
ScalarType time = integralCurve.GetTime(idx);
ParticleStatus status;
while (!ic.Done(idx))
while (!integralCurve.Done(idx))
{
status = integrator.Step(inpos, time, outpos);
status = integrator->Step(inpos, time, outpos);
// If the status is OK, we only need to check if the particle
// has completed the maximum steps required.
if (status == ParticleStatus::STATUS_OK)
{
ic.TakeStep(idx, outpos, status);
integralCurve.TakeStep(idx, outpos, status);
// This is to keep track of the particle's time.
// This is what the Evaluator uses to determine if the particle
// has exited temporal boundary.
ic.SetTime(idx, time);
integralCurve.SetTime(idx, time);
inpos = outpos;
}
// If the particle is at spatial or temporal boundary, take steps to just
......@@ -74,74 +78,53 @@ public:
if (status == ParticleStatus::AT_SPATIAL_BOUNDARY ||
status == ParticleStatus::AT_TEMPORAL_BOUNDARY)
{
vtkm::Id numSteps = ic.GetStep(idx);
status = integrator.PushOutOfBoundary(inpos, numSteps, time, status, outpos);
ic.TakeStep(idx, outpos, status);
ic.SetTime(idx, time);
vtkm::Id numSteps = integralCurve.GetStep(idx);
status = integrator->PushOutOfBoundary(inpos, numSteps, time, status, outpos);
integralCurve.TakeStep(idx, outpos, status);
integralCurve.SetTime(idx, time);
if (status == ParticleStatus::EXITED_SPATIAL_BOUNDARY)
ic.SetExitedSpatialBoundary(idx);
integralCurve.SetExitedSpatialBoundary(idx);
if (status == ParticleStatus::EXITED_TEMPORAL_BOUNDARY)
ic.SetExitedTemporalBoundary(idx);
integralCurve.SetExitedTemporalBoundary(idx);
}
// If the particle has exited spatial boundary, set corresponding status.
else if (status == ParticleStatus::EXITED_SPATIAL_BOUNDARY)
{
ic.TakeStep(idx, outpos, status);
ic.SetExitedSpatialBoundary(idx);
integralCurve.TakeStep(idx, outpos, status);
integralCurve.SetExitedSpatialBoundary(idx);
}
// If the particle has exited temporal boundary, set corresponding status.
else if (status == ParticleStatus::EXITED_TEMPORAL_BOUNDARY)
{
ic.TakeStep(idx, outpos, status);
ic.SetExitedTemporalBoundary(idx);
integralCurve.TakeStep(idx, outpos, status);
integralCurve.SetExitedTemporalBoundary(idx);
}
}