Shuenhoy/vtk-m
1
0
Fork 0
mirror of https://gitlab.kitware.com/vtk/vtk-m synced 2024-09-19 18:45:43 +00:00
Code Issues 2 Actions Packages Projects Releases Wiki Activity

Fixes for streamlines.

Browse Source
This commit is contained in:
dpugmire 2019-12-05 15:46:41 -05:00
parent bdab4a56d2
commit 8054dfdafb
8 changed files with 242 additions and 147 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
vtkm/filter/Streamline.hxx
View File

@ -72,8 +72,8 @@ inline VTKM_CONT vtkm::cont::DataSet Streamline::DoExecute(
res = this->Worklet.Run(rk4, seedArray, this->NumberOfSteps); res = this->Worklet.Run(rk4, seedArray, this->NumberOfSteps);
vtkm::cont::DataSet outData; vtkm::cont::DataSet outData;
vtkm::cont::CoordinateSystem outputCoords("coordinates", res.positions); vtkm::cont::CoordinateSystem outputCoords("coordinates", res.Positions);
outData.SetCellSet(res.polyLines); outData.SetCellSet(res.PolyLines);
outData.AddCoordinateSystem(outputCoords); outData.AddCoordinateSystem(outputCoords);
return outData; return outData;

47
vtkm/worklet/ParticleAdvection.h
View File

@ -80,44 +80,24 @@ public:
struct StreamlineResult struct StreamlineResult
{ {
StreamlineResult() StreamlineResult()
: positions() : Particles()
, polyLines() , Positions()
, status() , PolyLines()
, stepsTaken()
, times()
{ {
} }
StreamlineResult(const vtkm::cont::ArrayHandle<vtkm::Vec3f>& pos, StreamlineResult(const vtkm::cont::ArrayHandle<vtkm::Particle>& part,
const vtkm::cont::CellSetExplicit<>& lines, const vtkm::cont::ArrayHandle<vtkm::Vec3f>& pos,
const vtkm::cont::ArrayHandle<vtkm::Id>& stat, const vtkm::cont::CellSetExplicit<>& lines)
const vtkm::cont::ArrayHandle<vtkm::Id>& steps) : Particles(part)
: positions(pos) , Positions(pos)
, polyLines(lines) , PolyLines(lines)
, status(stat)
, stepsTaken(steps)
{ {
} }
StreamlineResult(const vtkm::cont::ArrayHandle<vtkm::Vec3f>& pos, vtkm::cont::ArrayHandle<vtkm::Particle> Particles;
const vtkm::cont::CellSetExplicit<>& lines, vtkm::cont::ArrayHandle<vtkm::Vec3f> Positions;
const vtkm::cont::ArrayHandle<vtkm::Id>& stat, vtkm::cont::CellSetExplicit<> PolyLines;
const vtkm::cont::ArrayHandle<vtkm::Id>& steps,
const vtkm::cont::ArrayHandle<vtkm::FloatDefault>& timeArray)
: positions(pos)
, polyLines(lines)
, status(stat)
, stepsTaken(steps)
, times(timeArray)
{
}
vtkm::cont::ArrayHandle<vtkm::Vec3f> positions;
vtkm::cont::CellSetExplicit<> polyLines;
vtkm::cont::ArrayHandle<vtkm::Id> status;
vtkm::cont::ArrayHandle<vtkm::Id> stepsTaken;
vtkm::cont::ArrayHandle<vtkm::FloatDefault> times;
}; };
class Streamline class Streamline
@ -137,8 +117,7 @@ public:
worklet.Run(it, particles, MaxSteps, positions, polyLines); worklet.Run(it, particles, MaxSteps, positions, polyLines);
std::cout << "FIX ME " << __FILE__ << " " << __LINE__ << std::endl; return StreamlineResult(particles, positions, polyLines);
return StreamlineResult();
} }
}; };
} }

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

@ -63,8 +63,8 @@ public:
integralCurve.PreStepUpdate(idx); integralCurve.PreStepUpdate(idx);
do do
{ {
vtkm::Particle p = integralCurve.GetParticle(idx); //vtkm::Particle p = integralCurve.GetParticle(idx);
std::cout << idx << ": " << inpos << " #" << p.NumSteps << " " << p.Status << std::endl; //std::cout<<idx<<": "<<inpos<<" #"<<p.NumSteps<<" "<<p.Status<<std::endl;
vtkm::Vec3f outpos; vtkm::Vec3f outpos;
auto status = integrator->Step(inpos, time, outpos); auto status = integrator->Step(inpos, time, outpos);
if (status.CheckOk()) if (status.CheckOk())
@ -96,9 +96,8 @@ public:
//Mark if any steps taken //Mark if any steps taken
integralCurve.UpdateTookSteps(idx, tookAnySteps); integralCurve.UpdateTookSteps(idx, tookAnySteps);
particle = integralCurve.GetParticle(idx); //particle = integralCurve.GetParticle(idx);
std::cout << idx << ": " << inpos << " #" << particle.NumSteps << " " << particle.Status //std::cout<<idx<<": "<<inpos<<" #"<<particle.NumSteps<<" "<<particle.Status<<std::endl;
<< std::endl;
} }
}; };
@ -143,16 +142,34 @@ public:
namespace detail namespace detail
{ {
class Subtract : public vtkm::worklet::WorkletMapField class GetSteps : public vtkm::worklet::WorkletMapField
{ {
public: public:
VTKM_CONT VTKM_CONT
Subtract() {} GetSteps() {}
using ControlSignature = void(FieldOut, FieldIn, FieldIn); using ControlSignature = void(FieldIn, FieldOut);
using ExecutionSignature = void(_1, _2, _3); using ExecutionSignature = void(_1, _2);
VTKM_EXEC void operator()(vtkm::Id& res, const vtkm::Id& x, const vtkm::Id& y) const VTKM_EXEC void operator()(const vtkm::Particle& p, vtkm::Id& numSteps) const
{ {
res = x - y; numSteps = p.NumSteps;
}
};
class ComputeNumPoints : public vtkm::worklet::WorkletMapField
{
public:
VTKM_CONT
ComputeNumPoints() {}
using ControlSignature = void(FieldIn, FieldIn, FieldOut);
using ExecutionSignature = void(_1, _2, _3);
// Offset is number of points in streamline.
// 1 (inital point) + number of steps taken (p.NumSteps - initalNumSteps)
VTKM_EXEC void operator()(const vtkm::Particle& p,
const vtkm::Id& initialNumSteps,
vtkm::Id& diff) const
{
diff = 1 + p.NumSteps - initialNumSteps;
} }
}; };
} // namespace detail } // namespace detail
@ -174,12 +191,20 @@ public:
vtkm::worklet::particleadvection::ParticleAdvectWorklet>; vtkm::worklet::particleadvection::ParticleAdvectWorklet>;
using StreamlineType = vtkm::worklet::particleadvection::StateRecordingParticles; using StreamlineType = vtkm::worklet::particleadvection::StateRecordingParticles;
//NEED TO DEFINE THESE vtkm::cont::ArrayHandle<vtkm::Id> initialStepsTaken;
vtkm::cont::ArrayHandle<vtkm::Id> stepsTaken, initialStepsTaken, offsets;
vtkm::Id numSeeds = static_cast<vtkm::Id>(particles.GetNumberOfValues()); vtkm::Id numSeeds = static_cast<vtkm::Id>(particles.GetNumberOfValues());
vtkm::cont::ArrayHandleIndex idxArray(numSeeds); vtkm::cont::ArrayHandleIndex idxArray(numSeeds);
vtkm::worklet::DispatcherMapField<detail::GetSteps> getStepDispatcher{ (detail::GetSteps{}) };
getStepDispatcher.Invoke(particles, initialStepsTaken);
#ifdef VTKM_CUDA
// This worklet needs some extra space on CUDA.
vtkm::cont::cuda::ScopedCudaStackSize stack(16 * 1024);
(void)stack;
#endif // VTKM_CUDA
//Run streamline worklet //Run streamline worklet
StreamlineType streamlines(particles, MaxSteps); StreamlineType streamlines(particles, MaxSteps);
ParticleWorkletDispatchType particleWorkletDispatch; ParticleWorkletDispatchType particleWorkletDispatch;
@ -188,16 +213,15 @@ public:
//Get the positions //Get the positions
streamlines.GetCompactedHistory(positions); streamlines.GetCompactedHistory(positions);
vtkm::cont::printSummary_ArrayHandle(positions, std::cout);
//Create the cells //Create the cells
vtkm::cont::ArrayHandle<vtkm::Id> stepsTakenNow; vtkm::cont::ArrayHandle<vtkm::Id> numPoints;
stepsTakenNow.Allocate(numSeeds); vtkm::worklet::DispatcherMapField<detail::ComputeNumPoints> computeNumPointsDispatcher{ (
vtkm::worklet::DispatcherMapField<detail::Subtract> subtractDispatcher{ (detail::Subtract{}) }; detail::ComputeNumPoints{}) };
subtractDispatcher.Invoke(stepsTakenNow, stepsTaken, initialStepsTaken); computeNumPointsDispatcher.Invoke(particles, initialStepsTaken, numPoints);
vtkm::cont::ArrayHandle<vtkm::Id> cellIndex; vtkm::cont::ArrayHandle<vtkm::Id> cellIndex;
vtkm::Id connectivityLen = vtkm::cont::Algorithm::ScanExclusive(stepsTakenNow, cellIndex); vtkm::Id connectivityLen = vtkm::cont::Algorithm::ScanExclusive(numPoints, cellIndex);
vtkm::cont::ArrayHandleCounting<vtkm::Id> connCount(0, 1, connectivityLen); vtkm::cont::ArrayHandleCounting<vtkm::Id> connCount(0, 1, connectivityLen);
vtkm::cont::ArrayHandle<vtkm::Id> connectivity; vtkm::cont::ArrayHandle<vtkm::Id> connectivity;
vtkm::cont::ArrayCopy(connCount, connectivity); vtkm::cont::ArrayCopy(connCount, connectivity);
@ -207,34 +231,9 @@ public:
vtkm::cont::make_ArrayHandleConstant<vtkm::UInt8>(vtkm::CELL_SHAPE_POLY_LINE, numSeeds); vtkm::cont::make_ArrayHandleConstant<vtkm::UInt8>(vtkm::CELL_SHAPE_POLY_LINE, numSeeds);
vtkm::cont::ArrayCopy(polyLineShape, cellTypes); vtkm::cont::ArrayCopy(polyLineShape, cellTypes);
auto numIndices = vtkm::cont::make_ArrayHandleCast(stepsTakenNow, vtkm::IdComponent()); auto numIndices = vtkm::cont::make_ArrayHandleCast(numPoints, vtkm::IdComponent());
//auto offsets = vtkm::cont::ConvertNumIndicesToOffsets(numIndices); auto offsets = vtkm::cont::ConvertNumIndicesToOffsets(numIndices);
polyLines.Fill(positions.GetNumberOfValues(), cellTypes, connectivity, offsets); polyLines.Fill(positions.GetNumberOfValues(), cellTypes, connectivity, offsets);
/*
using ParticleAdvectWorkletType = vtkm::worklet::particleadvection::StreamlineWorkletAOS;
using ParticleWorkletDispatchType =
typename vtkm::worklet::DispatcherMapField<ParticleAdvectWorkletType>;
vtkm::Id numSeeds = static_cast<vtkm::Id>(particles.GetNumberOfValues());
//Create and invoke the particle advection.
vtkm::cont::ArrayHandleConstant<vtkm::Id> maxSteps(MaxSteps, numSeeds);
vtkm::cont::ArrayHandleIndex idxArray(numSeeds);
#ifdef VTKM_CUDA
// This worklet needs some extra space on CUDA.
vtkm::cont::cuda::ScopedCudaStackSize stack(16 * 1024);
(void)stack;
#endif // VTKM_CUDA
//Invoke particle advection worklet
ParticleWorkletDispatchType particleWorkletDispatch;
particleWorkletDispatch.Invoke(idxArray, it, particles, maxSteps);
*/
} }
}; };
} }

31
vtkm/worklet/particleadvection/Particles.h
View File

@ -139,6 +139,7 @@ public:
: ParticleExecutionObject<Device>() : ParticleExecutionObject<Device>()
, History() , History()
, Length(0) , Length(0)
, StepCount()
, ValidPoint() , ValidPoint()
{ {
} }
@ -147,6 +148,7 @@ public:
StateRecordingParticleExecutionObject(vtkm::cont::ArrayHandle<vtkm::Particle> pArray, StateRecordingParticleExecutionObject(vtkm::cont::ArrayHandle<vtkm::Particle> pArray,
vtkm::cont::ArrayHandle<vtkm::Vec3f> historyArray, vtkm::cont::ArrayHandle<vtkm::Vec3f> historyArray,
vtkm::cont::ArrayHandle<vtkm::Id> validPointArray, vtkm::cont::ArrayHandle<vtkm::Id> validPointArray,
vtkm::cont::ArrayHandle<vtkm::Id> stepCountArray,
vtkm::Id maxSteps) vtkm::Id maxSteps)
: ParticleExecutionObject<Device>(pArray, maxSteps) : ParticleExecutionObject<Device>(pArray, maxSteps)
, Length(maxSteps + 1) , Length(maxSteps + 1)
@ -154,6 +156,7 @@ public:
vtkm::Id numPos = pArray.GetNumberOfValues(); vtkm::Id numPos = pArray.GetNumberOfValues();
History = historyArray.PrepareForOutput(numPos * Length, Device()); History = historyArray.PrepareForOutput(numPos * Length, Device());
ValidPoint = validPointArray.PrepareForInPlace(Device()); ValidPoint = validPointArray.PrepareForInPlace(Device());
StepCount = stepCountArray.PrepareForInPlace(Device());
} }
VTKM_EXEC VTKM_EXEC
@ -163,9 +166,10 @@ public:
if (p.NumSteps == 0) if (p.NumSteps == 0)
{ {
vtkm::Id loc = idx * Length; vtkm::Id loc = idx * Length;
std::cout << "PreStepUpdate " << idx << ": loc= " << loc << " " << p.Pos << std::endl; //std::cout<<"PreStepUpdate "<<idx<<": loc= "<<loc<<" "<<p.Pos<<std::endl;
this->History.Set(loc, p.Pos); this->History.Set(loc, p.Pos);
this->ValidPoint.Set(loc, 1); this->ValidPoint.Set(loc, 1);
this->StepCount.Set(idx, 1);
} }
} }
@ -174,12 +178,16 @@ public:
{ {
this->ParticleExecutionObject<Device>::StepUpdate(idx, time, pt); this->ParticleExecutionObject<Device>::StepUpdate(idx, time, pt);
vtkm::Particle p = this->ParticleExecutionObject<Device>::GetParticle(idx); //vtkm::Particle p = this->ParticleExecutionObject<Device>::GetParticle(idx);
vtkm::Id loc = idx * Length + p.NumSteps; //local step count.
std::cout << "StepUpdate " << idx << ": loc= " << loc << " " << pt << std::endl; vtkm::Id stepCount = this->StepCount.Get(idx);
vtkm::Id loc = idx * Length + stepCount;
//std::cout<<"StepUpdate "<<idx<<": loc= "<<loc<<" "<<pt<<std::endl;
this->History.Set(loc, pt); this->History.Set(loc, pt);
this->ValidPoint.Set(loc, 1); this->ValidPoint.Set(loc, 1);
this->StepCount.Set(idx, stepCount + 1);
} }
protected: protected:
@ -190,6 +198,7 @@ protected:
HistoryPortal History; HistoryPortal History;
vtkm::Id Length; vtkm::Id Length;
IdPortal StepCount;
IdPortal ValidPoint; IdPortal ValidPoint;
}; };
@ -201,7 +210,11 @@ public:
PrepareForExecution(Device) const PrepareForExecution(Device) const
{ {
return vtkm::worklet::particleadvection::StateRecordingParticleExecutionObject<Device>( return vtkm::worklet::particleadvection::StateRecordingParticleExecutionObject<Device>(
ParticleArray, HistoryArray, ValidPointArray, MaxSteps); this->ParticleArray,
this->HistoryArray,
this->ValidPointArray,
this->StepCountArray,
this->MaxSteps);
} }
VTKM_CONT VTKM_CONT
StateRecordingParticles(vtkm::cont::ArrayHandle<vtkm::Particle>& pArray, const vtkm::Id& maxSteps) StateRecordingParticles(vtkm::cont::ArrayHandle<vtkm::Particle>& pArray, const vtkm::Id& maxSteps)
@ -209,9 +222,14 @@ public:
, ParticleArray(pArray) , ParticleArray(pArray)
{ {
vtkm::Id numParticles = static_cast<vtkm::Id>(pArray.GetNumberOfValues()); vtkm::Id numParticles = static_cast<vtkm::Id>(pArray.GetNumberOfValues());
std::cout << "ctor size: " << this->MaxSteps + 1 << " " << numParticles << std::endl;
//Create ValidPointArray initialized to zero.
vtkm::cont::ArrayHandleConstant<vtkm::Id> tmp(0, (this->MaxSteps + 1) * numParticles); vtkm::cont::ArrayHandleConstant<vtkm::Id> tmp(0, (this->MaxSteps + 1) * numParticles);
vtkm::cont::ArrayCopy(tmp, this->ValidPointArray); vtkm::cont::ArrayCopy(tmp, this->ValidPointArray);
//Create StepCountArray initialized to zero.
vtkm::cont::ArrayHandleConstant<vtkm::Id> tmp2(0, numParticles);
vtkm::cont::ArrayCopy(tmp2, this->StepCountArray);
} }
VTKM_CONT VTKM_CONT
@ -236,6 +254,7 @@ protected:
vtkm::cont::ArrayHandle<vtkm::Vec3f> HistoryArray; vtkm::cont::ArrayHandle<vtkm::Vec3f> HistoryArray;
vtkm::Id MaxSteps; vtkm::Id MaxSteps;
vtkm::cont::ArrayHandle<vtkm::Particle> ParticleArray; vtkm::cont::ArrayHandle<vtkm::Particle> ParticleArray;
vtkm::cont::ArrayHandle<vtkm::Id> StepCountArray;
vtkm::cont::ArrayHandle<vtkm::Id> ValidPointArray; vtkm::cont::ArrayHandle<vtkm::Id> ValidPointArray;
struct IsOne struct IsOne

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

@ -11,6 +11,7 @@
#ifndef vtk_m_worklet_particleadvection_TemporalGridEvaluators_h #ifndef vtk_m_worklet_particleadvection_TemporalGridEvaluators_h
#define vtk_m_worklet_particleadvection_TemporalGridEvaluators_h #define vtk_m_worklet_particleadvection_TemporalGridEvaluators_h
#include <vtkm/worklet/particleadvection/GridEvaluatorStatus.h>
#include <vtkm/worklet/particleadvection/GridEvaluators.h> #include <vtkm/worklet/particleadvection/GridEvaluators.h>
namespace vtkm namespace vtkm
@ -68,10 +69,12 @@ public:
} }
template <typename Point> template <typename Point>
VTKM_EXEC EvaluatorStatus Evaluate(const Point& pos, vtkm::FloatDefault time, Point& out) const VTKM_EXEC GridEvaluatorStatus Evaluate(const Point& pos,
vtkm::FloatDefault time,
Point& out) const
{ {
// Validate time is in bounds for the current two slices. // Validate time is in bounds for the current two slices.
EvaluatorStatus status; GridEvaluatorStatus status;
if (!(time >= TimeOne && time <= TimeTwo)) if (!(time >= TimeOne && time <= TimeTwo))
{ {
@ -91,6 +94,8 @@ public:
// LERP between the two values of calculated fields to obtain the new value // LERP between the two values of calculated fields to obtain the new value
vtkm::FloatDefault proportion = (time - this->TimeOne) / this->TimeDiff; vtkm::FloatDefault proportion = (time - this->TimeOne) / this->TimeDiff;
out = vtkm::Lerp(one, two, proportion); out = vtkm::Lerp(one, two, proportion);
status.SetOk();
return status; return status;
} }

2
vtkm/worklet/testing/CMakeLists.txt
View File

@ -62,7 +62,7 @@ set(unit_tests
UnitTestStreamLineUniformGrid.cxx UnitTestStreamLineUniformGrid.cxx
UnitTestStreamSurface.cxx UnitTestStreamSurface.cxx
UnitTestSurfaceNormals.cxx UnitTestSurfaceNormals.cxx
# UnitTestTemporalAdvection.cxx UnitTestTemporalAdvection.cxx
UnitTestTetrahedralize.cxx UnitTestTetrahedralize.cxx
UnitTestThreshold.cxx UnitTestThreshold.cxx
UnitTestThresholdPoints.cxx UnitTestThresholdPoints.cxx

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

@ -534,23 +534,33 @@ void ValidateParticleAdvectionResult(const vtkm::worklet::ParticleAdvectionResul
VTKM_TEST_ASSERT(res.Particles.GetNumberOfValues() == nSeeds, VTKM_TEST_ASSERT(res.Particles.GetNumberOfValues() == nSeeds,
"Number of output particles does not match input."); "Number of output particles does not match input.");
for (vtkm::Id i = 0; i < nSeeds; i++) for (vtkm::Id i = 0; i < nSeeds; i++)
{
VTKM_TEST_ASSERT(res.Particles.GetPortalConstControl().Get(i).NumSteps <= maxSteps, VTKM_TEST_ASSERT(res.Particles.GetPortalConstControl().Get(i).NumSteps <= maxSteps,
"Too many steps taken in particle advection"); "Too many steps taken in particle advection");
VTKM_TEST_ASSERT(res.Particles.GetPortalConstControl().Get(i).Status.CheckOk(),
"Bad status in particle advection");
}
} }
void ValidateStreamlineResult(const vtkm::worklet::StreamlineResult& res, void ValidateStreamlineResult(const vtkm::worklet::StreamlineResult& res,
vtkm::Id nSeeds, vtkm::Id nSeeds,
vtkm::Id maxSteps) vtkm::Id maxSteps)
{ {
VTKM_TEST_ASSERT(res.polyLines.GetNumberOfCells() == nSeeds, VTKM_TEST_ASSERT(res.PolyLines.GetNumberOfCells() == nSeeds,
"Number of output streamlines does not match input."); "Number of output streamlines does not match input.");
for (vtkm::Id i = 0; i < nSeeds; i++) for (vtkm::Id i = 0; i < nSeeds; i++)
VTKM_TEST_ASSERT(res.stepsTaken.GetPortalConstControl().Get(i) <= maxSteps, {
VTKM_TEST_ASSERT(res.Particles.GetPortalConstControl().Get(i).NumSteps <= maxSteps,
"Too many steps taken in streamline"); "Too many steps taken in streamline");
VTKM_TEST_ASSERT(res.Particles.GetPortalConstControl().Get(i).Status.CheckOk(),
"Bad status in streamline");
}
VTKM_TEST_ASSERT(res.Particles.GetNumberOfValues() == nSeeds,
"Number of output particles does not match input.");
} }
void TestParticleWorklets() void TestParticleWorkletsWithDataSetTypes()
{ {
using FieldHandle = vtkm::cont::ArrayHandle<vtkm::Vec3f>; using FieldHandle = vtkm::cont::ArrayHandle<vtkm::Vec3f>;
using GridEvalType = vtkm::worklet::particleadvection::GridEvaluator<FieldHandle>; using GridEvalType = vtkm::worklet::particleadvection::GridEvaluator<FieldHandle>;
@ -641,9 +651,7 @@ void TestParticleWorklets()
auto seeds = vtkm::cont::make_ArrayHandle(pts2, vtkm::CopyFlag::On); auto seeds = vtkm::cont::make_ArrayHandle(pts2, vtkm::CopyFlag::On);
res = s.Run(rk4, seeds, maxSteps); res = s.Run(rk4, seeds, maxSteps);
} }
ValidateStreamlineResult(res, nSeeds, maxSteps);
std::cout << "FIX ME.. Add some testing " << __FILE__ << " " << __LINE__ << std::endl;
//ValidateStreamlineResult(res, nSeeds, maxSteps);
} }
} }
} }
@ -688,7 +696,7 @@ void TestParticleStatus()
VTKM_TEST_ASSERT(tookStep1 == false, "Particle took a step when it should not have."); VTKM_TEST_ASSERT(tookStep1 == false, "Particle took a step when it should not have.");
} }
void TestStreamline() void TestWorkletsBasic()
{ {
using FieldHandle = vtkm::cont::ArrayHandle<vtkm::Vec3f>; using FieldHandle = vtkm::cont::ArrayHandle<vtkm::Vec3f>;
using GridEvalType = vtkm::worklet::particleadvection::GridEvaluator<FieldHandle>; using GridEvalType = vtkm::worklet::particleadvection::GridEvaluator<FieldHandle>;
@ -699,14 +707,10 @@ void TestStreamline()
vtkm::Id nElements = dims[0] * dims[1] * dims[2] * 3; vtkm::Id nElements = dims[0] * dims[1] * dims[2] * 3;
std::vector<vtkm::Vec3f> field; std::vector<vtkm::Vec3f> field;
vtkm::Vec3f vecDir(1, 0, 0);
for (vtkm::Id i = 0; i < nElements; i++) for (vtkm::Id i = 0; i < nElements; i++)
{ field.push_back(vtkm::Normal(vecDir));
vtkm::FloatDefault x = vecData[i];
vtkm::FloatDefault y = vecData[++i];
vtkm::FloatDefault z = vecData[++i];
vtkm::Vec3f vec(x, y, z);
field.push_back(vtkm::Normal(vec));
}
vtkm::cont::ArrayHandle<vtkm::Vec3f> fieldArray; vtkm::cont::ArrayHandle<vtkm::Vec3f> fieldArray;
fieldArray = vtkm::cont::make_ArrayHandle(field); fieldArray = vtkm::cont::make_ArrayHandle(field);
@ -716,24 +720,112 @@ void TestStreamline()
GridEvalType eval(ds.GetCoordinateSystem(), ds.GetCellSet(), fieldArray); GridEvalType eval(ds.GetCoordinateSystem(), ds.GetCellSet(), fieldArray);
RK4Type rk4(eval, stepSize); RK4Type rk4(eval, stepSize);
std::vector<vtkm::Particle> pts; vtkm::Id maxSteps = 83;
pts.push_back(vtkm::Particle(vtkm::Vec3f(.5, .5, .5), 0)); std::vector<std::string> workletTypes = { "particleAdvection", "streamline" };
auto seedsArray = vtkm::cont::make_ArrayHandle(pts, vtkm::CopyFlag::On); vtkm::FloatDefault endT = stepSize * maxSteps;
vtkm::Id maxSteps = 100;
vtkm::worklet::Streamline s; for (auto w : workletTypes)
vtkm::worklet::StreamlineResult res; {
std::vector<vtkm::Particle> particles;
std::vector<vtkm::Vec3f> pts, samplePts, endPts;
vtkm::FloatDefault X = static_cast<vtkm::FloatDefault>(.1);
vtkm::FloatDefault Y = static_cast<vtkm::FloatDefault>(.1);
vtkm::FloatDefault Z = static_cast<vtkm::FloatDefault>(.1);
res = s.Run(rk4, seedsArray, maxSteps); for (int i = 0; i < 8; i++)
std::cout << "FIX ME.. Add some testing " << __FILE__ << " " << __LINE__ << std::endl; {
pts.push_back(vtkm::Vec3f(X, Y, Z));
Y += static_cast<vtkm::FloatDefault>(.1);
}
vtkm::Id id = 0;
for (std::size_t i = 0; i < pts.size(); i++, id++)
{
vtkm::Vec3f p = pts[i];
particles.push_back(vtkm::Particle(p, id));
samplePts.push_back(p);
for (vtkm::Id j = 0; j < maxSteps; j++)
{
p = p + vecDir * stepSize;
samplePts.push_back(p);
}
endPts.push_back(p);
}
auto seedsArray = vtkm::cont::make_ArrayHandle(particles, vtkm::CopyFlag::On);
if (w == "particleAdvection")
{
vtkm::worklet::ParticleAdvection pa;
vtkm::worklet::ParticleAdvectionResult res;
res = pa.Run(rk4, seedsArray, maxSteps);
vtkm::Id numRequiredPoints = static_cast<vtkm::Id>(endPts.size());
VTKM_TEST_ASSERT(res.Particles.GetNumberOfValues() == numRequiredPoints,
"Wrong number of points in particle advection result.");
auto portal = res.Particles.GetPortalConstControl();
for (vtkm::Id i = 0; i < res.Particles.GetNumberOfValues(); i++)
{
VTKM_TEST_ASSERT(portal.Get(i).Pos == endPts[i], "Particle advection point is wrong");
VTKM_TEST_ASSERT(portal.Get(i).NumSteps == maxSteps,
"Particle advection NumSteps is wrong");
VTKM_TEST_ASSERT(vtkm::Abs(portal.Get(i).Time - endT) < stepSize / 100,
"Particle advection Time is wrong");
VTKM_TEST_ASSERT(portal.Get(i).Status.CheckOk(), "Particle advection Status is wrong");
VTKM_TEST_ASSERT(portal.Get(i).Status.CheckTerminate(),
"Particle advection particle did not terminate");
}
}
else if (w == "streamline")
{
vtkm::worklet::Streamline s;
vtkm::worklet::StreamlineResult res;
res = s.Run(rk4, seedsArray, maxSteps);
vtkm::Id numRequiredPoints = static_cast<vtkm::Id>(samplePts.size());
VTKM_TEST_ASSERT(res.Positions.GetNumberOfValues() == numRequiredPoints,
"Wrong number of points in streamline result.");
//Make sure all the points match.
auto parPortal = res.Particles.GetPortalConstControl();
for (vtkm::Id i = 0; i < res.Particles.GetNumberOfValues(); i++)
{
VTKM_TEST_ASSERT(parPortal.Get(i).Pos == endPts[i], "Streamline end point is wrong");
VTKM_TEST_ASSERT(parPortal.Get(i).NumSteps == maxSteps, "Streamline NumSteps is wrong");
VTKM_TEST_ASSERT(vtkm::Abs(parPortal.Get(i).Time - endT) < stepSize / 100,
"Streamline Time is wrong");
VTKM_TEST_ASSERT(parPortal.Get(i).Status.CheckOk(), "Streamline Status is wrong");
VTKM_TEST_ASSERT(parPortal.Get(i).Status.CheckTerminate(),
"Streamline particle did not terminate");
}
auto posPortal = res.Positions.GetPortalConstControl();
for (vtkm::Id i = 0; i < res.Positions.GetNumberOfValues(); i++)
VTKM_TEST_ASSERT(posPortal.Get(i) == samplePts[i], "Streamline points do not match");
vtkm::Id numCells = res.PolyLines.GetNumberOfCells();
VTKM_TEST_ASSERT(numCells == static_cast<vtkm::Id>(pts.size()),
"Wrong number of polylines in streamline");
for (vtkm::Id i = 0; i < numCells; i++)
{
VTKM_TEST_ASSERT(res.PolyLines.GetCellShape(i) == vtkm::CELL_SHAPE_POLY_LINE,
"Wrong cell type in streamline.");
VTKM_TEST_ASSERT(res.PolyLines.GetNumberOfPointsInCell(i) ==
static_cast<vtkm::Id>(maxSteps + 1),
"Wrong number of points in streamline cell");
}
}
}
} }
void TestParticleAdvection() void TestParticleAdvection()
{ {
TestStreamline();
TestEvaluators(); TestEvaluators();
TestParticleWorklets();
TestParticleStatus(); TestParticleStatus();
TestWorkletsBasic();
TestParticleWorkletsWithDataSetTypes();
} }
int UnitTestParticleAdvection(int argc, char* argv[]) int UnitTestParticleAdvection(int argc, char* argv[])

65
vtkm/worklet/testing/UnitTestTemporalAdvection.cxx
View File

@ -39,7 +39,6 @@ vtkm::cont::DataSet CreateUniformDataSet(const vtkm::Bounds& bounds, const vtkm:
return ds; return ds;
} }
template <typename ScalarType>
class TestEvaluatorWorklet : public vtkm::worklet::WorkletMapField class TestEvaluatorWorklet : public vtkm::worklet::WorkletMapField
{ {
public: public:
@ -51,29 +50,30 @@ public:
using ExecutionSignature = void(_1, _2, _3, _4); using ExecutionSignature = void(_1, _2, _3, _4);
template <typename EvaluatorType> template <typename EvaluatorType>
VTKM_EXEC void operator()(vtkm::Vec<ScalarType, 3>& pointIn, VTKM_EXEC void operator()(vtkm::Particle& pointIn,
const EvaluatorType& evaluator, const EvaluatorType& evaluator,
vtkm::worklet::particleadvection::EvaluatorStatus& status, vtkm::worklet::particleadvection::GridEvaluatorStatus& status,
vtkm::Vec<ScalarType, 3>& pointOut) const vtkm::Vec3f& pointOut) const
{ {
status = evaluator.Evaluate(pointIn, 0.5f, pointOut); status = evaluator.Evaluate(pointIn.Pos, 0.5f, pointOut);
} }
}; };
template <typename EvalType, typename ScalarType> template <typename EvalType>
void ValidateEvaluator(const EvalType& eval, void ValidateEvaluator(const EvalType& eval,
const vtkm::cont::ArrayHandle<vtkm::Vec<ScalarType, 3>>& pointIns, const vtkm::cont::ArrayHandle<vtkm::Particle>& pointIns,
const vtkm::cont::ArrayHandle<vtkm::Vec<ScalarType, 3>>& validity, const vtkm::cont::ArrayHandle<vtkm::Vec3f>& validity,
const std::string& msg) const std::string& msg)
{ {
using EvalTester = TestEvaluatorWorklet<ScalarType>; using EvalTester = TestEvaluatorWorklet;
using EvalTesterDispatcher = vtkm::worklet::DispatcherMapField<EvalTester>; using EvalTesterDispatcher = vtkm::worklet::DispatcherMapField<EvalTester>;
using Status = vtkm::worklet::particleadvection::EvaluatorStatus; using Status = vtkm::worklet::particleadvection::GridEvaluatorStatus;
EvalTester evalTester; EvalTester evalTester;
EvalTesterDispatcher evalTesterDispatcher(evalTester); EvalTesterDispatcher evalTesterDispatcher(evalTester);
vtkm::Id numPoints = pointIns.GetNumberOfValues(); vtkm::Id numPoints = pointIns.GetNumberOfValues();
vtkm::cont::ArrayHandle<Status> evalStatus; vtkm::cont::ArrayHandle<Status> evalStatus;
vtkm::cont::ArrayHandle<vtkm::Vec<ScalarType, 3>> evalResults; vtkm::cont::ArrayHandle<vtkm::Vec3f> evalResults;
evalTesterDispatcher.Invoke(pointIns, eval, evalStatus, evalResults); evalTesterDispatcher.Invoke(pointIns, eval, evalStatus, evalResults);
auto statusPortal = evalStatus.GetPortalConstControl(); auto statusPortal = evalStatus.GetPortalConstControl();
auto resultsPortal = evalResults.GetPortalConstControl(); auto resultsPortal = evalResults.GetPortalConstControl();
@ -81,8 +81,8 @@ void ValidateEvaluator(const EvalType& eval,
for (vtkm::Id index = 0; index < numPoints; index++) for (vtkm::Id index = 0; index < numPoints; index++)
{ {
Status status = statusPortal.Get(index); Status status = statusPortal.Get(index);
vtkm::Vec<ScalarType, 3> result = resultsPortal.Get(index); vtkm::Vec3f result = resultsPortal.Get(index);
vtkm::Vec<ScalarType, 3> expected = validityPortal.Get(index); vtkm::Vec3f expected = validityPortal.Get(index);
VTKM_TEST_ASSERT(status.CheckOk(), "Error in evaluator for " + msg); VTKM_TEST_ASSERT(status.CheckOk(), "Error in evaluator for " + msg);
VTKM_TEST_ASSERT(result == expected, "Error in evaluator result for " + msg); VTKM_TEST_ASSERT(result == expected, "Error in evaluator result for " + msg);
} }
@ -100,45 +100,45 @@ void CreateConstantVectorField(vtkm::Id num,
vtkm::cont::ArrayCopy(vecConst, vecField); vtkm::cont::ArrayCopy(vecConst, vecField);
} }
template <typename ScalarType> vtkm::Vec3f RandomPt(const vtkm::Bounds& bounds)
vtkm::Vec<ScalarType, 3> RandomPoint(const vtkm::Bounds& bounds)
{ {
ScalarType rx = static_cast<ScalarType>(rand()) / static_cast<ScalarType>(RAND_MAX); vtkm::FloatDefault rx =
ScalarType ry = static_cast<ScalarType>(rand()) / static_cast<ScalarType>(RAND_MAX); static_cast<vtkm::FloatDefault>(rand()) / static_cast<vtkm::FloatDefault>(RAND_MAX);
ScalarType rz = static_cast<ScalarType>(rand()) / static_cast<ScalarType>(RAND_MAX); vtkm::FloatDefault ry =
static_cast<vtkm::FloatDefault>(rand()) / static_cast<vtkm::FloatDefault>(RAND_MAX);
vtkm::FloatDefault rz =
static_cast<vtkm::FloatDefault>(rand()) / static_cast<vtkm::FloatDefault>(RAND_MAX);
vtkm::Vec<ScalarType, 3> p; vtkm::Vec3f p;
p[0] = static_cast<ScalarType>(bounds.X.Min + rx * bounds.X.Length()); p[0] = static_cast<vtkm::FloatDefault>(bounds.X.Min + rx * bounds.X.Length());
p[1] = static_cast<ScalarType>(bounds.Y.Min + ry * bounds.Y.Length()); p[1] = static_cast<vtkm::FloatDefault>(bounds.Y.Min + ry * bounds.Y.Length());
p[2] = static_cast<ScalarType>(bounds.Z.Min + rz * bounds.Z.Length()); p[2] = static_cast<vtkm::FloatDefault>(bounds.Z.Min + rz * bounds.Z.Length());
return p; return p;
} }
template <typename ScalarType>
void GeneratePoints(const vtkm::Id numOfEntries, void GeneratePoints(const vtkm::Id numOfEntries,
const vtkm::Bounds& bounds, const vtkm::Bounds& bounds,
vtkm::cont::ArrayHandle<vtkm::Vec<ScalarType, 3>>& pointIns) vtkm::cont::ArrayHandle<vtkm::Particle>& pointIns)
{ {
pointIns.Allocate(numOfEntries); pointIns.Allocate(numOfEntries);
auto writePortal = pointIns.GetPortalControl(); auto writePortal = pointIns.GetPortalControl();
for (vtkm::Id index = 0; index < numOfEntries; index++) for (vtkm::Id index = 0; index < numOfEntries; index++)
{ {
vtkm::Vec<ScalarType, 3> value = RandomPoint<ScalarType>(bounds); vtkm::Particle particle(RandomPt(bounds), index);
writePortal.Set(index, value); writePortal.Set(index, particle);
} }
} }
template <typename ScalarType>
void GenerateValidity(const vtkm::Id numOfEntries, void GenerateValidity(const vtkm::Id numOfEntries,
vtkm::cont::ArrayHandle<vtkm::Vec<ScalarType, 3>>& validity, vtkm::cont::ArrayHandle<vtkm::Vec3f>& validity,
const vtkm::Vec<ScalarType, 3>& vecOne, const vtkm::Vec3f& vecOne,
const vtkm::Vec<ScalarType, 3>& vecTwo) const vtkm::Vec3f& vecTwo)
{ {
validity.Allocate(numOfEntries); validity.Allocate(numOfEntries);
auto writePortal = validity.GetPortalControl(); auto writePortal = validity.GetPortalControl();
for (vtkm::Id index = 0; index < numOfEntries; index++) for (vtkm::Id index = 0; index < numOfEntries; index++)
{ {
vtkm::Vec<ScalarType, 3> value = 0.5f * vecOne + (1.0f - 0.5f) * vecTwo; vtkm::Vec3f value = 0.5f * vecOne + (1.0f - 0.5f) * vecTwo;
writePortal.Set(index, value); writePortal.Set(index, value);
} }
} }
@ -170,7 +170,8 @@ void TestTemporalEvaluators()
// Test data : populate with meaningful values // Test data : populate with meaningful values
vtkm::Id numValues = 10; vtkm::Id numValues = 10;
vtkm::cont::ArrayHandle<PointType> pointIns, validity; vtkm::cont::ArrayHandle<vtkm::Particle> pointIns;
vtkm::cont::ArrayHandle<vtkm::Vec3f> validity;
GeneratePoints(numValues, bounds, pointIns); GeneratePoints(numValues, bounds, pointIns);
GenerateValidity(numValues, validity, X, Z); GenerateValidity(numValues, validity, X, Z);