Shuenhoy/vtk-m2
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Add more unittests. Put field into evaluator.

Browse Source
This commit is contained in:
Dave Pugmire 2017-08-01 16:27:03 -04:00
parent 90d0e0e966
commit 33ada346ae
5 changed files with 333 additions and 75 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
vtkm/worklet/ParticleAdvection.h
View File

@ -62,12 +62,10 @@ public:
template <typename IntegratorType,
typename FieldType,
typename PointStorage,
typename FieldStorage,
typename DeviceAdapter>
ParticleAdvectionResult<FieldType> Run(
const IntegratorType& it,
const vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>, PointStorage>& pts,
vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>, FieldStorage> fieldArray,
const vtkm::Id& nSteps,
const DeviceAdapter&)
{
@ -76,8 +74,8 @@ public:
DeviceAdapter>
worklet;
vtkm::cont::ArrayHandle<vtkm::Id, FieldStorage> stepsTaken, status;
worklet.Run(it, pts, fieldArray, nSteps, status, stepsTaken);
vtkm::cont::ArrayHandle<vtkm::Id> stepsTaken, status;
worklet.Run(it, pts, nSteps, status, stepsTaken);
//Create output.
ParticleAdvectionResult<FieldType> res(pts, status, stepsTaken);

133
vtkm/worklet/particleadvection/GridEvaluators.h
View File

@ -36,20 +36,6 @@ namespace worklet
namespace particleadvection
{
//Base class for all grid evaluators
class GridEvaluate
{
public:
enum Status
{
OK = 0,
OUTSIDE_SPATIAL,
OUTSIDE_TEMPORAL,
FAIL
};
};
// Constant vector
template <typename PortalType, typename FieldType>
class ConstantField
@ -63,9 +49,7 @@ public:
}
VTKM_EXEC_CONT
bool Evaluate(const vtkm::Vec<FieldType, 3>& pos,
const PortalType& vtkmNotUsed(vecData),
vtkm::Vec<FieldType, 3>& out) const
bool Evaluate(const vtkm::Vec<FieldType, 3>& pos, vtkm::Vec<FieldType, 3>& out) const
{
if (!bounds.Contains(pos))
return false;
@ -82,7 +66,6 @@ private:
};
// Circular Orbit.
template <typename PortalType>
class AnalyticalOrbitEvaluate
{
public:
@ -94,7 +77,6 @@ public:
template <typename FieldType>
VTKM_EXEC_CONT bool Evaluate(const vtkm::Vec<FieldType, 3>& pos,
const PortalType& vtkmNotUsed(vecData),
vtkm::Vec<FieldType, 3>& out) const
{
if (!bounds.Contains(pos))
@ -112,14 +94,52 @@ private:
vtkm::Bounds bounds;
};
//Uniform Grid Evaluator
template <typename PortalType, typename FieldType>
template <typename PortalType, typename FieldType, typename DeviceAdapterTag>
class UniformGridEvaluate
{
typedef vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>> FieldHandle;
public:
VTKM_CONT
UniformGridEvaluate() {}
VTKM_CONT
UniformGridEvaluate(const vtkm::cont::CoordinateSystem& coords,
const vtkm::cont::DynamicCellSet& cellSet,
const FieldHandle& vectorField)
{
vectors = vectorField.PrepareForInput(DeviceAdapterTag());
typedef vtkm::cont::ArrayHandleUniformPointCoordinates UniformType;
typedef vtkm::cont::CellSetStructured<3> StructuredType;
if (!coords.GetData().IsSameType(UniformType()))
throw vtkm::cont::ErrorInternal("Coordinates are not uniform.");
if (!cellSet.IsSameType(StructuredType()))
throw vtkm::cont::ErrorInternal("Cells are not 3D structured.");
bounds = coords.GetBounds();
vtkm::cont::CellSetStructured<3> cells;
cellSet.CopyTo(cells);
dims = cells.GetSchedulingRange(vtkm::TopologyElementTagPoint());
vtkm::Vec<FieldType, 3> castDims = dims;
spacing[0] = static_cast<FieldType>((bounds.X.Max - bounds.X.Min) / (castDims[0] - 1));
spacing[1] = static_cast<FieldType>((bounds.Y.Max - bounds.Y.Min) / (castDims[1] - 1));
spacing[2] = static_cast<FieldType>((bounds.Z.Max - bounds.Z.Min) / (castDims[2] - 1));
oldMin[0] =
static_cast<FieldType>(bounds.X.Min / ((bounds.X.Max - bounds.X.Min) / castDims[0]));
oldMin[1] =
static_cast<FieldType>(bounds.Y.Min / ((bounds.Y.Max - bounds.Y.Min) / castDims[1]));
oldMin[2] =
static_cast<FieldType>(bounds.Z.Min / ((bounds.Z.Max - bounds.Z.Min) / castDims[2]));
planeSize = dims[0] * dims[1];
rowSize = dims[0];
}
VTKM_CONT
UniformGridEvaluate(const vtkm::cont::DataSet& ds)
{
@ -151,9 +171,7 @@ public:
}
VTKM_EXEC_CONT
bool Evaluate(const vtkm::Vec<FieldType, 3>& pos,
const PortalType& vecData,
vtkm::Vec<FieldType, 3>& out) const
bool Evaluate(const vtkm::Vec<FieldType, 3>& pos, vtkm::Vec<FieldType, 3>& out) const
{
if (!bounds.Contains(pos))
return false;
@ -187,14 +205,14 @@ public:
// Get the vecdata at the eight corners
vtkm::Vec<FieldType, 3> v000, v001, v010, v011, v100, v101, v110, v111;
v000 = vecData.Get(idx000[2] * planeSize + idx000[1] * rowSize + idx000[0]);
v001 = vecData.Get(idx001[2] * planeSize + idx001[1] * rowSize + idx001[0]);
v010 = vecData.Get(idx010[2] * planeSize + idx010[1] * rowSize + idx010[0]);
v011 = vecData.Get(idx011[2] * planeSize + idx011[1] * rowSize + idx011[0]);
v100 = vecData.Get(idx100[2] * planeSize + idx100[1] * rowSize + idx100[0]);
v101 = vecData.Get(idx101[2] * planeSize + idx101[1] * rowSize + idx101[0]);
v110 = vecData.Get(idx110[2] * planeSize + idx110[1] * rowSize + idx110[0]);
v111 = vecData.Get(idx111[2] * planeSize + idx111[1] * rowSize + idx111[0]);
v000 = vectors.Get(idx000[2] * planeSize + idx000[1] * rowSize + idx000[0]);
v001 = vectors.Get(idx001[2] * planeSize + idx001[1] * rowSize + idx001[0]);
v010 = vectors.Get(idx010[2] * planeSize + idx010[1] * rowSize + idx010[0]);
v011 = vectors.Get(idx011[2] * planeSize + idx011[1] * rowSize + idx011[0]);
v100 = vectors.Get(idx100[2] * planeSize + idx100[1] * rowSize + idx100[0]);
v101 = vectors.Get(idx101[2] * planeSize + idx101[1] * rowSize + idx101[0]);
v110 = vectors.Get(idx110[2] * planeSize + idx110[1] * rowSize + idx110[0]);
v111 = vectors.Get(idx111[2] * planeSize + idx111[1] * rowSize + idx111[0]);
// Interpolation in X
vtkm::Vec<FieldType, 3> v00, v01, v10, v11;
@ -236,6 +254,7 @@ public:
private:
vtkm::Bounds bounds;
vtkm::Id3 dims;
PortalType vectors;
vtkm::Id planeSize;
vtkm::Id rowSize;
vtkm::Vec<FieldType, 3> spacing;
@ -244,10 +263,39 @@ private:
template <typename PortalType, typename FieldType>
template <typename PortalType, typename FieldType, typename DeviceAdapterTag>
class RectilinearGridEvaluate
{
typedef vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>> FieldHandle;
public:
VTKM_CONT
RectilinearGridEvaluate(const vtkm::cont::CoordinateSystem& coords,
const vtkm::cont::DynamicCellSet& cellSet,
const FieldHandle& vectorField)
{
typedef vtkm::cont::CellSetStructured<3> StructuredType;
if (!coords.GetData().IsSameType(RectilinearType()))
throw vtkm::cont::ErrorInternal("Coordinates are not rectilinear.");
if (!cellSet.IsSameType(StructuredType()))
throw vtkm::cont::ErrorInternal("Cells are not 3D structured.");
vectors = vectorField.PrepareForInput(DeviceAdapterTag());
bounds = coords.GetBounds();
vtkm::cont::CellSetStructured<3> cells;
cellSet.CopyTo(cells);
dims = cells.GetSchedulingRange(vtkm::TopologyElementTagPoint());
planeSize = dims[0] * dims[1];
rowSize = dims[0];
RectilinearType gridPoints = coords.GetData().Cast<RectilinearType>();
xAxis = gridPoints.GetPortalConstControl().GetFirstPortal();
yAxis = gridPoints.GetPortalConstControl().GetSecondPortal();
zAxis = gridPoints.GetPortalConstControl().GetThirdPortal();
}
VTKM_CONT
RectilinearGridEvaluate(const vtkm::cont::DataSet& dataset)
{
@ -278,9 +326,7 @@ public:
}
VTKM_EXEC_CONT
bool Evaluate(const vtkm::Vec<FieldType, 3>& pos,
const PortalType& vecData,
vtkm::Vec<FieldType, 3>& out) const
bool Evaluate(const vtkm::Vec<FieldType, 3>& pos, vtkm::Vec<FieldType, 3>& out) const
{
if (!bounds.Contains(pos))
return false;
@ -337,14 +383,14 @@ public:
// Get the vecdata at the eight corners
vtkm::Vec<FieldType, 3> v000, v001, v010, v011, v100, v101, v110, v111;
v000 = vecData.Get(idx000[2] * planeSize + idx000[1] * rowSize + idx000[0]);
v001 = vecData.Get(idx001[2] * planeSize + idx001[1] * rowSize + idx001[0]);
v010 = vecData.Get(idx010[2] * planeSize + idx010[1] * rowSize + idx010[0]);
v011 = vecData.Get(idx011[2] * planeSize + idx011[1] * rowSize + idx011[0]);
v100 = vecData.Get(idx100[2] * planeSize + idx100[1] * rowSize + idx100[0]);
v101 = vecData.Get(idx101[2] * planeSize + idx101[1] * rowSize + idx101[0]);
v110 = vecData.Get(idx110[2] * planeSize + idx110[1] * rowSize + idx110[0]);
v111 = vecData.Get(idx111[2] * planeSize + idx111[1] * rowSize + idx111[0]);
v000 = vectors.Get(idx000[2] * planeSize + idx000[1] * rowSize + idx000[0]);
v001 = vectors.Get(idx001[2] * planeSize + idx001[1] * rowSize + idx001[0]);
v010 = vectors.Get(idx010[2] * planeSize + idx010[1] * rowSize + idx010[0]);
v011 = vectors.Get(idx011[2] * planeSize + idx011[1] * rowSize + idx011[0]);
v100 = vectors.Get(idx100[2] * planeSize + idx100[1] * rowSize + idx100[0]);
v101 = vectors.Get(idx101[2] * planeSize + idx101[1] * rowSize + idx101[0]);
v110 = vectors.Get(idx110[2] * planeSize + idx110[1] * rowSize + idx110[0]);
v111 = vectors.Get(idx111[2] * planeSize + idx111[1] * rowSize + idx111[0]);
// Interpolation in X
vtkm::Vec<FieldType, 3> v00, v01, v10, v11;
@ -395,6 +441,7 @@ private:
typename AxisHandle::PortalConstControl zAxis;
vtkm::Bounds bounds;
vtkm::Id3 dims;
PortalType vectors;
vtkm::Id planeSize;
vtkm::Id rowSize;

9
vtkm/worklet/particleadvection/Integrators.h
View File

@ -49,15 +49,12 @@ public:
{
}
template <typename PortalType>
VTKM_EXEC bool Step(const vtkm::Vec<FieldType, 3>& pos,
const PortalType& field,
vtkm::Vec<FieldType, 3>& out) const
VTKM_EXEC bool Step(const vtkm::Vec<FieldType, 3>& pos, vtkm::Vec<FieldType, 3>& out) const
{
vtkm::Vec<FieldType, 3> k1, k2, k3, k4;
if (f.Evaluate(pos, field, k1) && f.Evaluate(pos + h_2 * k1, field, k2) &&
f.Evaluate(pos + h_2 * k2, field, k3) && f.Evaluate(pos + h * k3, field, k4))
if (f.Evaluate(pos, k1) && f.Evaluate(pos + h_2 * k1, k2) && f.Evaluate(pos + h_2 * k2, k3) &&
f.Evaluate(pos + h * k3, k4))
{
out = pos + h / 6.0f * (k1 + 2 * k2 + 2 * k3 + k4);
return true;

18
vtkm/worklet/particleadvection/ParticleAdvectionWorklets.h
View File

@ -42,10 +42,6 @@ template <typename IntegratorType, typename FieldType, typename DeviceAdapterTag
class ParticleAdvectWorklet : public vtkm::worklet::WorkletMapField
{
public:
typedef vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>> FieldHandle;
typedef typename FieldHandle::template ExecutionTypes<DeviceAdapterTag>::PortalConst
FieldPortalConstType;
typedef void ControlSignature(FieldIn<IdType> idx, ExecObject ic);
typedef void ExecutionSignature(_1, _2);
typedef _1 InputDomain;
@ -58,7 +54,7 @@ public:
while (!ic.Done(idx))
{
if (integrator.Step(p, field, p2))
if (integrator.Step(p, p2))
{
ic.TakeStep(idx, p2);
p = p2;
@ -70,14 +66,12 @@ public:
}
}
ParticleAdvectWorklet(const IntegratorType& it, const FieldPortalConstType& f)
ParticleAdvectWorklet(const IntegratorType& it)
: integrator(it)
, field(f)
{
}
IntegratorType integrator;
FieldPortalConstType field;
};
@ -85,9 +79,6 @@ template <typename IntegratorType, typename FieldType, typename DeviceAdapterTag
class ParticleAdvectionWorklet
{
public:
typedef vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>> FieldHandle;
typedef typename FieldHandle::template ExecutionTypes<DeviceAdapterTag>::PortalConst
FieldPortalConstType;
typedef vtkm::worklet::particleadvection::ParticleAdvectWorklet<IntegratorType,
FieldType,
DeviceAdapterTag>
@ -98,7 +89,6 @@ public:
template <typename PointStorage, typename FieldStorage>
void Run(const IntegratorType& it,
const vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>, PointStorage>& pts,
const vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>, FieldStorage> fieldArray,
const vtkm::Id& nSteps,
vtkm::cont::ArrayHandle<vtkm::Id, FieldStorage>& statusArray,
vtkm::cont::ArrayHandle<vtkm::Id, FieldStorage>& stepsTaken)
@ -106,7 +96,6 @@ public:
integrator = it;
seedArray = pts;
maxSteps = nSteps;
field = fieldArray.PrepareForInput(DeviceAdapterTag());
run(statusArray, stepsTaken);
}
@ -137,7 +126,7 @@ private:
vtkm::cont::ArrayHandleIndex idxArray(numSeeds);
ParticleType particles(seedArray, stepsTaken, statusArray, maxSteps);
ParticleAdvectWorkletType particleWorklet(integrator, field);
ParticleAdvectWorkletType particleWorklet(integrator);
ParticleWorkletDispatchType particleWorkletDispatch(particleWorklet);
particleWorkletDispatch.Invoke(idxArray, particles);
}
@ -146,7 +135,6 @@ private:
vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>> seedArray;
vtkm::cont::DataSet ds;
vtkm::Id maxSteps;
FieldPortalConstType field;
};

242
vtkm/worklet/testing/UnitTestParticleAdvection.cxx
View File

@ -20,6 +20,7 @@
#include <vtkm/cont/ArrayHandle.h>
#include <vtkm/cont/DataSet.h>
#include <vtkm/cont/DataSetBuilderRectilinear.h>
#include <vtkm/cont/DataSetBuilderUniform.h>
#include <vtkm/cont/DeviceAdapter.h>
#include <vtkm/cont/testing/Testing.h>
@ -29,6 +30,9 @@
#include <vtkm/worklet/particleadvection/Integrators.h>
#include <vtkm/worklet/particleadvection/Particles.h>
#include <typeinfo>
namespace
{
@ -99,7 +103,223 @@ vtkm::Float32 vecData[125 * 3] = {
};
}
void TestParticleAdvection()
template <typename FieldType>
void RandomPoint(const vtkm::Bounds& bounds, vtkm::Vec<FieldType, 3>& p)
{
FieldType rx = static_cast<FieldType>(rand()) / static_cast<FieldType>(RAND_MAX);
FieldType ry = static_cast<FieldType>(rand()) / static_cast<FieldType>(RAND_MAX);
FieldType rz = static_cast<FieldType>(rand()) / static_cast<FieldType>(RAND_MAX);
p[0] = static_cast<FieldType>(bounds.X.Min + rx * bounds.X.Length());
p[1] = static_cast<FieldType>(bounds.Y.Min + ry * bounds.Y.Length());
p[2] = static_cast<FieldType>(bounds.Z.Min + rz * bounds.Z.Length());
}
template <typename FieldType>
vtkm::cont::DataSet CreateUniformDataSet(const vtkm::Bounds& bounds, const vtkm::Id3& dims)
{
vtkm::Vec<FieldType, 3> origin(bounds.X.Min, bounds.Y.Min, bounds.Z.Min);
vtkm::Vec<FieldType, 3> spacing(bounds.X.Length() / (dims[0] - 1),
bounds.Y.Length() / (dims[1] - 1),
bounds.Z.Length() / (dims[2] - 1));
vtkm::cont::DataSetBuilderUniform dataSetBuilder;
vtkm::cont::DataSet ds = dataSetBuilder.Create(dims, origin, spacing);
return ds;
}
template <typename FieldType>
vtkm::cont::DataSet CreateRectilinearDataSet(const vtkm::Bounds& bounds, const vtkm::Id3& dims)
{
vtkm::cont::DataSetBuilderRectilinear dataSetBuilder;
std::vector<FieldType> xvals, yvals, zvals;
vtkm::Vec<FieldType, 3> spacing(bounds.X.Length() / (dims[0] - 1),
bounds.Y.Length() / (dims[1] - 1),
bounds.Z.Length() / (dims[2] - 1));
xvals.resize(dims[0]);
xvals[0] = bounds.X.Min;
for (vtkm::Id i = 1; i < dims[0]; i++)
xvals[i] = xvals[i - 1] + spacing[0];
yvals.resize(dims[1]);
yvals.resize(dims[1]);
yvals[0] = bounds.Y.Min;
for (vtkm::Id i = 1; i < dims[1]; i++)
yvals[i] = yvals[i - 1] + spacing[1];
zvals.resize(dims[2]);
zvals[0] = bounds.Z.Min;
for (vtkm::Id i = 1; i < dims[2]; i++)
zvals[i] = zvals[i - 1] + spacing[2];
vtkm::cont::DataSet ds = dataSetBuilder.Create(xvals, yvals, zvals);
return ds;
}
template <typename FieldType>
void CreateConstantVectorField(vtkm::Id num,
const vtkm::Vec<FieldType, 3>& vec,
vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>>& vecField)
{
typedef VTKM_DEFAULT_DEVICE_ADAPTER_TAG DeviceAdapter;
typedef typename vtkm::cont::DeviceAdapterAlgorithm<DeviceAdapter> DeviceAlgorithm;
vtkm::cont::ArrayHandleConstant<vtkm::Vec<FieldType, 3>> vecConst;
vecConst = vtkm::cont::make_ArrayHandleConstant(vec, num);
DeviceAlgorithm::Copy(vecConst, vecField);
}
template <typename EvalType, typename FieldType>
void ValidateEvaluator(const EvalType& eval,
const vtkm::Vec<FieldType, 3>& p,
const vtkm::Vec<FieldType, 3>& vec,
const std::string& msg)
{
vtkm::Vec<FieldType, 3> result;
bool val = eval.Evaluate(p, result);
VTKM_TEST_ASSERT(val, "Error in evaluator for " + msg);
VTKM_TEST_ASSERT(result == vec, "Error in evaluator result for " + msg);
}
template <typename IntegratorType, typename FieldType>
void ValidateIntegrator(const IntegratorType& integrator,
const vtkm::Vec<FieldType, 3>& p,
const vtkm::Vec<FieldType, 3>& q,
const std::string& msg)
{
vtkm::Vec<FieldType, 3> result;
bool val = integrator.Step(p, result);
VTKM_TEST_ASSERT(val, "Error in integrator for " + msg);
VTKM_TEST_ASSERT(result == q, "Error in integrator result for " + msg);
}
void TestEvaluators()
{
typedef VTKM_DEFAULT_DEVICE_ADAPTER_TAG DeviceAdapter;
typedef vtkm::Float32 FieldType;
typedef vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>> FieldHandle;
typedef FieldHandle::template ExecutionTypes<DeviceAdapter>::PortalConst FieldPortalConstType;
typedef typename vtkm::cont::DeviceAdapterAlgorithm<DeviceAdapter> DeviceAlgorithm;
//Constant field evaluator and RK4 integrator.
typedef vtkm::worklet::particleadvection::ConstantField<FieldPortalConstType, FieldType>
CEvalType;
typedef vtkm::worklet::particleadvection::RK4Integrator<CEvalType, FieldType> RK4CType;
//Uniform grid evaluator and RK4 integrator.
typedef vtkm::worklet::particleadvection::UniformGridEvaluate<FieldPortalConstType,
FieldType,
DeviceAdapter>
UniformEvalType;
typedef vtkm::worklet::particleadvection::RK4Integrator<UniformEvalType, FieldType>
RK4UniformType;
//Rectilinear grid evaluator and RK4 integrator.
typedef vtkm::worklet::particleadvection::RectilinearGridEvaluate<FieldPortalConstType,
FieldType,
DeviceAdapter>
RectilinearEvalType;
typedef vtkm::worklet::particleadvection::RK4Integrator<RectilinearEvalType, FieldType>
RK4RectilinearType;
std::vector<vtkm::Vec<FieldType, 3>> vecs;
vecs.push_back(vtkm::Vec<FieldType, 3>(1, 0, 0));
vecs.push_back(vtkm::Vec<FieldType, 3>(0, 1, 0));
vecs.push_back(vtkm::Vec<FieldType, 3>(0, 0, 1));
vecs.push_back(vtkm::Vec<FieldType, 3>(1, 1, 0));
vecs.push_back(vtkm::Vec<FieldType, 3>(0, 1, 1));
vecs.push_back(vtkm::Vec<FieldType, 3>(1, 0, 1));
vecs.push_back(vtkm::Vec<FieldType, 3>(1, 1, 1));
std::vector<vtkm::Bounds> bounds;
bounds.push_back(vtkm::Bounds(0, 10, 0, 10, 0, 10));
bounds.push_back(vtkm::Bounds(-1, 1, -1, 1, -1, 1));
bounds.push_back(vtkm::Bounds(0, 1, 0, 1, -1, 1));
std::vector<vtkm::Id3> dims;
dims.push_back(vtkm::Id3(5, 5, 5));
dims.push_back(vtkm::Id3(10, 5, 5));
dims.push_back(vtkm::Id3(10, 5, 5));
std::vector<FieldType> steps;
steps.push_back(0.01f);
// steps.push_back(0.05f);
srand(314);
for (std::size_t d = 0; d < dims.size(); d++)
{
vtkm::Id3 dim = dims[d];
for (std::size_t i = 0; i < vecs.size(); i++)
{
vtkm::Vec<FieldType, 3> vec = vecs[i];
for (std::size_t j = 0; j < bounds.size(); j++)
{
//Create uniform and rectilinear data.
vtkm::cont::DataSet uniformData, rectData;
vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>> vecField;
uniformData = CreateUniformDataSet<FieldType>(bounds[j], dim);
rectData = CreateRectilinearDataSet<FieldType>(bounds[j], dim);
CreateConstantVectorField(dim[0] * dim[1] * dim[2], vec, vecField);
for (std::size_t s = 0; s < steps.size(); s++)
{
FieldType stepSize = steps[s];
//Constant vector field evaluator.
CEvalType constEval(bounds[j], vecs[i]);
RK4CType constRK4(constEval, stepSize);
//Uniform vector field evaluator
UniformEvalType uniformEval(
uniformData.GetCoordinateSystem(), uniformData.GetCellSet(0), vecField);
RK4UniformType uniformRK4(uniformEval, stepSize);
//Rectilinear grid evaluator.
RectilinearEvalType rectEval(
rectData.GetCoordinateSystem(), rectData.GetCellSet(0), vecField);
RK4RectilinearType rectRK4(rectEval, stepSize);
//Create a bunch of random points in the bounds.
for (int k = 0; k < 38; k++)
{
//Generate points 2 steps inside the bounding box.
vtkm::Bounds interiorBounds = bounds[j];
interiorBounds.X.Min += 2 * stepSize;
interiorBounds.Y.Min += 2 * stepSize;
interiorBounds.Z.Min += 2 * stepSize;
interiorBounds.X.Max -= 2 * stepSize;
interiorBounds.Y.Max -= 2 * stepSize;
interiorBounds.Z.Max -= 2 * stepSize;
vtkm::Vec<FieldType, 3> p;
RandomPoint<FieldType>(interiorBounds, p);
//Test the result for the evaluator
ValidateEvaluator(constEval, p, vec, "constant vector evaluator");
ValidateEvaluator(uniformEval, p, vec, "uniform evaluator");
ValidateEvaluator(rectEval, p, vec, "rectilinear evaluator");
//Test taking one step.
vtkm::Vec<FieldType, 3> result = p + vec * stepSize;
ValidateIntegrator(constRK4, p, result, "constant vector RK4");
ValidateIntegrator(uniformRK4, p, result, "uniform RK4");
ValidateIntegrator(rectRK4, p, result, "rectilinear RK4");
}
}
}
}
}
}
void TestParticleWorklets()
{
typedef VTKM_DEFAULT_DEVICE_ADAPTER_TAG DeviceAdapter;
typedef vtkm::Float32 FieldType;
@ -126,11 +346,16 @@ void TestParticleAdvection()
}
vtkm::cont::DataSet ds = dataSetBuilder.Create(dims);
typedef vtkm::worklet::particleadvection::UniformGridEvaluate<FieldPortalConstType, FieldType>
typedef vtkm::worklet::particleadvection::UniformGridEvaluate<FieldPortalConstType,
FieldType,
DeviceAdapter>
RGEvalType;
typedef vtkm::worklet::particleadvection::RK4Integrator<RGEvalType, FieldType> RK4RGType;
RGEvalType eval(ds);
vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>> fieldArray;
fieldArray = vtkm::cont::make_ArrayHandle(field);
RGEvalType eval(ds.GetCoordinateSystem(), ds.GetCellSet(0), fieldArray);
RK4RGType rk4(eval, stepSize);
std::vector<vtkm::Vec<FieldType, 3>> pts;
@ -141,16 +366,19 @@ void TestParticleAdvection()
vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>> seeds;
seeds = vtkm::cont::make_ArrayHandle(pts);
vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>> fieldArray;
fieldArray = vtkm::cont::make_ArrayHandle(field);
vtkm::worklet::ParticleAdvection particleAdvection;
vtkm::worklet::ParticleAdvectionResult<FieldType> res;
res = particleAdvection.Run(rk4, seeds, fieldArray, 1000, DeviceAdapter());
res = particleAdvection.Run(rk4, seeds, 1000, DeviceAdapter());
VTKM_TEST_ASSERT(res.positions.GetNumberOfValues() == seeds.GetNumberOfValues(),
"Number of output particles does not match input.");
}
void TestParticleAdvection()
{
TestEvaluators();
TestParticleWorklets();
}
int UnitTestParticleAdvection(int, char* [])
{
return vtkm::cont::testing::Testing::Run(TestParticleAdvection);