Shuenhoy/vtk-m
1
0
Fork 0
mirror of https://gitlab.kitware.com/vtk/vtk-m synced 2024-09-16 17:22:55 +00:00
Code Issues 2 Actions Packages Projects Releases Wiki Activity

Merge topic 'pa_storage'

Browse Source 
e4923823b Add comment for extra casts for compiler warnings.
98d4df79a compiler warnings.
a11bc1769 compiler warnings.
5646527b6 cleanup for LCS filter.
29c21cdbc Code cleanup. Update Lagranian filter to use vtkm::Particle
136aba429 Cleanup for streamsurface.
85fc73fc4 Updates for streamline/pathline filter.
33ba234d5 code cleanup
...

Acked-by: Kitware Robot <kwrobot@kitware.com>
Merge-request: !1854
This commit is contained in:
Dave Pugmire 2019-12-16 02:48:17 +00:00 committed by Kitware Robot
commit dd9dc19410
25 changed files with 1045 additions and 793 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

17
vtkm/Bitset.h
View File

@ -22,18 +22,27 @@ namespace vtkm
/// \brief A bitmap to serve different needs.
/// Ex. Editing particular bits in a byte(s), checkint if particular bit values
/// are present or not. Once Cuda supports std::bitset, we should use the
/// standard one if possible
/// standard one if possible. Additional cast in logical operations are required
/// to avoid compiler warnings when using 16 or 8 bit MaskType.
template <typename MaskType>
struct Bitset
{
VTKM_EXEC_CONT void set(vtkm::Id bitIndex)
{
this->Mask = this->Mask | (static_cast<MaskType>(1) << bitIndex);
this->Mask = static_cast<MaskType>(this->Mask | (static_cast<MaskType>(1) << bitIndex));
}
VTKM_EXEC_CONT void set(vtkm::Id bitIndex, bool val)
{
if (val)
this->set(bitIndex);
else
this->reset(bitIndex);
}
VTKM_EXEC_CONT void reset(vtkm::Id bitIndex)
{
this->Mask = this->Mask & ~(static_cast<MaskType>(1) << bitIndex);
this->Mask = static_cast<MaskType>(this->Mask & ~(static_cast<MaskType>(1) << bitIndex));
}
VTKM_EXEC_CONT void toggle(vtkm::Id bitIndex)
@ -41,7 +50,7 @@ struct Bitset
this->Mask = this->Mask ^ (static_cast<MaskType>(0) << bitIndex);
}
VTKM_EXEC_CONT bool test(vtkm::Id bitIndex)
VTKM_EXEC_CONT bool test(vtkm::Id bitIndex) const
{
return ((this->Mask & (static_cast<MaskType>(1) << bitIndex)) != 0);
}

1
vtkm/CMakeLists.txt
View File

@ -36,6 +36,7 @@ set(headers
Matrix.h
NewtonsMethod.h
Pair.h
Particle.h
Range.h
RangeId.h
RangeId2.h

110
vtkm/Particle.h Normal file
View File

@ -0,0 +1,110 @@
//============================================================================
// 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_Particle_h
#define vtk_m_Particle_h
#include <vtkm/Bitset.h>
namespace vtkm
{
//Bit field describing the status:
class ParticleStatus : public vtkm::Bitset<vtkm::UInt8>
{
public:
VTKM_EXEC_CONT ParticleStatus()
{
this->SetOk();
this->ClearTerminate();
}
VTKM_EXEC_CONT void SetOk() { this->set(this->SUCCESS_BIT); }
VTKM_EXEC_CONT bool CheckOk() const { return this->test(this->SUCCESS_BIT); }
VTKM_EXEC_CONT void SetFail() { this->reset(this->SUCCESS_BIT); }
VTKM_EXEC_CONT bool CheckFail() const { return !this->test(this->SUCCESS_BIT); }
VTKM_EXEC_CONT void SetTerminate() { this->set(this->TERMINATE_BIT); }
VTKM_EXEC_CONT void ClearTerminate() { this->reset(this->TERMINATE_BIT); }
VTKM_EXEC_CONT bool CheckTerminate() const { return this->test(this->TERMINATE_BIT); }
VTKM_EXEC_CONT void SetSpatialBounds() { this->set(this->SPATIAL_BOUNDS_BIT); }
VTKM_EXEC_CONT void ClearSpatialBounds() { this->reset(this->SPATIAL_BOUNDS_BIT); }
VTKM_EXEC_CONT bool CheckSpatialBounds() const { return this->test(this->SPATIAL_BOUNDS_BIT); }
VTKM_EXEC_CONT void SetTemporalBounds() { this->set(this->TEMPORAL_BOUNDS_BIT); }
VTKM_EXEC_CONT void ClearTemporalBounds() { this->reset(this->TEMPORAL_BOUNDS_BIT); }
VTKM_EXEC_CONT bool CheckTemporalBounds() const { return this->test(this->TEMPORAL_BOUNDS_BIT); }
VTKM_EXEC_CONT void SetTookAnySteps() { this->set(this->TOOK_ANY_STEPS_BIT); }
VTKM_EXEC_CONT void ClearTookAnySteps() { this->reset(this->TOOK_ANY_STEPS_BIT); }
VTKM_EXEC_CONT bool CheckTookAnySteps() const { return this->test(this->TOOK_ANY_STEPS_BIT); }
private:
static constexpr vtkm::Id SUCCESS_BIT = 0;
static constexpr vtkm::Id TERMINATE_BIT = 1;
static constexpr vtkm::Id SPATIAL_BOUNDS_BIT = 2;
static constexpr vtkm::Id TEMPORAL_BOUNDS_BIT = 3;
static constexpr vtkm::Id TOOK_ANY_STEPS_BIT = 4;
};
inline VTKM_CONT std::ostream& operator<<(std::ostream& s, const vtkm::ParticleStatus& status)
{
s << "[" << status.CheckOk() << " " << status.CheckTerminate() << " "
<< status.CheckSpatialBounds() << " " << status.CheckTemporalBounds() << "]";
return s;
}
class Particle
{
public:
VTKM_EXEC_CONT
Particle()
: Pos()
, ID(-1)
, NumSteps(0)
, Status()
, Time(0)
{
}
VTKM_EXEC_CONT
Particle(const vtkm::Particle& p)
: Pos(p.Pos)
, ID(p.ID)
, NumSteps(p.NumSteps)
, Status(p.Status)
, Time(p.Time)
{
}
VTKM_EXEC_CONT
Particle(const vtkm::Vec3f& p,
const vtkm::Id& id,
const vtkm::Id& numSteps = 0,
const vtkm::ParticleStatus& status = vtkm::ParticleStatus(),
const vtkm::FloatDefault& time = 0)
: Pos(p)
, ID(id)
, NumSteps(numSteps)
, Status(status)
, Time(time)
{
}
vtkm::Vec3f Pos;
vtkm::Id ID;
vtkm::Id NumSteps;
vtkm::ParticleStatus Status;
vtkm::FloatDefault Time;
};
}
#endif // vtk_m_Particle_h

69
vtkm/filter/Lagrangian.hxx
View File

@ -32,8 +32,8 @@
#include <string.h>
static vtkm::Id cycle = 0;
static vtkm::cont::ArrayHandle<vtkm::Vec3f_64> BasisParticles;
static vtkm::cont::ArrayHandle<vtkm::Vec3f_64> BasisParticlesOriginal;
static vtkm::cont::ArrayHandle<vtkm::Particle> BasisParticles;
static vtkm::cont::ArrayHandle<vtkm::Particle> BasisParticlesOriginal;
static vtkm::cont::ArrayHandle<vtkm::Id> BasisParticlesValidity;
namespace
@ -41,8 +41,8 @@ namespace
class ValidityCheck : public vtkm::worklet::WorkletMapField
{
public:
using ControlSignature = void(FieldIn end_point, FieldIn steps, FieldInOut output);
using ExecutionSignature = void(_1, _2, _3);
using ControlSignature = void(FieldIn end_point, FieldInOut output);
using ExecutionSignature = void(_1, _2);
using InputDomain = _1;
ValidityCheck(vtkm::Bounds b)
@ -50,16 +50,15 @@ public:
{
}
template <typename PosType, typename StepType, typename ValidityType>
VTKM_EXEC void operator()(const PosType& end_point,
const StepType& steps,
ValidityType& res) const
template <typename ValidityType>
VTKM_EXEC void operator()(const vtkm::Particle& end_point, ValidityType& res) const
{
vtkm::Id steps = end_point.NumSteps;
if (steps > 0 && res == 1)
{
if (end_point[0] >= bounds.X.Min && end_point[0] <= bounds.X.Max &&
end_point[1] >= bounds.Y.Min && end_point[1] <= bounds.Y.Max &&
end_point[2] >= bounds.Z.Min && end_point[2] <= bounds.Z.Max)
if (end_point.Pos[0] >= bounds.X.Min && end_point.Pos[0] <= bounds.X.Max &&
end_point.Pos[1] >= bounds.Y.Min && end_point.Pos[1] <= bounds.Y.Max &&
end_point.Pos[2] >= bounds.Z.Min && end_point.Pos[2] <= bounds.Z.Max)
{
res = 1;
}
@ -178,19 +177,24 @@ inline void Lagrangian::InitializeUniformSeeds(const vtkm::cont::DataSet& input)
auto portal1 = BasisParticles.GetPortalControl();
auto portal2 = BasisParticlesValidity.GetPortalControl();
vtkm::Id count = 0;
vtkm::Id count = 0, id = 0;
for (int x = 0; x < this->SeedRes[0]; x++)
{
vtkm::FloatDefault xi = static_cast<vtkm::FloatDefault>(x * x_spacing);
for (int y = 0; y < this->SeedRes[1]; y++)
{
vtkm::FloatDefault yi = static_cast<vtkm::FloatDefault>(y * y_spacing);
for (int z = 0; z < this->SeedRes[2]; z++)
{
vtkm::FloatDefault zi = static_cast<vtkm::FloatDefault>(z * z_spacing);
portal1.Set(count,
vtkm::Vec3f_64(bounds.X.Min + (x * x_spacing),
bounds.Y.Min + (y * y_spacing),
bounds.Z.Min + (z * z_spacing)));
vtkm::Particle(Vec3f(static_cast<vtkm::FloatDefault>(bounds.X.Min) + xi,
static_cast<vtkm::FloatDefault>(bounds.Y.Min) + yi,
static_cast<vtkm::FloatDefault>(bounds.Z.Min) + zi),
id));
portal2.Set(count, 1);
count++;
id++;
}
}
}
@ -223,7 +227,7 @@ inline VTKM_CONT vtkm::cont::DataSet Lagrangian::DoExecute(
throw vtkm::cont::ErrorFilterExecution(
"Write frequency can not be 0. Use SetWriteFrequency().");
}
vtkm::cont::ArrayHandle<vtkm::Vec<T, 3>> basisParticleArray;
vtkm::cont::ArrayHandle<vtkm::Particle> basisParticleArray;
vtkm::cont::ArrayCopy(BasisParticles, basisParticleArray);
cycle += 1;
@ -240,15 +244,13 @@ inline VTKM_CONT vtkm::cont::DataSet Lagrangian::DoExecute(
GridEvalType gridEval(coords, cells, field);
RK4Type rk4(gridEval, static_cast<vtkm::Float32>(this->stepSize));
res = particleadvection.Run(rk4, basisParticleArray, 1); // Taking a single step
auto particle_positions = res.positions;
auto particle_stepstaken = res.stepsTaken;
auto particles = res.Particles;
auto particlePortal = particles.GetPortalControl();
auto start_position = BasisParticlesOriginal.GetPortalControl();
auto end_position = particle_positions.GetPortalControl();
auto portal_stepstaken = particle_stepstaken.GetPortalControl();
auto portal_validity = BasisParticlesValidity.GetPortalControl();
vtkm::cont::DataSet outputData;
@ -262,25 +264,26 @@ inline VTKM_CONT vtkm::cont::DataSet Lagrangian::DoExecute(
std::vector<vtkm::UInt8> shapes;
std::vector<vtkm::IdComponent> numIndices;
for (vtkm::Id index = 0; index < res.positions.GetNumberOfValues(); index++)
for (vtkm::Id index = 0; index < particlePortal.GetNumberOfValues(); index++)
{
auto start_point = start_position.Get(index);
auto end_point = end_position.Get(index);
auto steps = portal_stepstaken.Get(index);
auto end_point = particlePortal.Get(index).Pos;
auto steps = particlePortal.Get(index).NumSteps;
if (steps > 0 && portal_validity.Get(index) == 1)
{
if (end_point[0] >= bounds.X.Min && end_point[0] <= bounds.X.Max &&
end_point[1] >= bounds.Y.Min && end_point[1] <= bounds.Y.Max &&
end_point[2] >= bounds.Z.Min && end_point[2] <= bounds.Z.Max)
if (bounds.Contains(end_point))
{
connectivity.push_back(connectivity_index);
connectivity.push_back(connectivity_index + 1);
connectivity_index += 2;
pointCoordinates.push_back(
vtkm::Vec<T, 3>((float)start_point[0], (float)start_point[1], (float)start_point[2]));
pointCoordinates.push_back(
vtkm::Vec<T, 3>((float)end_point[0], (float)end_point[1], (float)end_point[2]));
vtkm::Vec3f(static_cast<vtkm::FloatDefault>(start_point.Pos[0]),
static_cast<vtkm::FloatDefault>(start_point.Pos[1]),
static_cast<vtkm::FloatDefault>(start_point.Pos[2])));
pointCoordinates.push_back(vtkm::Vec3f(static_cast<vtkm::FloatDefault>(end_point[0]),
static_cast<vtkm::FloatDefault>(end_point[1]),
static_cast<vtkm::FloatDefault>(end_point[2])));
shapes.push_back(vtkm::CELL_SHAPE_LINE);
numIndices.push_back(2);
}
@ -308,14 +311,14 @@ inline VTKM_CONT vtkm::cont::DataSet Lagrangian::DoExecute(
}
else
{
vtkm::cont::ArrayCopy(particle_positions, BasisParticles);
vtkm::cont::ArrayCopy(particles, BasisParticles);
}
}
else
{
ValidityCheck check(bounds);
this->Invoke(check, particle_positions, particle_stepstaken, BasisParticlesValidity);
vtkm::cont::ArrayCopy(particle_positions, BasisParticles);
this->Invoke(check, particles, BasisParticlesValidity);
vtkm::cont::ArrayCopy(particles, BasisParticles);
}
return outputData;

22
vtkm/filter/LagrangianStructures.hxx
View File

@ -13,6 +13,7 @@
#include <vtkm/cont/ArrayCopy.h>
#include <vtkm/cont/ArrayHandleIndex.h>
#include <vtkm/cont/ErrorFilterExecution.h>
#include <vtkm/cont/Invoker.h>
#include <vtkm/worklet/particleadvection/GridEvaluators.h>
#include <vtkm/worklet/particleadvection/Integrators.h>
#include <vtkm/worklet/particleadvection/Particles.h>
@ -24,6 +25,23 @@ namespace vtkm
namespace filter
{
namespace detail
{
class ExtractParticlePosition : public vtkm::worklet::WorkletMapField
{
public:
using ControlSignature = void(FieldIn particle, FieldOut position);
using ExecutionSignature = void(_1, _2);
using InputDomain = _1;
VTKM_EXEC void operator()(const vtkm::Particle& particle, vtkm::Vec3f& pt) const
{
pt = particle.Pos;
}
};
} //detail
//-----------------------------------------------------------------------------
inline VTKM_CONT LagrangianStructures::LagrangianStructures()
: vtkm::filter::FilterDataSetWithField<LagrangianStructures>()
@ -104,7 +122,9 @@ inline VTKM_CONT vtkm::cont::DataSet LagrangianStructures::DoExecute(
vtkm::cont::ArrayHandle<vtkm::Vec3f> advectionPoints;
vtkm::cont::ArrayCopy(lcsInputPoints, advectionPoints);
advectionResult = particles.Run(integrator, advectionPoints, numberOfSteps);
lcsOutputPoints = advectionResult.positions;
vtkm::cont::Invoker invoke;
invoke(detail::ExtractParticlePosition{}, advectionResult.Particles, lcsOutputPoints);
}
// FTLE output field
vtkm::cont::ArrayHandle<vtkm::FloatDefault> outputField;

4
vtkm/filter/Pathline.h
View File

@ -47,7 +47,7 @@ public:
void SetNumberOfSteps(vtkm::Id n) { this->NumberOfSteps = n; }
VTKM_CONT
void SetSeeds(vtkm::cont::ArrayHandle<vtkm::Vec3f>& seeds);
void SetSeeds(vtkm::cont::ArrayHandle<vtkm::Particle>& seeds);
template <typename T, typename StorageType, typename DerivedPolicy>
VTKM_CONT vtkm::cont::DataSet DoExecute(
@ -71,7 +71,7 @@ private:
vtkm::FloatDefault NextTime;
vtkm::cont::DataSet NextDataSet;
vtkm::Id NumberOfSteps;
vtkm::cont::ArrayHandle<vtkm::Vec3f> Seeds;
vtkm::cont::ArrayHandle<vtkm::Particle> Seeds;
};
}
} // namespace vtkm::filter

8
vtkm/filter/Pathline.hxx
View File

@ -32,7 +32,7 @@ inline VTKM_CONT Pathline::Pathline()
}
//-----------------------------------------------------------------------------
inline VTKM_CONT void Pathline::SetSeeds(vtkm::cont::ArrayHandle<vtkm::Vec3f>& seeds)
inline VTKM_CONT void Pathline::SetSeeds(vtkm::cont::ArrayHandle<vtkm::Particle>& seeds)
{
this->Seeds = seeds;
}
@ -77,13 +77,13 @@ inline VTKM_CONT vtkm::cont::DataSet Pathline::DoExecute(
vtkm::worklet::Streamline streamline;
vtkm::worklet::StreamlineResult res;
vtkm::cont::ArrayHandle<vtkm::Vec<T, 3>> seedArray;
vtkm::cont::ArrayHandle<vtkm::Particle> seedArray;
vtkm::cont::ArrayCopy(this->Seeds, seedArray);
res = Worklet.Run(rk4, seedArray, this->NumberOfSteps);
vtkm::cont::DataSet outData;
vtkm::cont::CoordinateSystem outputCoords("coordinates", res.positions);
outData.SetCellSet(res.polyLines);
vtkm::cont::CoordinateSystem outputCoords("coordinates", res.Positions);
outData.SetCellSet(res.PolyLines);
outData.AddCoordinateSystem(outputCoords);
return outData;

4
vtkm/filter/StreamSurface.h
View File

@ -40,7 +40,7 @@ public:
void SetNumberOfSteps(vtkm::Id n) { this->NumberOfSteps = n; }
VTKM_CONT
void SetSeeds(vtkm::cont::ArrayHandle<vtkm::Vec3f>& seeds) { this->Seeds = seeds; }
void SetSeeds(vtkm::cont::ArrayHandle<vtkm::Particle>& seeds) { this->Seeds = seeds; }
template <typename T, typename StorageType, typename DerivedPolicy>
VTKM_CONT vtkm::cont::DataSet DoExecute(
@ -59,7 +59,7 @@ public:
private:
vtkm::Id NumberOfSteps;
vtkm::cont::ArrayHandle<vtkm::Vec3f> Seeds;
vtkm::cont::ArrayHandle<vtkm::Particle> Seeds;
vtkm::FloatDefault StepSize;
vtkm::worklet::StreamSurface Worklet;
};

6
vtkm/filter/StreamSurface.hxx
View File

@ -60,15 +60,15 @@ inline VTKM_CONT vtkm::cont::DataSet StreamSurface::DoExecute(
vtkm::worklet::Streamline streamline;
vtkm::cont::ArrayHandle<vtkm::Vec<T, 3>> seedArray;
vtkm::cont::ArrayHandle<vtkm::Particle> seedArray;
vtkm::cont::ArrayCopy(this->Seeds, seedArray);
auto res = streamline.Run(rk4, seedArray, this->NumberOfSteps);
//compute surface from streamlines
vtkm::cont::ArrayHandle<vtkm::Vec3f> srfPoints;
vtkm::cont::CellSetSingleType<> srfCells;
vtkm::cont::CoordinateSystem slCoords("coordinates", res.positions);
this->Worklet.Run(slCoords, res.polyLines, srfPoints, srfCells);
vtkm::cont::CoordinateSystem slCoords("coordinates", res.Positions);
this->Worklet.Run(slCoords, res.PolyLines, srfPoints, srfCells);
vtkm::cont::DataSet outData;
vtkm::cont::CoordinateSystem outputCoords("coordinates", srfPoints);

4
vtkm/filter/Streamline.h
View File

@ -39,7 +39,7 @@ public:
void SetNumberOfSteps(vtkm::Id n) { this->NumberOfSteps = n; }
VTKM_CONT
void SetSeeds(vtkm::cont::ArrayHandle<vtkm::Vec3f>& seeds);
void SetSeeds(vtkm::cont::ArrayHandle<vtkm::Particle>& seeds);
template <typename T, typename StorageType, typename DerivedPolicy>
VTKM_CONT vtkm::cont::DataSet DoExecute(
@ -59,7 +59,7 @@ public:
private:
vtkm::Id NumberOfSteps;
vtkm::FloatDefault StepSize;
vtkm::cont::ArrayHandle<vtkm::Vec3f> Seeds;
vtkm::cont::ArrayHandle<vtkm::Particle> Seeds;
vtkm::worklet::Streamline Worklet;
};
}

8
vtkm/filter/Streamline.hxx
View File

@ -30,7 +30,7 @@ inline VTKM_CONT Streamline::Streamline()
}
//-----------------------------------------------------------------------------
inline VTKM_CONT void Streamline::SetSeeds(vtkm::cont::ArrayHandle<vtkm::Vec3f>& seeds)
inline VTKM_CONT void Streamline::SetSeeds(vtkm::cont::ArrayHandle<vtkm::Particle>& seeds)
{
this->Seeds = seeds;
}
@ -67,13 +67,13 @@ inline VTKM_CONT vtkm::cont::DataSet Streamline::DoExecute(
vtkm::worklet::StreamlineResult res;
vtkm::cont::ArrayHandle<vtkm::Vec<T, 3>> seedArray;
vtkm::cont::ArrayHandle<vtkm::Particle> seedArray;
vtkm::cont::ArrayCopy(this->Seeds, seedArray);
res = this->Worklet.Run(rk4, seedArray, this->NumberOfSteps);
vtkm::cont::DataSet outData;
vtkm::cont::CoordinateSystem outputCoords("coordinates", res.positions);
outData.SetCellSet(res.polyLines);
vtkm::cont::CoordinateSystem outputCoords("coordinates", res.Positions);
outData.SetCellSet(res.PolyLines);
outData.AddCoordinateSystem(outputCoords);
return outData;

17
vtkm/filter/testing/UnitTestStreamSurfaceFilter.cxx
View File

@ -37,17 +37,16 @@ void TestStreamSurface()
{
std::cout << "Testing Stream Surface Filter" << std::endl;
using VecType = vtkm::Vec3f;
const vtkm::Id3 dims(5, 5, 5);
const vtkm::Vec3f vecX(1, 0, 0);
vtkm::cont::DataSet ds = CreateDataSet(dims, vecX);
vtkm::cont::ArrayHandle<VecType> seedArray;
std::vector<VecType> seeds(4);
seeds[0] = VecType(.1f, 1.0f, .2f);
seeds[1] = VecType(.1f, 2.0f, .1f);
seeds[2] = VecType(.1f, 3.0f, .3f);
seeds[3] = VecType(.1f, 3.5f, .2f);
vtkm::cont::ArrayHandle<vtkm::Particle> seedArray;
std::vector<vtkm::Particle> seeds(4);
seeds[0] = vtkm::Particle(vtkm::Vec3f(.1f, 1.0f, .2f), 0);
seeds[1] = vtkm::Particle(vtkm::Vec3f(.1f, 2.0f, .1f), 1);
seeds[2] = vtkm::Particle(vtkm::Vec3f(.1f, 3.0f, .3f), 2);
seeds[3] = vtkm::Particle(vtkm::Vec3f(.1f, 3.5f, .2f), 3);
seedArray = vtkm::cont::make_ArrayHandle(seeds);
@ -65,10 +64,10 @@ void TestStreamSurface()
"Wrong number of coordinate systems in the output dataset");
vtkm::cont::CoordinateSystem coords = output.GetCoordinateSystem();
VTKM_TEST_ASSERT(coords.GetNumberOfPoints() == 80, "Wrong number of coordinates");
VTKM_TEST_ASSERT(coords.GetNumberOfPoints() == 84, "Wrong number of coordinates");
vtkm::cont::DynamicCellSet dcells = output.GetCellSet();
VTKM_TEST_ASSERT(dcells.GetNumberOfCells() == 114, "Wrong number of cells");
VTKM_TEST_ASSERT(dcells.GetNumberOfCells() == 120, "Wrong number of cells");
}
}

24
vtkm/filter/testing/UnitTestStreamlineFilter.cxx
View File

@ -39,11 +39,11 @@ void TestStreamline()
const vtkm::Vec3f vecX(1, 0, 0);
vtkm::cont::DataSet ds = CreateDataSet(dims, vecX);
vtkm::cont::ArrayHandle<vtkm::Vec3f> seedArray;
std::vector<vtkm::Vec3f> seeds(3);
seeds[0] = vtkm::Vec3f(.2f, 1.0f, .2f);
seeds[1] = vtkm::Vec3f(.2f, 2.0f, .2f);
seeds[2] = vtkm::Vec3f(.2f, 3.0f, .2f);
vtkm::cont::ArrayHandle<vtkm::Particle> seedArray;
std::vector<vtkm::Particle> seeds(3);
seeds[0] = vtkm::Particle(vtkm::Vec3f(.2f, 1.0f, .2f), 0);
seeds[1] = vtkm::Particle(vtkm::Vec3f(.2f, 2.0f, .2f), 1);
seeds[2] = vtkm::Particle(vtkm::Vec3f(.2f, 3.0f, .2f), 2);
seedArray = vtkm::cont::make_ArrayHandle(seeds);
@ -61,7 +61,7 @@ void TestStreamline()
"Wrong number of coordinate systems in the output dataset");
vtkm::cont::CoordinateSystem coords = output.GetCoordinateSystem();
VTKM_TEST_ASSERT(coords.GetNumberOfPoints() == 60, "Wrong number of coordinates");
VTKM_TEST_ASSERT(coords.GetNumberOfPoints() == 63, "Wrong number of coordinates");
vtkm::cont::DynamicCellSet dcells = output.GetCellSet();
VTKM_TEST_ASSERT(dcells.GetNumberOfCells() == 3, "Wrong number of cells");
@ -76,11 +76,11 @@ void TestPathline()
vtkm::cont::DataSet ds1 = CreateDataSet(dims, vecX);
vtkm::cont::DataSet ds2 = CreateDataSet(dims, vecY);
vtkm::cont::ArrayHandle<vtkm::Vec3f> seedArray;
std::vector<vtkm::Vec3f> seeds(3);
seeds[0] = vtkm::Vec3f(.2f, 1.0f, .2f);
seeds[1] = vtkm::Vec3f(.2f, 2.0f, .2f);
seeds[2] = vtkm::Vec3f(.2f, 3.0f, .2f);
vtkm::cont::ArrayHandle<vtkm::Particle> seedArray;
std::vector<vtkm::Particle> seeds(3);
seeds[0] = vtkm::Particle(vtkm::Vec3f(.2f, 1.0f, .2f), 0);
seeds[1] = vtkm::Particle(vtkm::Vec3f(.2f, 2.0f, .2f), 1);
seeds[2] = vtkm::Particle(vtkm::Vec3f(.2f, 3.0f, .2f), 2);
seedArray = vtkm::cont::make_ArrayHandle(seeds);
@ -98,7 +98,7 @@ void TestPathline()
//Validate the result is correct.
vtkm::cont::CoordinateSystem coords = output.GetCoordinateSystem();
VTKM_TEST_ASSERT(coords.GetNumberOfPoints() == 60, "Wrong number of coordinates");
VTKM_TEST_ASSERT(coords.GetNumberOfPoints() == 63, "Wrong number of coordinates");
vtkm::cont::DynamicCellSet dcells = output.GetCellSet();
VTKM_TEST_ASSERT(dcells.GetNumberOfCells() == 3, "Wrong number of cells");

265
vtkm/worklet/ParticleAdvection.h
View File

@ -14,6 +14,7 @@
#include <vtkm/Types.h>
#include <vtkm/cont/ArrayCopy.h>
#include <vtkm/cont/ArrayHandle.h>
#include <vtkm/cont/Invoker.h>
#include <vtkm/worklet/particleadvection/ParticleAdvectionWorklets.h>
namespace vtkm
@ -21,40 +22,45 @@ namespace vtkm
namespace worklet
{
namespace detail
{
class CopyToParticle : public vtkm::worklet::WorkletMapField
{
public:
using ControlSignature =
void(FieldIn pt, FieldIn id, FieldIn time, FieldIn step, FieldOut particle);
using ExecutionSignature = void(_1, _2, _3, _4, _5);
using InputDomain = _1;
VTKM_EXEC void operator()(const vtkm::Vec3f& pt,
const vtkm::Id& id,
const vtkm::FloatDefault& time,
const vtkm::Id& step,
vtkm::Particle& particle) const
{
particle.Pos = pt;
particle.ID = id;
particle.Time = time;
particle.NumSteps = step;
particle.Status.SetOk();
}
};
} //detail
struct ParticleAdvectionResult
{
ParticleAdvectionResult()
: positions()
, status()
, stepsTaken()
, times()
: Particles()
{
}
ParticleAdvectionResult(const vtkm::cont::ArrayHandle<vtkm::Vec3f>& pos,
const vtkm::cont::ArrayHandle<vtkm::Id>& stat,
const vtkm::cont::ArrayHandle<vtkm::Id>& steps)
: positions(pos)
, status(stat)
, stepsTaken(steps)
ParticleAdvectionResult(const vtkm::cont::ArrayHandle<vtkm::Particle>& p)
: Particles(p)
{
}
ParticleAdvectionResult(const vtkm::cont::ArrayHandle<vtkm::Vec3f>& pos,
const vtkm::cont::ArrayHandle<vtkm::Id>& stat,
const vtkm::cont::ArrayHandle<vtkm::Id>& steps,
const vtkm::cont::ArrayHandle<vtkm::FloatDefault>& timeArray)
: positions(pos)
, status(stat)
, stepsTaken(steps)
, times(timeArray)
{
}
vtkm::cont::ArrayHandle<vtkm::Vec3f> positions;
vtkm::cont::ArrayHandle<vtkm::Id> status;
vtkm::cont::ArrayHandle<vtkm::Id> stepsTaken;
vtkm::cont::ArrayHandle<vtkm::FloatDefault> times;
vtkm::cont::ArrayHandle<vtkm::Particle> Particles;
};
class ParticleAdvection
@ -62,120 +68,66 @@ class ParticleAdvection
public:
ParticleAdvection() {}
template <typename IntegratorType, typename FieldType, typename PointStorage>
template <typename IntegratorType, typename ParticleStorage>
ParticleAdvectionResult Run(const IntegratorType& it,
vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>, PointStorage>& pts,
const vtkm::Id& nSteps)
{
vtkm::Id numSeeds = static_cast<vtkm::Id>(pts.GetNumberOfValues());
vtkm::cont::ArrayHandle<vtkm::Id> stepsTaken;
vtkm::cont::ArrayHandle<vtkm::FloatDefault> timeArray;
//Allocate status and steps arrays.
vtkm::cont::ArrayHandleConstant<vtkm::Id> init(0, numSeeds);
vtkm::cont::ArrayCopy(init, stepsTaken);
//Allocate memory to store the time for temporal integration.
vtkm::cont::ArrayHandleConstant<vtkm::FloatDefault> time(0, numSeeds);
vtkm::cont::ArrayCopy(time, timeArray);
return Run(it, pts, stepsTaken, timeArray, nSteps);
}
template <typename IntegratorType, typename FieldType, typename PointStorage>
ParticleAdvectionResult Run(const IntegratorType& it,
vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>, PointStorage>& pts,
vtkm::cont::ArrayHandle<vtkm::Id>& inputSteps,
const vtkm::Id& nSteps)
{
vtkm::Id numSeeds = static_cast<vtkm::Id>(pts.GetNumberOfValues());
vtkm::cont::ArrayHandle<vtkm::FloatDefault> timeArray;
//Allocate memory to store the time for temporal integration.
vtkm::cont::ArrayHandleConstant<vtkm::FloatDefault> time(0, numSeeds);
timeArray.Allocate(numSeeds);
vtkm::cont::ArrayCopy(time, timeArray);
return Run(it, pts, inputSteps, timeArray, nSteps);
}
template <typename IntegratorType>
ParticleAdvectionResult Run(const IntegratorType& it,
vtkm::cont::ArrayHandle<vtkm::Vec3f>& pts,
vtkm::cont::ArrayHandle<vtkm::Id>& inputSteps,
vtkm::cont::ArrayHandle<vtkm::FloatDefault>& inputTime,
const vtkm::Id& nSteps)
vtkm::cont::ArrayHandle<vtkm::Particle, ParticleStorage>& particles,
vtkm::Id MaxSteps)
{
vtkm::worklet::particleadvection::ParticleAdvectionWorklet<IntegratorType> worklet;
vtkm::Id numSeeds = static_cast<vtkm::Id>(pts.GetNumberOfValues());
vtkm::cont::ArrayHandle<vtkm::Id> status;
//Allocate status arrays.
vtkm::cont::ArrayHandleConstant<vtkm::Id> statusOK(static_cast<vtkm::Id>(1), numSeeds);
status.Allocate(numSeeds);
vtkm::cont::ArrayCopy(statusOK, status);
worklet.Run(it, pts, nSteps, status, inputSteps, inputTime);
//Create output.
return ParticleAdvectionResult(pts, status, inputSteps, inputTime);
worklet.Run(it, particles, MaxSteps);
return ParticleAdvectionResult(particles);
}
template <typename IntegratorType, typename FieldType, typename PointStorage>
template <typename IntegratorType, typename PointStorage>
ParticleAdvectionResult Run(const IntegratorType& it,
vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>, PointStorage>& pts,
vtkm::cont::ArrayHandle<vtkm::Id>& inputSteps,
vtkm::cont::ArrayHandle<vtkm::FloatDefault>& inputTime,
const vtkm::Id& nSteps)
const vtkm::cont::ArrayHandle<vtkm::Vec3f, PointStorage>& points,
vtkm::Id MaxSteps)
{
vtkm::cont::ArrayHandle<vtkm::Vec3f> ptsCopy;
vtkm::cont::ArrayCopy(pts, ptsCopy);
return Run(it, ptsCopy, inputSteps, inputTime, nSteps);
vtkm::worklet::particleadvection::ParticleAdvectionWorklet<IntegratorType> worklet;
vtkm::cont::ArrayHandle<vtkm::Particle> particles;
vtkm::cont::ArrayHandle<vtkm::Id> step, ids;
vtkm::cont::ArrayHandle<vtkm::FloatDefault> time;
vtkm::cont::Invoker invoke;
vtkm::Id numPts = points.GetNumberOfValues();
vtkm::cont::ArrayHandleConstant<vtkm::Id> s(0, numPts);
vtkm::cont::ArrayHandleConstant<vtkm::FloatDefault> t(0, numPts);
vtkm::cont::ArrayHandleCounting<vtkm::Id> id(0, 1, numPts);
//Copy input to vtkm::Particle
vtkm::cont::ArrayCopy(s, step);
vtkm::cont::ArrayCopy(t, time);
vtkm::cont::ArrayCopy(id, ids);
invoke(detail::CopyToParticle{}, points, ids, time, step, particles);
worklet.Run(it, particles, MaxSteps);
return ParticleAdvectionResult(particles);
}
};
struct StreamlineResult
{
StreamlineResult()
: positions()
, polyLines()
, status()
, stepsTaken()
, times()
: Particles()
, Positions()
, PolyLines()
{
}
StreamlineResult(const vtkm::cont::ArrayHandle<vtkm::Vec3f>& pos,
const vtkm::cont::CellSetExplicit<>& lines,
const vtkm::cont::ArrayHandle<vtkm::Id>& stat,
const vtkm::cont::ArrayHandle<vtkm::Id>& steps)
: positions(pos)
, polyLines(lines)
, status(stat)
, stepsTaken(steps)
StreamlineResult(const vtkm::cont::ArrayHandle<vtkm::Particle>& part,
const vtkm::cont::ArrayHandle<vtkm::Vec3f>& pos,
const vtkm::cont::CellSetExplicit<>& lines)
: Particles(part)
, Positions(pos)
, PolyLines(lines)
{
}
StreamlineResult(const vtkm::cont::ArrayHandle<vtkm::Vec3f>& pos,
const vtkm::cont::CellSetExplicit<>& lines,
const vtkm::cont::ArrayHandle<vtkm::Id>& stat,
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;
vtkm::cont::ArrayHandle<vtkm::Particle> Particles;
vtkm::cont::ArrayHandle<vtkm::Vec3f> Positions;
vtkm::cont::CellSetExplicit<> PolyLines;
};
class Streamline
@ -183,84 +135,19 @@ class Streamline
public:
Streamline() {}
template <typename IntegratorType, typename FieldType, typename PointStorage>
template <typename IntegratorType, typename ParticleStorage>
StreamlineResult Run(const IntegratorType& it,
vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>, PointStorage>& seedArray,
const vtkm::Id& nSteps)
{
vtkm::Id numSeeds = seedArray.GetNumberOfValues();
//Allocate status and steps arrays.
vtkm::cont::ArrayHandle<vtkm::Id> status, steps;
vtkm::cont::ArrayHandleConstant<vtkm::Id> statusOK(static_cast<vtkm::Id>(1), numSeeds);
vtkm::cont::ArrayCopy(statusOK, status);
vtkm::cont::ArrayHandleConstant<vtkm::Id> zero(0, numSeeds);
vtkm::cont::ArrayCopy(zero, steps);
//Allocate memory to store the time for temporal integration.
vtkm::cont::ArrayHandle<vtkm::FloatDefault> timeArray;
vtkm::cont::ArrayHandleConstant<vtkm::FloatDefault> time(0, numSeeds);
vtkm::cont::ArrayCopy(time, timeArray);
return Run(it, seedArray, steps, timeArray, nSteps);
}
template <typename IntegratorType, typename FieldType, typename PointStorage>
StreamlineResult Run(const IntegratorType& it,
vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>, PointStorage>& seedArray,
vtkm::cont::ArrayHandle<vtkm::Id>& inputSteps,
const vtkm::Id& nSteps)
{
vtkm::Id numSeeds = seedArray.GetNumberOfValues();
//Allocate and initializr status array.
vtkm::cont::ArrayHandle<vtkm::Id> status;
vtkm::cont::ArrayHandleConstant<vtkm::Id> statusOK(static_cast<vtkm::Id>(1), numSeeds);
vtkm::cont::ArrayCopy(statusOK, status);
//Allocate memory to store the time for temporal integration.
vtkm::cont::ArrayHandle<vtkm::FloatDefault> timeArray;
vtkm::cont::ArrayHandleConstant<vtkm::FloatDefault> time(0, numSeeds);
vtkm::cont::ArrayCopy(time, timeArray);
return Run(it, seedArray, inputSteps, timeArray, nSteps);
}
template <typename IntegratorType>
StreamlineResult Run(const IntegratorType& it,
vtkm::cont::ArrayHandle<vtkm::Vec3f>& seedArray,
vtkm::cont::ArrayHandle<vtkm::Id>& inputSteps,
vtkm::cont::ArrayHandle<vtkm::FloatDefault>& inputTime,
const vtkm::Id& nSteps)
vtkm::cont::ArrayHandle<vtkm::Particle, ParticleStorage>& particles,
vtkm::Id MaxSteps)
{
vtkm::worklet::particleadvection::StreamlineWorklet<IntegratorType> worklet;
vtkm::cont::ArrayHandle<vtkm::Vec3f> positions;
vtkm::cont::CellSetExplicit<> polyLines;
//Allocate and initialize status array.
vtkm::Id numSeeds = seedArray.GetNumberOfValues();
vtkm::cont::ArrayHandle<vtkm::Id> status;
vtkm::cont::ArrayHandleConstant<vtkm::Id> statusOK(static_cast<vtkm::Id>(1), numSeeds);
vtkm::cont::ArrayCopy(statusOK, status);
worklet.Run(it, particles, MaxSteps, positions, polyLines);
worklet.Run(it, seedArray, nSteps, positions, polyLines, status, inputSteps, inputTime);
return StreamlineResult(positions, polyLines, status, inputSteps, inputTime);
}
template <typename IntegratorType, typename FieldType, typename PointStorage>
StreamlineResult Run(const IntegratorType& it,
vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>, PointStorage>& seedArray,
vtkm::cont::ArrayHandle<vtkm::Id>& inputSteps,
vtkm::cont::ArrayHandle<vtkm::FloatDefault>& inputTime,
const vtkm::Id& nSteps)
{
vtkm::cont::ArrayHandle<vtkm::Vec3f> seedCopy;
vtkm::cont::ArrayCopy(seedArray, seedCopy);
return Run(it, seedCopy, inputSteps, inputTime, nSteps);
return StreamlineResult(particles, positions, polyLines);
}
};
}

3
vtkm/worklet/particleadvection/CMakeLists.txt
View File

@ -10,9 +10,10 @@
set(headers
CellInterpolationHelper.h
EvaluatorStatus.h
GridEvaluators.h
GridEvaluatorStatus.h
Integrators.h
IntegratorStatus.h
Particles.h
ParticleAdvectionWorklets.h
TemporalGridEvaluators.h

34
vtkm/worklet/particleadvection/EvaluatorStatus.h
View File

@ -1,34 +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_EvaluatorStatus_h
#define vtk_m_worklet_particleadvection_EvaluatorStatus_h
namespace vtkm
{
namespace worklet
{
namespace particleadvection
{
enum class EvaluatorStatus
{
SUCCESS = 0,
OUTSIDE_SPATIAL_BOUNDS,
OUTSIDE_TEMPORAL_BOUNDS
};
} //namespace particleadvection
} //namespace worklet
} //namespace vtkm
#endif // vtk_m_worklet_particleadvection_EvaluatorStatus_h

67
vtkm/worklet/particleadvection/GridEvaluatorStatus.h Normal file
View File

@ -0,0 +1,67 @@
//============================================================================
// 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_GridEvaluatorStatus_h
#define vtk_m_worklet_particleadvection_GridEvaluatorStatus_h
#include <vtkm/Bitset.h>
#include <vtkm/Types.h>
#include <vtkm/VectorAnalysis.h>
#include <vtkm/cont/ArrayHandle.h>
#include <vtkm/cont/CellLocator.h>
#include <vtkm/cont/CellLocatorRectilinearGrid.h>
#include <vtkm/cont/CellLocatorUniformBins.h>
#include <vtkm/cont/CellLocatorUniformGrid.h>
#include <vtkm/cont/CellSetStructured.h>
#include <vtkm/cont/DataSet.h>
#include <vtkm/cont/DeviceAdapter.h>
#include <vtkm/worklet/particleadvection/CellInterpolationHelper.h>
#include <vtkm/worklet/particleadvection/Integrators.h>
namespace vtkm
{
namespace worklet
{
namespace particleadvection
{
class GridEvaluatorStatus : public vtkm::Bitset<vtkm::UInt8>
{
public:
VTKM_EXEC_CONT GridEvaluatorStatus(){};
VTKM_EXEC_CONT GridEvaluatorStatus(bool ok, bool spatial, bool temporal)
{
this->set(this->SUCCESS_BIT, ok);
this->set(this->SPATIAL_BOUNDS_BIT, spatial);
this->set(this->TEMPORAL_BOUNDS_BIT, temporal);
};
VTKM_EXEC_CONT void SetOk() { this->set(this->SUCCESS_BIT); }
VTKM_EXEC_CONT bool CheckOk() const { return this->test(this->SUCCESS_BIT); }
VTKM_EXEC_CONT void SetFail() { this->reset(this->SUCCESS_BIT); }
VTKM_EXEC_CONT bool CheckFail() const { return !this->test(this->SUCCESS_BIT); }
VTKM_EXEC_CONT void SetSpatialBounds() { this->set(this->SPATIAL_BOUNDS_BIT); }
VTKM_EXEC_CONT bool CheckSpatialBounds() const { return this->test(this->SPATIAL_BOUNDS_BIT); }
VTKM_EXEC_CONT void SetTemporalBounds() { this->set(this->TEMPORAL_BOUNDS_BIT); }
VTKM_EXEC_CONT bool CheckTemporalBounds() const { return this->test(this->TEMPORAL_BOUNDS_BIT); }
private:
static constexpr vtkm::Id SUCCESS_BIT = 0;
static constexpr vtkm::Id SPATIAL_BOUNDS_BIT = 1;
static constexpr vtkm::Id TEMPORAL_BOUNDS_BIT = 2;
};
}
}
}
#endif // vtk_m_worklet_particleadvection_GridEvaluatorStatus_h

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

@ -11,6 +11,7 @@
#ifndef vtk_m_worklet_particleadvection_GridEvaluators_h
#define vtk_m_worklet_particleadvection_GridEvaluators_h
#include <vtkm/Bitset.h>
#include <vtkm/Types.h>
#include <vtkm/VectorAnalysis.h>
#include <vtkm/cont/ArrayHandle.h>
@ -23,7 +24,7 @@
#include <vtkm/cont/DeviceAdapter.h>
#include <vtkm/worklet/particleadvection/CellInterpolationHelper.h>
#include <vtkm/worklet/particleadvection/EvaluatorStatus.h>
#include <vtkm/worklet/particleadvection/GridEvaluatorStatus.h>
#include <vtkm/worklet/particleadvection/Integrators.h>
namespace vtkm
@ -81,22 +82,28 @@ public:
}
template <typename Point>
VTKM_EXEC EvaluatorStatus Evaluate(const Point& pos,
vtkm::FloatDefault vtkmNotUsed(time),
Point& out) const
VTKM_EXEC GridEvaluatorStatus Evaluate(const Point& pos,
vtkm::FloatDefault vtkmNotUsed(time),
Point& out) const
{
return this->Evaluate(pos, out);
}
template <typename Point>
VTKM_EXEC EvaluatorStatus Evaluate(const Point& point, Point& out) const
VTKM_EXEC GridEvaluatorStatus Evaluate(const Point& point, Point& out) const
{
vtkm::Id cellId;
Point parametric;
vtkm::exec::FunctorBase tmp;
GridEvaluatorStatus status;
Locator->FindCell(point, cellId, parametric, tmp);
if (cellId == -1)
return EvaluatorStatus::OUTSIDE_SPATIAL_BOUNDS;
{
status.SetFail();
status.SetSpatialBounds();
return status;
}
vtkm::UInt8 cellShape;
vtkm::IdComponent nVerts;
@ -108,7 +115,8 @@ public:
fieldValues.Append(Field.Get(ptIndices[i]));
out = vtkm::exec::CellInterpolate(fieldValues, parametric, cellShape, tmp);
return EvaluatorStatus::SUCCESS;
status.SetOk();
return status;
}
private:

79
vtkm/worklet/particleadvection/IntegratorStatus.h Normal file
View File

@ -0,0 +1,79 @@
//=============================================================================
//
// 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_Integrator_Status_h
#define vtk_m_worklet_particleadvection_Integrator_Status_h
#include <iomanip>
#include <limits>
#include <vtkm/Bitset.h>
#include <vtkm/TypeTraits.h>
#include <vtkm/Types.h>
#include <vtkm/VectorAnalysis.h>
namespace vtkm
{
namespace worklet
{
namespace particleadvection
{
class IntegratorStatus : public vtkm::Bitset<vtkm::UInt8>
{
public:
VTKM_EXEC_CONT IntegratorStatus() {}
VTKM_EXEC_CONT IntegratorStatus(const bool& ok, const bool& spatial, const bool& temporal)
{
this->set(this->SUCCESS_BIT, ok);
this->set(this->SPATIAL_BOUNDS_BIT, spatial);
this->set(this->TEMPORAL_BOUNDS_BIT, temporal);
}
VTKM_EXEC_CONT IntegratorStatus(const GridEvaluatorStatus& es)
{
this->set(this->SUCCESS_BIT, es.CheckOk());
this->set(this->SPATIAL_BOUNDS_BIT, es.CheckSpatialBounds());
this->set(this->TEMPORAL_BOUNDS_BIT, es.CheckTemporalBounds());
}
VTKM_EXEC_CONT void SetOk() { this->set(this->SUCCESS_BIT); }
VTKM_EXEC_CONT bool CheckOk() const { return this->test(this->SUCCESS_BIT); }
VTKM_EXEC_CONT void SetFail() { this->reset(this->SUCCESS_BIT); }
VTKM_EXEC_CONT bool CheckFail() const { return !this->test(this->SUCCESS_BIT); }
VTKM_EXEC_CONT void SetSpatialBounds() { this->set(this->SPATIAL_BOUNDS_BIT); }
VTKM_EXEC_CONT bool CheckSpatialBounds() const { return this->test(this->SPATIAL_BOUNDS_BIT); }
VTKM_EXEC_CONT void SetTemporalBounds() { this->set(this->TEMPORAL_BOUNDS_BIT); }
VTKM_EXEC_CONT bool CheckTemporalBounds() const { return this->test(this->TEMPORAL_BOUNDS_BIT); }
private:
static constexpr vtkm::Id SUCCESS_BIT = 0;
static constexpr vtkm::Id SPATIAL_BOUNDS_BIT = 1;
static constexpr vtkm::Id TEMPORAL_BOUNDS_BIT = 2;
};
inline VTKM_CONT std::ostream& operator<<(std::ostream& s, const IntegratorStatus& status)
{
s << "[" << status.CheckOk() << " " << status.CheckSpatialBounds() << " "
<< status.CheckTemporalBounds() << "]";
return s;
}
}
}
}
#endif // vtk_m_worklet_particleadvection_IntegratorStatus_h

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

@ -13,8 +13,10 @@
#ifndef vtk_m_worklet_particleadvection_Integrators_h
#define vtk_m_worklet_particleadvection_Integrators_h
#include <iomanip>
#include <limits>
#include <vtkm/Bitset.h>
#include <vtkm/TypeTraits.h>
#include <vtkm/Types.h>
#include <vtkm/VectorAnalysis.h>
@ -22,7 +24,8 @@
#include <vtkm/cont/DataSet.h>
#include <vtkm/cont/VirtualObjectHandle.h>
#include <vtkm/worklet/particleadvection/EvaluatorStatus.h>
#include <vtkm/worklet/particleadvection/GridEvaluators.h>
#include <vtkm/worklet/particleadvection/IntegratorStatus.h>
#include <vtkm/worklet/particleadvection/Particles.h>
namespace vtkm
@ -31,13 +34,6 @@ namespace worklet
{
namespace particleadvection
{
enum class IntegratorStatus
{
SUCCESS = 0,
OUTSIDE_SPATIAL_BOUNDS,
OUTSIDE_TEMPORAL_BOUNDS,
FAIL
};
class Integrator : public vtkm::cont::ExecutionObjectBase
{
@ -73,22 +69,6 @@ public:
vtkm::FloatDefault& time,
vtkm::Vec3f& outpos) const = 0;
VTKM_EXEC
IntegratorStatus ConvertToIntegratorStatus(EvaluatorStatus status) const
{
switch (status)
{
case EvaluatorStatus::SUCCESS:
return IntegratorStatus::SUCCESS;
case EvaluatorStatus::OUTSIDE_SPATIAL_BOUNDS:
return IntegratorStatus::OUTSIDE_SPATIAL_BOUNDS;
case EvaluatorStatus::OUTSIDE_TEMPORAL_BOUNDS:
return IntegratorStatus::OUTSIDE_TEMPORAL_BOUNDS;
default:
return IntegratorStatus::FAIL;
}
}
protected:
vtkm::FloatDefault StepLength = 1.0f;
vtkm::FloatDefault Tolerance = 0.001f;
@ -135,20 +115,30 @@ protected:
{
// If particle is out of either spatial or temporal boundary to begin with,
// then return the corresponding status.
IntegratorStatus status;
if (!this->Evaluator.IsWithinSpatialBoundary(inpos))
return IntegratorStatus::OUTSIDE_SPATIAL_BOUNDS;
{
status.SetFail();
status.SetSpatialBounds();
return status;
}
if (!this->Evaluator.IsWithinTemporalBoundary(time))
return IntegratorStatus::OUTSIDE_TEMPORAL_BOUNDS;
{
status.SetFail();
status.SetTemporalBounds();
return status;
}
vtkm::Vec3f velocity;
IntegratorStatus status = CheckStep(inpos, this->StepLength, time, velocity);
if (status == IntegratorStatus::SUCCESS)
status = CheckStep(inpos, this->StepLength, time, velocity);
if (status.CheckOk())
{
outpos = inpos + StepLength * velocity;
time += StepLength;
}
else
outpos = inpos;
return status;
}
@ -160,12 +150,12 @@ protected:
if (!this->Evaluator.IsWithinSpatialBoundary(inpos))
{
outpos = inpos;
return IntegratorStatus::OUTSIDE_SPATIAL_BOUNDS;
return IntegratorStatus(false, true, false);
}
if (!this->Evaluator.IsWithinTemporalBoundary(time))
{
outpos = inpos;
return IntegratorStatus::OUTSIDE_TEMPORAL_BOUNDS;
return IntegratorStatus(false, false, true);
}
//Stepping by this->StepLength goes beyond the bounds of the dataset.
@ -180,9 +170,11 @@ protected:
vtkm::Vec3f currPos(inpos), velocity;
vtkm::FloatDefault currTime = time;
this->Evaluator.Evaluate(currPos, time, velocity);
auto evalStatus = this->Evaluator.Evaluate(currPos, time, velocity);
if (evalStatus.CheckFail())
return IntegratorStatus(evalStatus);
const vtkm::FloatDefault eps = 10 * vtkm::Epsilon<vtkm::FloatDefault>();
const vtkm::FloatDefault eps = vtkm::Epsilon<vtkm::FloatDefault>();
vtkm::FloatDefault div = 1;
for (int i = 0; i < 50; i++)
{
@ -190,7 +182,7 @@ protected:
vtkm::FloatDefault stepCurr = stepShort + (this->StepLength / div);
//See if we can step by stepCurr.
IntegratorStatus status = this->CheckStep(inpos, stepCurr, time, velocity);
if (status == IntegratorStatus::SUCCESS)
if (status.CheckOk())
{
currPos = inpos + velocity * stepShort;
stepShort = stepCurr;
@ -204,11 +196,13 @@ protected:
}
//Take Euler step.
time = currTime + stepShort;
this->Evaluator.Evaluate(currPos, time, velocity);
currTime = time + stepShort;
evalStatus = this->Evaluator.Evaluate(currPos, currTime, velocity);
if (evalStatus.CheckFail())
return IntegratorStatus(evalStatus);
outpos = currPos + stepLong * velocity;
return IntegratorStatus::OUTSIDE_SPATIAL_BOUNDS;
return IntegratorStatus(true, true, !this->Evaluator.IsWithinTemporalBoundary(time));
}
VTKM_EXEC
@ -301,22 +295,22 @@ public:
vtkm::Vec3f k1 = vtkm::TypeTraits<vtkm::Vec3f>::ZeroInitialization();
vtkm::Vec3f k2 = k1, k3 = k1, k4 = k1;
EvaluatorStatus status;
status = this->Evaluator.Evaluate(inpos, time, k1);
if (status != EvaluatorStatus::SUCCESS)
return this->ConvertToIntegratorStatus(status);
status = this->Evaluator.Evaluate(inpos + var1 * k1, var2, k2);
if (status != EvaluatorStatus::SUCCESS)
return this->ConvertToIntegratorStatus(status);
status = this->Evaluator.Evaluate(inpos + var1 * k2, var2, k3);
if (status != EvaluatorStatus::SUCCESS)
return this->ConvertToIntegratorStatus(status);
status = this->Evaluator.Evaluate(inpos + stepLength * k3, var3, k4);
if (status != EvaluatorStatus::SUCCESS)
return this->ConvertToIntegratorStatus(status);
GridEvaluatorStatus evalStatus;
evalStatus = this->Evaluator.Evaluate(inpos, time, k1);
if (evalStatus.CheckFail())
return IntegratorStatus(evalStatus);
evalStatus = this->Evaluator.Evaluate(inpos + var1 * k1, var2, k2);
if (evalStatus.CheckFail())
return IntegratorStatus(evalStatus);
evalStatus = this->Evaluator.Evaluate(inpos + var1 * k2, var2, k3);
if (evalStatus.CheckFail())
return IntegratorStatus(evalStatus);
evalStatus = this->Evaluator.Evaluate(inpos + stepLength * k3, var3, k4);
if (evalStatus.CheckFail())
return IntegratorStatus(evalStatus);
velocity = (k1 + 2 * k2 + 2 * k3 + k4) / 6.0f;
return IntegratorStatus::SUCCESS;
return IntegratorStatus(true, false, evalStatus.CheckTemporalBounds());
}
};
@ -378,8 +372,8 @@ public:
vtkm::FloatDefault time,
vtkm::Vec3f& velocity) const
{
EvaluatorStatus status = this->Evaluator.Evaluate(inpos, time, velocity);
return this->ConvertToIntegratorStatus(status);
GridEvaluatorStatus status = this->Evaluator.Evaluate(inpos, time, velocity);
return IntegratorStatus(status);
}
};

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

@ -20,9 +20,11 @@
#include <vtkm/cont/CellSetExplicit.h>
#include <vtkm/cont/ExecutionObjectBase.h>
#include <vtkm/Particle.h>
#include <vtkm/worklet/DispatcherMapField.h>
#include <vtkm/worklet/particleadvection/Integrators.h>
#include <vtkm/worklet/particleadvection/Particles.h>
//#include <vtkm/worklet/particleadvection/ParticlesAOS.h>
namespace vtkm
{
@ -34,67 +36,68 @@ namespace particleadvection
class ParticleAdvectWorklet : public vtkm::worklet::WorkletMapField
{
public:
using ControlSignature = void(FieldIn idx, ExecObject integrator, ExecObject integralCurve);
using ExecutionSignature = void(_1, _2, _3);
using ControlSignature = void(FieldIn idx,
ExecObject integrator,
ExecObject integralCurve,
FieldIn maxSteps);
using ExecutionSignature = void(_1 idx, _2 integrator, _3 integralCurve, _4 maxSteps);
using InputDomain = _1;
template <typename IntegratorType, typename IntegralCurveType>
VTKM_EXEC void operator()(const vtkm::Id& idx,
const IntegratorType* integrator,
IntegralCurveType& integralCurve) const
IntegralCurveType& integralCurve,
const vtkm::Id& maxSteps) const
{
vtkm::Vec3f inpos = integralCurve.GetPos(idx);
vtkm::Vec3f outpos;
vtkm::FloatDefault time = integralCurve.GetTime(idx);
IntegratorStatus status;
vtkm::Particle particle = integralCurve.GetParticle(idx);
vtkm::Vec3f inpos = particle.Pos;
vtkm::FloatDefault time = particle.Time;
bool tookAnySteps = false;
while (!integralCurve.Done(idx))
//the integrator status needs to be more robust:
// 1. you could have success AND at temporal boundary.
// 2. could you have success AND at spatial?
// 3. all three?
integralCurve.PreStepUpdate(idx);
do
{
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 == IntegratorStatus::SUCCESS)
//vtkm::Particle p = integralCurve.GetParticle(idx);
//std::cout<<idx<<": "<<inpos<<" #"<<p.NumSteps<<" "<<p.Status<<std::endl;
vtkm::Vec3f outpos;
auto status = integrator->Step(inpos, time, outpos);
if (status.CheckOk())
{
integralCurve.TakeStep(idx, outpos);
// 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.
integralCurve.SetTime(idx, time);
inpos = outpos;
integralCurve.StepUpdate(idx, time, outpos);
tookAnySteps = true;
inpos = outpos;
}
// If the particle is at spatial or temporal boundary, take steps to just
// push it a little out of the boundary so that it will start advection in
// another domain, or in another time slice. Taking small steps enables
// reducing the error introduced at spatial or temporal boundaries.
else if (status == IntegratorStatus::OUTSIDE_TEMPORAL_BOUNDS)
//We can't take a step inside spatial boundary.
//Try and take a step just past the boundary.
else if (status.CheckSpatialBounds())
{
integralCurve.SetExitTemporalBoundary(idx);
break;
}
else if (status == IntegratorStatus::OUTSIDE_SPATIAL_BOUNDS)
{
status = integrator->SmallStep(inpos, time, outpos);
integralCurve.TakeStep(idx, outpos);
integralCurve.SetTime(idx, time);
if (status == IntegratorStatus::OUTSIDE_TEMPORAL_BOUNDS)
IntegratorStatus status2 = integrator->SmallStep(inpos, time, outpos);
if (status2.CheckOk())
{
integralCurve.SetExitTemporalBoundary(idx);
break;
}
else if (status == IntegratorStatus::OUTSIDE_SPATIAL_BOUNDS)
{
integralCurve.SetExitSpatialBoundary(idx);
break;
}
else if (status == IntegratorStatus::FAIL)
{
integralCurve.SetError(idx);
break;
integralCurve.StepUpdate(idx, time, outpos);
tookAnySteps = true;
//we took a step, so use this status to consider below.
status = status2;
}
}
}
integralCurve.SetTookAnySteps(idx, tookAnySteps);
integralCurve.StatusUpdate(idx, status, maxSteps);
} while (integralCurve.CanContinue(idx));
//Mark if any steps taken
integralCurve.UpdateTookSteps(idx, tookAnySteps);
//particle = integralCurve.GetParticle(idx);
//std::cout<<idx<<": "<<inpos<<" #"<<particle.NumSteps<<" "<<particle.Status<<std::endl;
}
};
@ -108,21 +111,19 @@ public:
~ParticleAdvectionWorklet() {}
void Run(const IntegratorType& integrator,
vtkm::cont::ArrayHandle<vtkm::Vec3f>& seedArray,
vtkm::Id maxSteps,
vtkm::cont::ArrayHandle<vtkm::Id>& statusArray,
vtkm::cont::ArrayHandle<vtkm::Id>& stepsTaken,
vtkm::cont::ArrayHandle<vtkm::FloatDefault>& timeArray)
vtkm::cont::ArrayHandle<vtkm::Particle>& particles,
vtkm::Id& MaxSteps)
{
using ParticleAdvectWorkletType = vtkm::worklet::particleadvection::ParticleAdvectWorklet;
using ParticleWorkletDispatchType =
typename vtkm::worklet::DispatcherMapField<ParticleAdvectWorkletType>;
using ParticleType = vtkm::worklet::particleadvection::Particles;
vtkm::Id numSeeds = static_cast<vtkm::Id>(seedArray.GetNumberOfValues());
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);
ParticleType particles(seedArray, stepsTaken, statusArray, timeArray, maxSteps);
#ifdef VTKM_CUDA
// This worklet needs some extra space on CUDA.
@ -130,123 +131,110 @@ public:
(void)stack;
#endif // VTKM_CUDA
ParticleType particlesObj(particles, MaxSteps);
//Invoke particle advection worklet
ParticleWorkletDispatchType particleWorkletDispatch;
particleWorkletDispatch.Invoke(idxArray, integrator, particles);
particleWorkletDispatch.Invoke(idxArray, integrator, particlesObj, maxSteps);
}
};
namespace detail
{
class Subtract : public vtkm::worklet::WorkletMapField
class GetSteps : public vtkm::worklet::WorkletMapField
{
public:
VTKM_CONT
Subtract() {}
using ControlSignature = void(FieldOut, FieldIn, FieldIn);
using ExecutionSignature = void(_1, _2, _3);
VTKM_EXEC void operator()(vtkm::Id& res, const vtkm::Id& x, const vtkm::Id& y) const
GetSteps() {}
using ControlSignature = void(FieldIn, FieldOut);
using ExecutionSignature = void(_1, _2);
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
template <typename IntegratorType>
class StreamlineWorklet
{
public:
template <typename PointStorage, typename FieldStorage>
template <typename PointStorage, typename PointStorage2>
void Run(const IntegratorType& it,
const vtkm::cont::ArrayHandle<vtkm::Vec3f, PointStorage>& pts,
const vtkm::Id& nSteps,
vtkm::cont::ArrayHandle<vtkm::Vec3f, PointStorage>& positions,
vtkm::cont::CellSetExplicit<>& polyLines,
vtkm::cont::ArrayHandle<vtkm::Id, FieldStorage>& statusArray,
vtkm::cont::ArrayHandle<vtkm::Id, FieldStorage>& stepsTaken,
vtkm::cont::ArrayHandle<vtkm::FloatDefault, FieldStorage>& timeArray)
vtkm::cont::ArrayHandle<vtkm::Particle, PointStorage>& particles,
vtkm::Id& MaxSteps,
vtkm::cont::ArrayHandle<vtkm::Vec3f, PointStorage2>& positions,
vtkm::cont::CellSetExplicit<>& polyLines)
{
integrator = it;
seedArray = pts;
maxSteps = nSteps;
run(positions, polyLines, statusArray, stepsTaken, timeArray);
}
struct IsOne
{
template <typename T>
VTKM_EXEC_CONT bool operator()(const T& x) const
{
return x == T(1);
}
};
private:
void run(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>& timeArray)
{
using ParticleWorkletDispatchType = typename vtkm::worklet::DispatcherMapField<
vtkm::worklet::particleadvection::ParticleAdvectWorklet>;
using StreamlineType = vtkm::worklet::particleadvection::StateRecordingParticles;
vtkm::cont::ArrayHandle<vtkm::Id> initialStepsTaken;
vtkm::Id numSeeds = static_cast<vtkm::Id>(particles.GetNumberOfValues());
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(4 * 1024);
vtkm::cont::cuda::ScopedCudaStackSize stack(16 * 1024);
(void)stack;
#endif // VTKM_CUDA
vtkm::cont::ArrayHandle<vtkm::Id> initialStepsTaken;
vtkm::cont::ArrayCopy(stepsTaken, initialStepsTaken);
vtkm::Id numSeeds = static_cast<vtkm::Id>(seedArray.GetNumberOfValues());
//Run streamline worklet
StreamlineType streamlines(particles, MaxSteps);
ParticleWorkletDispatchType particleWorkletDispatch;
vtkm::cont::ArrayHandleConstant<vtkm::Id> maxSteps(MaxSteps, numSeeds);
particleWorkletDispatch.Invoke(idxArray, it, streamlines, maxSteps);
vtkm::cont::ArrayHandleIndex idxArray(numSeeds);
//Get the positions
streamlines.GetCompactedHistory(positions);
vtkm::cont::ArrayHandle<vtkm::Id> validPoint;
std::vector<vtkm::Id> vpa(static_cast<std::size_t>(numSeeds * maxSteps), 0);
validPoint = vtkm::cont::make_ArrayHandle(vpa);
//Create the cells
vtkm::cont::ArrayHandle<vtkm::Id> numPoints;
vtkm::worklet::DispatcherMapField<detail::ComputeNumPoints> computeNumPointsDispatcher{ (
detail::ComputeNumPoints{}) };
computeNumPointsDispatcher.Invoke(particles, initialStepsTaken, numPoints);
//Compact history into positions.
vtkm::cont::ArrayHandle<vtkm::Vec3f> history;
StreamlineType streamlines(
seedArray, history, stepsTaken, status, timeArray, validPoint, maxSteps);
particleWorkletDispatch.Invoke(idxArray, integrator, streamlines);
vtkm::cont::Algorithm::CopyIf(history, validPoint, positions, IsOne());
vtkm::cont::ArrayHandle<vtkm::Id> stepsTakenNow;
stepsTakenNow.Allocate(numSeeds);
vtkm::worklet::DispatcherMapField<detail::Subtract> subtractDispatcher{ (detail::Subtract{}) };
subtractDispatcher.Invoke(stepsTakenNow, stepsTaken, initialStepsTaken);
//Create cells.
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::ArrayHandle<vtkm::Id> connectivity;
vtkm::cont::ArrayCopy(connCount, connectivity);
vtkm::cont::ArrayHandle<vtkm::UInt8> cellTypes;
cellTypes.Allocate(numSeeds);
vtkm::cont::ArrayHandleConstant<vtkm::UInt8> polyLineShape(vtkm::CELL_SHAPE_POLY_LINE,
numSeeds);
auto polyLineShape =
vtkm::cont::make_ArrayHandleConstant<vtkm::UInt8>(vtkm::CELL_SHAPE_POLY_LINE, numSeeds);
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);
polyLines.Fill(positions.GetNumberOfValues(), cellTypes, connectivity, offsets);
}
IntegratorType integrator;
vtkm::cont::ArrayHandle<vtkm::Vec3f> seedArray;
vtkm::Id maxSteps;
};
}
}

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

@ -11,11 +11,11 @@
#ifndef vtk_m_worklet_particleadvection_Particles_h
#define vtk_m_worklet_particleadvection_Particles_h
class ParticleExecutionObjectType;
#include <vtkm/Particle.h>
#include <vtkm/Types.h>
#include <vtkm/cont/DeviceAdapter.h>
#include <vtkm/cont/ExecutionObjectBase.h>
#include <vtkm/worklet/particleadvection/IntegratorStatus.h>
namespace vtkm
{
@ -23,171 +23,86 @@ namespace worklet
{
namespace particleadvection
{
enum class ParticleStatus
{
SUCCESS = 1,
TERMINATED = 1 << 1,
EXIT_SPATIAL_BOUNDARY = 1 << 2,
EXIT_TEMPORAL_BOUNDARY = 1 << 3,
FAIL = 1 << 4,
TOOK_ANY_STEPS = 1 << 5
};
template <typename Device>
class ParticleExecutionObject
{
public:
VTKM_EXEC_CONT
ParticleExecutionObject()
: Pos()
, Status()
, Steps()
, Time()
: Particles()
, MaxSteps(0)
{
}
ParticleExecutionObject(vtkm::cont::ArrayHandle<vtkm::Vec3f> posArray,
vtkm::cont::ArrayHandle<vtkm::Id> stepsArray,
vtkm::cont::ArrayHandle<vtkm::Id> statusArray,
vtkm::cont::ArrayHandle<vtkm::FloatDefault> timeArray,
vtkm::Id maxSteps)
ParticleExecutionObject(vtkm::cont::ArrayHandle<vtkm::Particle> particleArray, vtkm::Id maxSteps)
{
Pos = posArray.PrepareForInPlace(Device());
Steps = stepsArray.PrepareForInPlace(Device());
Status = statusArray.PrepareForInPlace(Device());
Time = timeArray.PrepareForInPlace(Device());
Particles = particleArray.PrepareForInPlace(Device());
MaxSteps = maxSteps;
}
VTKM_EXEC
void TakeStep(const vtkm::Id& idx, const vtkm::Vec3f& pt)
{
// Irrespective of what the advected status of the particle is,
// we need to set the output position as the last step taken by
// the particle, and increase the number of steps take by 1.
this->Pos.Set(idx, pt);
vtkm::Id nSteps = Steps.Get(idx);
this->Steps.Set(idx, ++nSteps);
vtkm::Particle GetParticle(const vtkm::Id& idx) { return this->Particles.Get(idx); }
// Check if the particle has completed the maximum steps required.
// If yes, set it to terminated.
if (nSteps == MaxSteps)
SetTerminated(idx);
VTKM_EXEC
void PreStepUpdate(const vtkm::Id& vtkmNotUsed(idx)) {}
VTKM_EXEC
void StepUpdate(const vtkm::Id& idx, vtkm::FloatDefault time, const vtkm::Vec3f& pt)
{
vtkm::Particle p = this->GetParticle(idx);
p.Pos = pt;
p.Time = time;
p.NumSteps++;
this->Particles.Set(idx, p);
}
/* Set/Change Status */
VTKM_EXEC
void SetOK(const vtkm::Id& idx)
void StatusUpdate(const vtkm::Id& idx,
const vtkm::worklet::particleadvection::IntegratorStatus& status,
vtkm::Id maxSteps)
{
Clear(idx);
Status.Set(idx, ParticleStatus::SUCCESS);
vtkm::Particle p = this->GetParticle(idx);
if (p.NumSteps == maxSteps)
p.Status.SetTerminate();
if (status.CheckFail())
p.Status.SetFail();
if (status.CheckSpatialBounds())
p.Status.SetSpatialBounds();
if (status.CheckTemporalBounds())
p.Status.SetTemporalBounds();
this->Particles.Set(idx, p);
}
VTKM_EXEC
void SetTerminated(const vtkm::Id& idx)
bool CanContinue(const vtkm::Id& idx)
{
ClearBit(idx, ParticleStatus::SUCCESS);
SetBit(idx, ParticleStatus::TERMINATED);
vtkm::Particle p = this->GetParticle(idx);
return (p.Status.CheckOk() && !p.Status.CheckTerminate() && !p.Status.CheckSpatialBounds() &&
!p.Status.CheckTemporalBounds());
}
VTKM_EXEC
void SetTookAnySteps(const vtkm::Id& idx, const bool& val)
void UpdateTookSteps(const vtkm::Id& idx, bool val)
{
vtkm::Particle p = this->GetParticle(idx);
if (val)
SetBit(idx, ParticleStatus::TOOK_ANY_STEPS);
p.Status.SetTookAnySteps();
else
ClearBit(idx, ParticleStatus::TOOK_ANY_STEPS);
p.Status.ClearTookAnySteps();
this->Particles.Set(idx, p);
}
VTKM_EXEC
void SetExitSpatialBoundary(const vtkm::Id& idx)
{
ClearBit(idx, ParticleStatus::SUCCESS);
SetBit(idx, ParticleStatus::EXIT_SPATIAL_BOUNDARY);
}
VTKM_EXEC
void SetExitTemporalBoundary(const vtkm::Id& idx)
{
ClearBit(idx, ParticleStatus::SUCCESS);
SetBit(idx, ParticleStatus::EXIT_TEMPORAL_BOUNDARY);
}
VTKM_EXEC
void SetError(const vtkm::Id& idx)
{
ClearBit(idx, ParticleStatus::SUCCESS);
SetBit(idx, ParticleStatus::FAIL);
}
/* Check Status */
VTKM_EXEC
bool OK(const vtkm::Id& idx) { return CheckBit(idx, ParticleStatus::SUCCESS); }
VTKM_EXEC
bool Terminated(const vtkm::Id& idx) { return CheckBit(idx, ParticleStatus::TERMINATED); }
VTKM_EXEC
bool ExitSpatialBoundary(const vtkm::Id& idx)
{
return CheckBit(idx, ParticleStatus::EXIT_SPATIAL_BOUNDARY);
}
VTKM_EXEC
bool ExitTemporalBoundary(const vtkm::Id& idx)
{
return CheckBit(idx, ParticleStatus::EXIT_TEMPORAL_BOUNDARY);
}
VTKM_EXEC
bool Error(const vtkm::Id& idx) { return CheckBit(idx, ParticleStatus::FAIL); }
VTKM_EXEC
bool Integrateable(const vtkm::Id& idx)
{
return OK(idx) && !(Terminated(idx) || ExitSpatialBoundary(idx) || ExitTemporalBoundary(idx));
}
VTKM_EXEC
bool Done(const vtkm::Id& idx) { return !Integrateable(idx); }
/* Bit Operations */
VTKM_EXEC
void Clear(const vtkm::Id& idx) { Status.Set(idx, 0); }
VTKM_EXEC
void SetBit(const vtkm::Id& idx, const ParticleStatus& b)
{
Status.Set(idx, Status.Get(idx) | static_cast<vtkm::Id>(b));
}
VTKM_EXEC
void ClearBit(const vtkm::Id& idx, const ParticleStatus& b)
{
Status.Set(idx, Status.Get(idx) & ~static_cast<vtkm::Id>(b));
}
VTKM_EXEC
bool CheckBit(const vtkm::Id& idx, const ParticleStatus& b) const
{
return (Status.Get(idx) & static_cast<vtkm::Id>(b)) != 0;
}
VTKM_EXEC
vtkm::Vec3f GetPos(const vtkm::Id& idx) const { return Pos.Get(idx); }
VTKM_EXEC
vtkm::Id GetStep(const vtkm::Id& idx) const { return Steps.Get(idx); }
VTKM_EXEC
vtkm::Id GetStatus(const vtkm::Id& idx) const { return Status.Get(idx); }
VTKM_EXEC
vtkm::FloatDefault GetTime(const vtkm::Id& idx) const { return Time.Get(idx); }
VTKM_EXEC
void SetTime(const vtkm::Id& idx, vtkm::FloatDefault time) const { Time.Set(idx, time); }
protected:
using IdPortal =
typename vtkm::cont::ArrayHandle<vtkm::Id>::template ExecutionTypes<Device>::Portal;
using PositionPortal =
typename vtkm::cont::ArrayHandle<vtkm::Vec3f>::template ExecutionTypes<Device>::Portal;
using FloatPortal =
typename vtkm::cont::ArrayHandle<vtkm::FloatDefault>::template ExecutionTypes<Device>::Portal;
using ParticlePortal =
typename vtkm::cont::ArrayHandle<vtkm::Particle>::template ExecutionTypes<Device>::Portal;
PositionPortal Pos;
IdPortal Status;
IdPortal Steps;
FloatPortal Time;
ParticlePortal Particles;
vtkm::Id MaxSteps;
};
class Particles : public vtkm::cont::ExecutionObjectBase
{
public:
@ -195,21 +110,13 @@ public:
VTKM_CONT vtkm::worklet::particleadvection::ParticleExecutionObject<Device> PrepareForExecution(
Device) const
{
return vtkm::worklet::particleadvection::ParticleExecutionObject<Device>(
this->PosArray, this->StepsArray, this->StatusArray, this->TimeArray, this->MaxSteps);
return vtkm::worklet::particleadvection::ParticleExecutionObject<Device>(this->ParticleArray,
this->MaxSteps);
}
VTKM_CONT
Particles(vtkm::cont::ArrayHandle<vtkm::Vec3f>& posArray,
vtkm::cont::ArrayHandle<vtkm::Id>& stepsArray,
vtkm::cont::ArrayHandle<vtkm::Id>& statusArray,
vtkm::cont::ArrayHandle<vtkm::FloatDefault>& timeArray,
const vtkm::Id& maxSteps)
: PosArray(posArray)
, StepsArray(stepsArray)
, StatusArray(statusArray)
, TimeArray(timeArray)
Particles(vtkm::cont::ArrayHandle<vtkm::Particle>& pArray, vtkm::Id& maxSteps)
: ParticleArray(pArray)
, MaxSteps(maxSteps)
{
}
@ -217,13 +124,7 @@ public:
Particles() {}
protected:
bool fromArray = false;
protected:
vtkm::cont::ArrayHandle<vtkm::Vec3f> PosArray;
vtkm::cont::ArrayHandle<vtkm::Id> StepsArray;
vtkm::cont::ArrayHandle<vtkm::Id> StatusArray;
vtkm::cont::ArrayHandle<vtkm::FloatDefault> TimeArray;
vtkm::cont::ArrayHandle<vtkm::Particle> ParticleArray;
vtkm::Id MaxSteps;
};
@ -237,93 +138,129 @@ public:
: ParticleExecutionObject<Device>()
, History()
, Length(0)
, StepCount()
, ValidPoint()
{
}
StateRecordingParticleExecutionObject(vtkm::cont::ArrayHandle<vtkm::Vec3f> posArray,
StateRecordingParticleExecutionObject(vtkm::cont::ArrayHandle<vtkm::Particle> pArray,
vtkm::cont::ArrayHandle<vtkm::Vec3f> historyArray,
vtkm::cont::ArrayHandle<vtkm::Id> stepsArray,
vtkm::cont::ArrayHandle<vtkm::Id> statusArray,
vtkm::cont::ArrayHandle<vtkm::FloatDefault> timeArray,
vtkm::cont::ArrayHandle<vtkm::Id> validPointArray,
vtkm::cont::ArrayHandle<vtkm::Id> stepCountArray,
vtkm::Id maxSteps)
: ParticleExecutionObject<Device>(posArray, stepsArray, statusArray, timeArray, maxSteps)
: ParticleExecutionObject<Device>(pArray, maxSteps)
, Length(maxSteps + 1)
{
Length = maxSteps;
vtkm::Id numPos = posArray.GetNumberOfValues();
vtkm::Id numPos = pArray.GetNumberOfValues();
History = historyArray.PrepareForOutput(numPos * Length, Device());
ValidPoint = validPointArray.PrepareForInPlace(Device());
StepCount = stepCountArray.PrepareForInPlace(Device());
}
VTKM_EXEC_CONT
void TakeStep(const vtkm::Id& idx, const vtkm::Vec3f& pt)
VTKM_EXEC
void PreStepUpdate(const vtkm::Id& idx)
{
this->ParticleExecutionObject<Device>::TakeStep(idx, pt);
vtkm::Particle p = this->ParticleExecutionObject<Device>::GetParticle(idx);
if (p.NumSteps == 0)
{
vtkm::Id loc = idx * Length;
this->History.Set(loc, p.Pos);
this->ValidPoint.Set(loc, 1);
this->StepCount.Set(idx, 1);
}
}
//TakeStep incremented the step, so we want the PREV step value.
vtkm::Id nSteps = this->Steps.Get(idx) - 1;
VTKM_EXEC
void StepUpdate(const vtkm::Id& idx, vtkm::FloatDefault time, const vtkm::Vec3f& pt)
{
this->ParticleExecutionObject<Device>::StepUpdate(idx, time, pt);
// Update the step for streamline storing portals.
// This includes updating the history and the valid points.
vtkm::Id loc = idx * Length + nSteps;
//local step count.
vtkm::Id stepCount = this->StepCount.Get(idx);
vtkm::Id loc = idx * Length + stepCount;
this->History.Set(loc, pt);
this->ValidPoint.Set(loc, 1);
}
vtkm::Vec3f GetHistory(const vtkm::Id& idx, const vtkm::Id& step) const
{
return this->History.Get(idx * this->Length + step);
this->StepCount.Set(idx, stepCount + 1);
}
protected:
using IdPortal =
typename vtkm::cont::ArrayHandle<vtkm::Id>::template ExecutionTypes<Device>::Portal;
using PositionPortal =
using HistoryPortal =
typename vtkm::cont::ArrayHandle<vtkm::Vec3f>::template ExecutionTypes<Device>::Portal;
PositionPortal History;
HistoryPortal History;
vtkm::Id Length;
IdPortal StepCount;
IdPortal ValidPoint;
};
class StateRecordingParticles : vtkm::cont::ExecutionObjectBase
{
public:
//Helper functor for compacting history
struct IsOne
{
template <typename T>
VTKM_EXEC_CONT bool operator()(const T& x) const
{
return x == T(1);
}
};
template <typename Device>
VTKM_CONT vtkm::worklet::particleadvection::StateRecordingParticleExecutionObject<Device>
PrepareForExecution(Device) const
{
return vtkm::worklet::particleadvection::StateRecordingParticleExecutionObject<Device>(
PosArray, HistoryArray, StepsArray, StatusArray, TimeArray, ValidPointArray, MaxSteps);
this->ParticleArray,
this->HistoryArray,
this->ValidPointArray,
this->StepCountArray,
this->MaxSteps);
}
VTKM_CONT
StateRecordingParticles(vtkm::cont::ArrayHandle<vtkm::Particle>& pArray, const vtkm::Id& maxSteps)
: MaxSteps(maxSteps)
, ParticleArray(pArray)
{
vtkm::Id numParticles = static_cast<vtkm::Id>(pArray.GetNumberOfValues());
//Create ValidPointArray initialized to zero.
vtkm::cont::ArrayHandleConstant<vtkm::Id> tmp(0, (this->MaxSteps + 1) * numParticles);
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
StateRecordingParticles(vtkm::cont::ArrayHandle<vtkm::Vec3f>& posArray,
StateRecordingParticles(vtkm::cont::ArrayHandle<vtkm::Particle>& pArray,
vtkm::cont::ArrayHandle<vtkm::Vec3f>& historyArray,
vtkm::cont::ArrayHandle<vtkm::Id>& stepsArray,
vtkm::cont::ArrayHandle<vtkm::Id>& statusArray,
vtkm::cont::ArrayHandle<vtkm::FloatDefault>& timeArray,
vtkm::cont::ArrayHandle<vtkm::Id>& validPointArray,
const vtkm::Id& maxSteps)
vtkm::Id& maxSteps)
{
PosArray = posArray;
ParticleArray = pArray;
HistoryArray = historyArray;
StepsArray = stepsArray;
StatusArray = statusArray;
TimeArray = timeArray;
ValidPointArray = validPointArray;
MaxSteps = maxSteps;
}
VTKM_CONT
void GetCompactedHistory(vtkm::cont::ArrayHandle<vtkm::Vec3f>& positions)
{
vtkm::cont::Algorithm::CopyIf(this->HistoryArray, this->ValidPointArray, positions, IsOne());
}
protected:
vtkm::cont::ArrayHandle<vtkm::Id> StepsArray;
vtkm::cont::ArrayHandle<vtkm::Id> StatusArray;
vtkm::cont::ArrayHandle<vtkm::FloatDefault> TimeArray;
vtkm::cont::ArrayHandle<vtkm::Id> ValidPointArray;
vtkm::cont::ArrayHandle<vtkm::Vec3f> HistoryArray;
vtkm::cont::ArrayHandle<vtkm::Vec3f> PosArray;
vtkm::Id MaxSteps;
vtkm::cont::ArrayHandle<vtkm::Particle> ParticleArray;
vtkm::cont::ArrayHandle<vtkm::Id> StepCountArray;
vtkm::cont::ArrayHandle<vtkm::Id> ValidPointArray;
};

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

@ -11,6 +11,7 @@
#ifndef 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>
namespace vtkm
@ -68,23 +69,34 @@ public:
}
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.
GridEvaluatorStatus status;
if (!(time >= TimeOne && time <= TimeTwo))
return EvaluatorStatus::OUTSIDE_TEMPORAL_BOUNDS;
EvaluatorStatus eval;
{
status.SetFail();
status.SetTemporalBounds();
return status;
}
Point one, two;
eval = this->EvaluatorOne.Evaluate(pos, one);
if (eval != EvaluatorStatus::SUCCESS)
return eval;
eval = this->EvaluatorTwo.Evaluate(pos, two);
if (eval != EvaluatorStatus::SUCCESS)
return eval;
status = this->EvaluatorOne.Evaluate(pos, one);
if (status.CheckFail())
return status;
status = this->EvaluatorTwo.Evaluate(pos, two);
if (status.CheckFail())
return status;
// LERP between the two values of calculated fields to obtain the new value
vtkm::FloatDefault proportion = (time - this->TimeOne) / this->TimeDiff;
out = vtkm::Lerp(one, two, proportion);
return EvaluatorStatus::SUCCESS;
status.SetOk();
return status;
}
private:

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

@ -23,6 +23,12 @@
#include <vtkm/io/writer/VTKDataSetWriter.h>
//timers
#include <chrono>
#include <ctime>
namespace
{
vtkm::FloatDefault vecData[125 * 3] = {
@ -285,27 +291,27 @@ public:
using ExecutionSignature = void(_1, _2, _3, _4);
template <typename EvaluatorType>
VTKM_EXEC void operator()(vtkm::Vec3f& pointIn,
VTKM_EXEC void operator()(vtkm::Particle& pointIn,
const EvaluatorType& evaluator,
vtkm::worklet::particleadvection::EvaluatorStatus& status,
vtkm::worklet::particleadvection::GridEvaluatorStatus& status,
vtkm::Vec3f& pointOut) const
{
status = evaluator.Evaluate(pointIn, pointOut);
status = evaluator.Evaluate(pointIn.Pos, pointOut);
}
};
template <typename EvalType>
void ValidateEvaluator(const EvalType& eval,
const std::vector<vtkm::Vec3f>& pointIns,
const std::vector<vtkm::Particle>& pointIns,
const vtkm::Vec3f& vec,
const std::string& msg)
{
using EvalTester = TestEvaluatorWorklet;
using EvalTesterDispatcher = vtkm::worklet::DispatcherMapField<EvalTester>;
using Status = vtkm::worklet::particleadvection::EvaluatorStatus;
using Status = vtkm::worklet::particleadvection::GridEvaluatorStatus;
EvalTester evalTester;
EvalTesterDispatcher evalTesterDispatcher(evalTester);
vtkm::cont::ArrayHandle<vtkm::Vec3f> pointsHandle = vtkm::cont::make_ArrayHandle(pointIns);
vtkm::cont::ArrayHandle<vtkm::Particle> pointsHandle = vtkm::cont::make_ArrayHandle(pointIns);
vtkm::Id numPoints = pointsHandle.GetNumberOfValues();
vtkm::cont::ArrayHandle<Status> evalStatus;
vtkm::cont::ArrayHandle<vtkm::Vec3f> evalResults;
@ -316,8 +322,8 @@ void ValidateEvaluator(const EvalType& eval,
{
Status status = statusPortal.Get(index);
vtkm::Vec3f result = resultsPortal.Get(index);
VTKM_TEST_ASSERT(status == Status::SUCCESS, "Error in evaluator for " + msg);
VTKM_TEST_ASSERT(test_equal(result, vec), "Error in evaluator result for " + msg);
VTKM_TEST_ASSERT(status.CheckOk(), "Error in evaluator for " + msg);
VTKM_TEST_ASSERT(result == vec, "Error in evaluator result for " + msg);
}
pointsHandle.ReleaseResources();
evalStatus.ReleaseResources();
@ -335,20 +341,22 @@ public:
using ExecutionSignature = void(_1, _2, _3, _4);
template <typename IntegratorType>
VTKM_EXEC void operator()(vtkm::Vec3f& pointIn,
VTKM_EXEC void operator()(vtkm::Particle& pointIn,
const IntegratorType* integrator,
vtkm::worklet::particleadvection::IntegratorStatus& status,
vtkm::Vec3f& pointOut) const
{
vtkm::FloatDefault time = 0;
status = integrator->Step(pointIn, time, pointOut);
status = integrator->Step(pointIn.Pos, time, pointOut);
if (status.CheckSpatialBounds())
status = integrator->SmallStep(pointIn.Pos, time, pointOut);
}
};
template <typename IntegratorType>
void ValidateIntegrator(const IntegratorType& integrator,
const std::vector<vtkm::Vec3f>& pointIns,
const std::vector<vtkm::Particle>& pointIns,
const std::vector<vtkm::Vec3f>& expStepResults,
const std::string& msg)
{
@ -356,7 +364,7 @@ void ValidateIntegrator(const IntegratorType& integrator,
using IntegratorTesterDispatcher = vtkm::worklet::DispatcherMapField<IntegratorTester>;
using Status = vtkm::worklet::particleadvection::IntegratorStatus;
IntegratorTesterDispatcher integratorTesterDispatcher;
vtkm::cont::ArrayHandle<vtkm::Vec3f> pointsHandle = vtkm::cont::make_ArrayHandle(pointIns);
auto pointsHandle = vtkm::cont::make_ArrayHandle(pointIns);
vtkm::Id numPoints = pointsHandle.GetNumberOfValues();
vtkm::cont::ArrayHandle<Status> stepStatus;
vtkm::cont::ArrayHandle<vtkm::Vec3f> stepResults;
@ -368,12 +376,12 @@ void ValidateIntegrator(const IntegratorType& integrator,
{
Status status = statusPortal.Get(index);
vtkm::Vec3f result = resultsPortal.Get(index);
VTKM_TEST_ASSERT(status != Status::FAIL, "Error in evaluator for " + msg);
if (status == Status::OUTSIDE_SPATIAL_BOUNDS)
VTKM_TEST_ASSERT(test_equal(result, pointsPortal.Get(index)),
VTKM_TEST_ASSERT(status.CheckOk(), "Error in evaluator for " + msg);
if (status.CheckSpatialBounds())
VTKM_TEST_ASSERT(result == pointsPortal.Get(index).Pos,
"Error in evaluator result for [OUTSIDE SPATIAL]" + msg);
else
VTKM_TEST_ASSERT(test_equal(result, expStepResults[(size_t)index]),
VTKM_TEST_ASSERT(result == expStepResults[(size_t)index],
"Error in evaluator result for " + msg);
}
pointsHandle.ReleaseResources();
@ -382,19 +390,17 @@ void ValidateIntegrator(const IntegratorType& integrator,
}
template <typename IntegratorType>
void ValidateIntegratorForBoundary(const vtkm::Vec3f& vector,
const vtkm::Bounds& bounds,
void ValidateIntegratorForBoundary(const vtkm::Bounds& bounds,
const IntegratorType& integrator,
const std::vector<vtkm::Vec3f>& pointIns,
const std::vector<vtkm::Particle>& pointIns,
const std::string& msg)
{
using IntegratorTester = TestIntegratorWorklet;
using IntegratorTesterDispatcher = vtkm::worklet::DispatcherMapField<IntegratorTester>;
using Status = vtkm::worklet::particleadvection::IntegratorStatus;
vtkm::FloatDefault tolerance = vtkm::Epsilon<vtkm::FloatDefault>() * 100;
IntegratorTesterDispatcher integratorTesterDispatcher;
vtkm::cont::ArrayHandle<vtkm::Vec3f> pointsHandle = vtkm::cont::make_ArrayHandle(pointIns);
auto pointsHandle = vtkm::cont::make_ArrayHandle(pointIns);
vtkm::Id numPoints = pointsHandle.GetNumberOfValues();
vtkm::cont::ArrayHandle<Status> stepStatus;
vtkm::cont::ArrayHandle<vtkm::Vec3f> stepResults;
@ -404,15 +410,18 @@ void ValidateIntegratorForBoundary(const vtkm::Vec3f& vector,
for (vtkm::Id index = 0; index < numPoints; index++)
{
Status status = statusPortal.Get(index);
VTKM_TEST_ASSERT(status.CheckSpatialBounds(), "Error in evaluator for " + msg);
vtkm::Vec3f result = resultsPortal.Get(index);
if (vector[0] == 1.)
VTKM_TEST_ASSERT((bounds.X.Max - result[0]) < tolerance, "X Tolerance not satisfied.");
if (vector[1] == 1.)
VTKM_TEST_ASSERT((bounds.Y.Max - result[1]) < tolerance, "Y Tolerance not satisfied.");
if (vector[2] == 1.)
VTKM_TEST_ASSERT((bounds.Z.Max - result[2]) < tolerance, "Z Tolerance not satisfied.");
VTKM_TEST_ASSERT(status == Status::OUTSIDE_SPATIAL_BOUNDS, "Error in evaluator for " + msg);
if (bounds.Contains(result))
{
std::cout << "Failure. " << bounds << std::endl;
std::cout << std::setprecision(12) << index << ": " << bounds << " res= " << result
<< std::endl;
}
//VTKM_TEST_ASSERT(!bounds.Contains(result), "Tolerance not satisfied.");
}
pointsHandle.ReleaseResources();
stepStatus.ReleaseResources();
stepResults.ReleaseResources();
@ -425,18 +434,19 @@ void TestEvaluators()
using RK4Type = vtkm::worklet::particleadvection::RK4Integrator<GridEvalType>;
std::vector<vtkm::Vec3f> vecs;
vecs.push_back(vtkm::Vec3f(1, 0, 0));
vecs.push_back(vtkm::Vec3f(0, 1, 0));
vecs.push_back(vtkm::Vec3f(0, 0, 1));
vecs.push_back(vtkm::Vec3f(1, 1, 0));
vecs.push_back(vtkm::Vec3f(0, 1, 1));
vecs.push_back(vtkm::Vec3f(1, 0, 1));
vecs.push_back(vtkm::Vec3f(1, 1, 1));
vtkm::FloatDefault vals[3] = { -1., 0., 1. };
for (int i = 0; i < 3; i++)
for (int j = 0; j < 3; j++)
for (int k = 0; k < 3; k++)
if (!(i == 1 && j == 1 && k == 1)) //don't add a [0,0,0] vec.
vecs.push_back(vtkm::Vec3f(vals[i], vals[j], vals[k]));
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));
bounds.push_back(vtkm::Bounds(0, 1000, 0, 1, -1, 1000));
bounds.push_back(vtkm::Bounds(0, 1000, -100, 0, -1, 1000));
std::vector<vtkm::Id3> dims;
dims.push_back(vtkm::Id3(5, 5, 5));
@ -456,8 +466,9 @@ void TestEvaluators()
vtkm::cont::ArrayHandle<vtkm::Vec3f> vecField;
CreateConstantVectorField(dim[0] * dim[1] * dim[2], vec, vecField);
vtkm::FloatDefault stepSize = 0.01f;
std::vector<vtkm::Vec3f> pointIns;
//vtkm::FloatDefault stepSize = 0.01f;
vtkm::FloatDefault stepSize = 0.1f;
std::vector<vtkm::Particle> pointIns;
std::vector<vtkm::Vec3f> stepResult;
//Create a bunch of random points in the bounds.
srand(314);
@ -472,23 +483,34 @@ void TestEvaluators()
for (int k = 0; k < 38; k++)
{
auto p = RandomPoint(interiorBounds);
pointIns.push_back(p);
pointIns.push_back(vtkm::Particle(p, k));
stepResult.push_back(p + vec * stepSize);
}
std::vector<vtkm::Vec3f> boundaryPoints;
const vtkm::Range xRange(bound.X.Max - stepSize / 2., bound.X.Max);
const vtkm::Range yRange(bound.Y.Max - stepSize / 2., bound.Y.Max);
const vtkm::Range zRange(bound.Z.Max - stepSize / 2., bound.Z.Max);
vtkm::Range xRange, yRange, zRange;
if (vec[0] > 0)
xRange = vtkm::Range(bound.X.Max - stepSize / 2., bound.X.Max);
else
xRange = vtkm::Range(bound.X.Min, bound.X.Min + stepSize / 2.);
if (vec[1] > 0)
yRange = vtkm::Range(bound.Y.Max - stepSize / 2., bound.Y.Max);
else
yRange = vtkm::Range(bound.Y.Min, bound.Y.Min + stepSize / 2.);
if (vec[2] > 0)
zRange = vtkm::Range(bound.Z.Max - stepSize / 2., bound.Z.Max);
else
zRange = vtkm::Range(bound.Z.Min, bound.Z.Min + stepSize / 2.);
vtkm::Bounds forBoundary(xRange, yRange, zRange);
// Generate a bunch of boundary points towards the face of the direction
// of the velocity field
// All velocities are in the +ve direction.
std::vector<vtkm::Particle> boundaryPoints;
for (int k = 0; k < 10; k++)
{
// Generate bunch of boundary points towards the face of the direction
// of the velocity field
// All velocities are in the +ve direction.
auto p = RandomPoint(forBoundary);
pointIns.push_back(p);
}
boundaryPoints.push_back(vtkm::Particle(RandomPoint(forBoundary), k));
for (auto& ds : dataSets)
{
@ -498,7 +520,7 @@ void TestEvaluators()
RK4Type rk4(gridEval, stepSize);
ValidateIntegrator(rk4, pointIns, stepResult, "constant vector RK4");
ValidateIntegratorForBoundary(vec, bound, rk4, boundaryPoints, "constant vector RK4");
ValidateIntegratorForBoundary(bound, rk4, boundaryPoints, "constant vector RK4");
}
}
}
@ -509,26 +531,36 @@ void ValidateParticleAdvectionResult(const vtkm::worklet::ParticleAdvectionResul
vtkm::Id nSeeds,
vtkm::Id maxSteps)
{
VTKM_TEST_ASSERT(res.positions.GetNumberOfValues() == nSeeds,
VTKM_TEST_ASSERT(res.Particles.GetNumberOfValues() == nSeeds,
"Number of output particles does not match input.");
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 particle advection");
VTKM_TEST_ASSERT(res.Particles.GetPortalConstControl().Get(i).Status.CheckOk(),
"Bad status in particle advection");
}
}
void ValidateStreamlineResult(const vtkm::worklet::StreamlineResult& res,
vtkm::Id nSeeds,
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.");
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");
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 GridEvalType = vtkm::worklet::particleadvection::GridEvaluator<FieldHandle>;
@ -568,13 +600,16 @@ void TestParticleWorklets()
dataSets.push_back(CreateWeirdnessFromStructuredDataSet(dataSets[0], DataSetOption::EXPLICIT));
//Generate three random points.
std::vector<vtkm::Vec3f> pts;
pts.push_back(RandomPoint(bound));
pts.push_back(RandomPoint(bound));
pts.push_back(RandomPoint(bound));
std::vector<vtkm::Particle> pts;
pts.push_back(vtkm::Particle(RandomPoint(bound), 0));
pts.push_back(vtkm::Particle(RandomPoint(bound), 1));
pts.push_back(vtkm::Particle(RandomPoint(bound), 2));
std::vector<vtkm::Particle> pts2 = pts;
vtkm::Id nSeeds = static_cast<vtkm::Id>(pts.size());
std::vector<vtkm::Id> stepsTaken = { 10, 20, 600 };
for (std::size_t i = 0; i < stepsTaken.size(); i++)
pts2[i].NumSteps = stepsTaken[i];
for (auto& ds : dataSets)
{
@ -586,17 +621,20 @@ void TestParticleWorklets()
//Streamline worklet with and without steps taken.
for (int i = 0; i < 4; i++)
{
auto seedsArray = vtkm::cont::make_ArrayHandle(pts, vtkm::CopyFlag::On);
auto stepsTakenArray = vtkm::cont::make_ArrayHandle(stepsTaken, vtkm::CopyFlag::On);
if (i < 2)
{
vtkm::worklet::ParticleAdvection pa;
vtkm::worklet::ParticleAdvectionResult res;
if (i == 0)
res = pa.Run(rk4, seedsArray, maxSteps);
{
auto seeds = vtkm::cont::make_ArrayHandle(pts, vtkm::CopyFlag::On);
res = pa.Run(rk4, seeds, maxSteps);
}
else
res = pa.Run(rk4, seedsArray, stepsTakenArray, maxSteps);
{
auto seeds = vtkm::cont::make_ArrayHandle(pts2, vtkm::CopyFlag::On);
res = pa.Run(rk4, seeds, maxSteps);
}
ValidateParticleAdvectionResult(res, nSeeds, maxSteps);
}
else
@ -604,9 +642,15 @@ void TestParticleWorklets()
vtkm::worklet::Streamline s;
vtkm::worklet::StreamlineResult res;
if (i == 2)
res = s.Run(rk4, seedsArray, maxSteps);
{
auto seeds = vtkm::cont::make_ArrayHandle(pts, vtkm::CopyFlag::On);
res = s.Run(rk4, seeds, maxSteps);
}
else
res = s.Run(rk4, seedsArray, stepsTakenArray, maxSteps);
{
auto seeds = vtkm::cont::make_ArrayHandle(pts2, vtkm::CopyFlag::On);
res = s.Run(rk4, seeds, maxSteps);
}
ValidateStreamlineResult(res, nSeeds, maxSteps);
}
}
@ -639,26 +683,152 @@ void TestParticleStatus()
RK4Type rk4(eval, stepSize);
vtkm::worklet::ParticleAdvection pa;
std::vector<vtkm::Vec3f> pts;
pts.push_back(vtkm::Vec3f(.5, .5, .5));
pts.push_back(vtkm::Vec3f(-1, -1, -1));
std::vector<vtkm::Particle> pts;
pts.push_back(vtkm::Particle(vtkm::Vec3f(.5, .5, .5), 0));
pts.push_back(vtkm::Particle(vtkm::Vec3f(-1, -1, -1), 1));
auto seedsArray = vtkm::cont::make_ArrayHandle(pts, vtkm::CopyFlag::On);
auto res = pa.Run(rk4, seedsArray, maxSteps);
auto statusPortal = res.status.GetPortalConstControl();
pa.Run(rk4, seedsArray, maxSteps);
auto portal = seedsArray.GetPortalConstControl();
vtkm::Id tookStep0 = statusPortal.Get(0) &
static_cast<vtkm::Id>(vtkm::worklet::particleadvection::ParticleStatus::TOOK_ANY_STEPS);
vtkm::Id tookStep1 = statusPortal.Get(1) &
static_cast<vtkm::Id>(vtkm::worklet::particleadvection::ParticleStatus::TOOK_ANY_STEPS);
VTKM_TEST_ASSERT(tookStep0 != 0, "Particle failed to take any steps");
VTKM_TEST_ASSERT(tookStep1 == 0, "Particle took a step when it should not have.");
bool tookStep0 = portal.Get(0).Status.CheckTookAnySteps();
bool tookStep1 = portal.Get(1).Status.CheckTookAnySteps();
VTKM_TEST_ASSERT(tookStep0 == true, "Particle failed to take any steps");
VTKM_TEST_ASSERT(tookStep1 == false, "Particle took a step when it should not have.");
}
void TestWorkletsBasic()
{
using FieldHandle = vtkm::cont::ArrayHandle<vtkm::Vec3f>;
using GridEvalType = vtkm::worklet::particleadvection::GridEvaluator<FieldHandle>;
using RK4Type = vtkm::worklet::particleadvection::RK4Integrator<GridEvalType>;
vtkm::FloatDefault stepSize = 0.01f;
const vtkm::Id3 dims(5, 5, 5);
vtkm::Id nElements = dims[0] * dims[1] * dims[2] * 3;
std::vector<vtkm::Vec3f> field;
vtkm::Vec3f vecDir(1, 0, 0);
for (vtkm::Id i = 0; i < nElements; i++)
field.push_back(vtkm::Normal(vecDir));
vtkm::cont::ArrayHandle<vtkm::Vec3f> fieldArray;
fieldArray = vtkm::cont::make_ArrayHandle(field);
vtkm::Bounds bounds(0, 1, 0, 1, 0, 1);
auto ds = CreateUniformDataSet(bounds, dims);
GridEvalType eval(ds.GetCoordinateSystem(), ds.GetCellSet(), fieldArray);
RK4Type rk4(eval, stepSize);
vtkm::Id maxSteps = 83;
std::vector<std::string> workletTypes = { "particleAdvection", "streamline" };
vtkm::FloatDefault endT = stepSize * static_cast<vtkm::FloatDefault>(maxSteps);
for (auto w : workletTypes)
{
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);
for (int i = 0; i < 8; i++)
{
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[static_cast<std::size_t>(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[static_cast<std::size_t>(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[static_cast<std::size_t>(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()
{
TestEvaluators();
TestParticleWorklets();
TestParticleStatus();
TestWorkletsBasic();
TestParticleWorkletsWithDataSetTypes();
}
int UnitTestParticleAdvection(int argc, char* argv[])

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

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