Shuenhoy/vtk-m
1
0
Fork 0
mirror of https://gitlab.kitware.com/vtk/vtk-m synced 2024-10-08 11:29:02 +00:00
Code Issues 2 Actions Packages Projects Releases Wiki Activity

Merge topic 'flow_restrucutre'

Browse Source 
67b7543a3 Adding documentation for flow filter restructure
dbc873efa Changes to address feedback from MR
67716402b Correct export in class declaration
6d1d4f90a Fixing linking issues for flow Analysis class
adcb42455 Removing unnecessary file
78ca3f301 Fixing linking issues for flow Analysis class
0e1ade83a Fixing linking issues for flow Analysis class
12a3bc94e Adding test dependency of filter_flow on tests
...

Acked-by: Kitware Robot <kwrobot@kitware.com>
Merge-request: !3087
This commit is contained in:
Abhishek Yenpure 2023-08-18 02:17:20 +00:00 committed by Kitware Robot
commit 969e55048b
42 changed files with 1946 additions and 1413 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

48
docs/changelog/make-flow-filters-modular-and-extensible-using-traits.md Normal file
View File

@ -0,0 +1,48 @@
## Make flow filters modular and extensible using traits
Many flow filters have common underpinnings in term of the components they use.
E.g., the choice and handling for solvers, analysis, termination, vector field, etc.
However, having these components baked hard in the infrastructure makes extensibility chanllenging,
which leads to developers implementing bespoke solutions.
This change establishes an infrastructure for easy specification and development of flow filter.
To that end, two new abstractions are introduced along with thier basic implementations : `Analysis` and `Termination`
* `Analysis` defines how each step of the particle needs to be analyzed
* `Termination` defines the termination criteria for every particle
The two, in addition to the existing abstractions for `Particle` and `Field` can be used to specify
novel flow filters. This is accomplished by defining a new trait for the new filter using implementations
for these abstractions.
E.g., for specifying the streamline filter for a general case the following trait can be used
```cpp
template <>
struct FlowTraits<Streamline>
{
using ParticleType = vtkm::Particle;
using TerminationType = vtkm::worklet::flow::NormalTermination;
using AnalysisType = vtkm::worklet::flow::StreamlineAnalysis<ParticleType>;
using ArrayType = vtkm::cont::ArrayHandle<vtkm::Vec3f>;
using FieldType = vtkm::worklet::flow::VelocityField<ArrayType>;
};
```
Similarly, to produce a flow map, the following trait can be used
```cpp
template <>
struct FlowTraits<ParticleAdvection>
{
using ParticleType = vtkm::Particle;
using TerminationType = vtkm::worklet::flow::NormalTermination;
using AnalysisType = vtkm::worklet::flow::NoAnalysis<ParticleType>;
using ArrayType = vtkm::cont::ArrayHandle<vtkm::Vec3f>;
using FieldType = vtkm::worklet::flow::VelocityField<ArrayType>;
};
```
These traits are enough for the infrastrucutre to use the correct code paths to produce the desired
result.
Along with changing the existing filters to use this new way of specification of components,
a new filter `WarpXStreamline` has been added to enable streamline analysis for charged particles for
the WarpX simulation.

34
vtkm/Particle.h
View File

@ -143,7 +143,7 @@ public:
VTKM_EXEC_CONT
vtkm::Vec3f Velocity(const vtkm::VecVariable<vtkm::Vec3f, 2>& vectors,
const vtkm::FloatDefault& vtkmNotUsed(length))
const vtkm::FloatDefault& vtkmNotUsed(length)) const
{
// Velocity is evaluated from the Velocity field
// and is not influenced by the particle
@ -213,6 +213,29 @@ public:
{
}
VTKM_EXEC_CONT
ChargedParticle(const vtkm::ChargedParticle& other)
: Position(other.Position)
, ID(other.ID)
, NumSteps(other.NumSteps)
, Status(other.Status)
, Time(other.Time)
, Mass(other.Mass)
, Charge(other.Charge)
, Weighting(other.Weighting)
, Momentum(other.Momentum)
{
}
vtkm::ChargedParticle& operator=(const vtkm::ChargedParticle&) = default;
VTKM_EXEC_CONT
~ChargedParticle() noexcept
{
// This must not be defaulted, since defaulted virtual destructors are
// troublesome with CUDA __host__ __device__ markup.
}
VTKM_EXEC_CONT const vtkm::Vec3f& GetPosition() const { return this->Position; }
VTKM_EXEC_CONT void SetPosition(const vtkm::Vec3f& position) { this->Position = position; }
@ -230,7 +253,7 @@ public:
VTKM_EXEC_CONT void SetTime(vtkm::FloatDefault time) { this->Time = time; }
VTKM_EXEC_CONT
vtkm::Float64 Gamma(vtkm::Vec3f momentum, bool reciprocal = false) const
vtkm::Float64 Gamma(const vtkm::Vec3f& momentum, bool reciprocal = false) const
{
constexpr vtkm::FloatDefault c2 = SPEED_OF_LIGHT * SPEED_OF_LIGHT;
const vtkm::Float64 fMom2 = vtkm::MagnitudeSquared(momentum);
@ -244,7 +267,7 @@ public:
VTKM_EXEC_CONT
vtkm::Vec3f Velocity(const vtkm::VecVariable<vtkm::Vec3f, 2>& vectors,
const vtkm::FloatDefault& length)
const vtkm::FloatDefault& length) const
{
VTKM_ASSERT(vectors.GetNumberOfComponents() == 2);
@ -265,8 +288,6 @@ public:
const vtkm::Vec3f mom_plus = mom_minus + vtkm::Cross(mom_prime, s);
const vtkm::Vec3f mom_new = mom_plus + 0.5 * this->Charge * eField * length;
//TODO : Sould this be a const method?
// If yes, need a better way to update momentum
this->Momentum = mom_new;
// momentum = velocity * mass * gamma;
@ -302,7 +323,7 @@ private:
vtkm::Float64 Mass;
vtkm::Float64 Charge;
vtkm::Float64 Weighting;
vtkm::Vec3f Momentum;
mutable vtkm::Vec3f Momentum;
constexpr static vtkm::FloatDefault SPEED_OF_LIGHT =
static_cast<vtkm::FloatDefault>(2.99792458e8);
@ -327,7 +348,6 @@ public:
} //namespace vtkm
namespace mangled_diy_namespace
{
template <>

13
vtkm/filter/flow/CMakeLists.txt
View File

@ -20,23 +20,28 @@ set(flow_headers
PathParticle.h
Streamline.h
StreamSurface.h
WarpXStreamline.h
)
set(flow_sources
internal/Messenger.cxx
FilterParticleAdvection.cxx
ParticleAdvection.cxx
Pathline.cxx
PathParticle.cxx
Streamline.cxx
)
set(flow_device_sources
worklet/Analysis.cxx
Lagrangian.cxx
LagrangianStructures.cxx
FilterParticleAdvectionSteadyState.cxx
FilterParticleAdvectionUnsteadyState.cxx
StreamSurface.cxx
Lagrangian.cxx
LagrangianStructures.cxx
ParticleAdvection.cxx
Pathline.cxx
PathParticle.cxx
Streamline.cxx
WarpXStreamline.cxx
)
vtkm_library(

21
vtkm/filter/flow/FilterParticleAdvection.h
View File

@ -12,6 +12,7 @@
#define vtk_m_filter_flow_FilterParticleAdvection_h
#include <vtkm/Particle.h>
#include <vtkm/cont/ErrorFilterExecution.h>
#include <vtkm/filter/FilterField.h>
#include <vtkm/filter/flow/FlowTypes.h>
#include <vtkm/filter/flow/vtkm_filter_flow_export.h>
@ -68,23 +69,6 @@ public:
this->SolverType = vtkm::filter::flow::IntegrationSolverType::EULER_TYPE;
}
VTKM_CONT
void SetVectorFieldType(vtkm::filter::flow::VectorFieldType vecFieldType)
{
this->VecFieldType = vecFieldType;
}
///@{
/// Choose the field to operate on. Note, if
/// `this->UseCoordinateSystemAsField` is true, then the active field is not used.
VTKM_CONT void SetEField(const std::string& name) { this->SetActiveField(0, name); }
VTKM_CONT void SetBField(const std::string& name) { this->SetActiveField(1, name); }
VTKM_CONT std::string GetEField() const { return this->GetActiveFieldName(0); }
VTKM_CONT std::string GetBField() const { return this->GetActiveFieldName(1); }
VTKM_CONT
bool GetUseThreadedAlgorithm() { return this->UseThreadedAlgorithm; }
@ -104,10 +88,9 @@ public:
protected:
VTKM_CONT virtual void ValidateOptions() const;
VTKM_CONT virtual vtkm::filter::flow::FlowResultType GetResultType() const = 0;
bool BlockIdsSet = false;
std::vector<vtkm::Id> BlockIds;
vtkm::Id NumberOfSteps = 0;
vtkm::cont::UnknownArrayHandle Seeds;
vtkm::filter::flow::IntegrationSolverType SolverType =

144
vtkm/filter/flow/FilterParticleAdvectionSteadyState.cxx
View File

@ -9,9 +9,9 @@
//============================================================================
#include <vtkm/filter/flow/FilterParticleAdvectionSteadyState.h>
#include <vtkm/filter/flow/internal/DataSetIntegratorSteadyState.h>
#include <vtkm/filter/flow/internal/BoundsMap.h>
#include <vtkm/filter/flow/internal/DataSetIntegratorSteadyState.h>
#include <vtkm/filter/flow/internal/ParticleAdvector.h>
namespace vtkm
@ -21,87 +21,89 @@ namespace filter
namespace flow
{
namespace
template <typename Derived>
VTKM_CONT typename FilterParticleAdvectionSteadyState<Derived>::FieldType
FilterParticleAdvectionSteadyState<Derived>::GetField(const vtkm::cont::DataSet& data) const
{
VTKM_CONT std::vector<vtkm::filter::flow::internal::DataSetIntegratorSteadyState>
CreateDataSetIntegrators(
const vtkm::cont::PartitionedDataSet& input,
const vtkm::cont::Variant<std::string, std::pair<std::string, std::string>>& activeField,
const vtkm::filter::flow::internal::BoundsMap& boundsMap,
const vtkm::filter::flow::IntegrationSolverType solverType,
const vtkm::filter::flow::VectorFieldType vecFieldType,
const vtkm::filter::flow::FlowResultType resultType)
{
using DSIType = vtkm::filter::flow::internal::DataSetIntegratorSteadyState;
const Derived* inst = static_cast<const Derived*>(this);
return inst->GetField(data);
}
template <typename Derived>
VTKM_CONT typename FilterParticleAdvectionSteadyState<Derived>::TerminationType
FilterParticleAdvectionSteadyState<Derived>::GetTermination(const vtkm::cont::DataSet& data) const
{
const Derived* inst = static_cast<const Derived*>(this);
return inst->GetTermination(data);
}
template <typename Derived>
VTKM_CONT typename FilterParticleAdvectionSteadyState<Derived>::AnalysisType
FilterParticleAdvectionSteadyState<Derived>::GetAnalysis(const vtkm::cont::DataSet& data) const
{
const Derived* inst = static_cast<const Derived*>(this);
return inst->GetAnalysis(data);
}
template <typename Derived>
VTKM_CONT vtkm::cont::PartitionedDataSet
FilterParticleAdvectionSteadyState<Derived>::DoExecutePartitions(
const vtkm::cont::PartitionedDataSet& input)
{
//using ParticleType = FilterParticleAdvectionSteadyState<Derived>::ParticleType;
//using FieldType = FilterParticleAdvectionSteadyState<Derived>::FieldType;
//using TerminationType = FilterParticleAdvectionSteadyState<Derived>::TerminationType;
//using AnalysisType = FilterParticleAdvectionSteadyState<Derived>::AnalysisType;
using DSIType = vtkm::filter::flow::internal::
DataSetIntegratorSteadyState<ParticleType, FieldType, TerminationType, AnalysisType>;
this->ValidateOptions();
vtkm::filter::flow::internal::BoundsMap boundsMap(input);
std::vector<DSIType> dsi;
for (vtkm::Id i = 0; i < input.GetNumberOfPartitions(); i++)
{
vtkm::Id blockId = boundsMap.GetLocalBlockId(i);
auto ds = input.GetPartition(i);
if (activeField.IsType<DSIType::VelocityFieldNameType>())
{
const auto& fieldNm = activeField.Get<DSIType::VelocityFieldNameType>();
if (!ds.HasPointField(fieldNm) && !ds.HasCellField(fieldNm))
throw vtkm::cont::ErrorFilterExecution("Unsupported field assocation");
}
else if (activeField.IsType<DSIType::ElectroMagneticFieldNameType>())
{
const auto& fieldNm = activeField.Get<DSIType::ElectroMagneticFieldNameType>();
const auto& electric = fieldNm.first;
const auto& magnetic = fieldNm.second;
if (!ds.HasPointField(electric) && !ds.HasCellField(electric))
throw vtkm::cont::ErrorFilterExecution("Unsupported field assocation");
if (!ds.HasPointField(magnetic) && !ds.HasCellField(magnetic))
throw vtkm::cont::ErrorFilterExecution("Unsupported field assocation");
}
dsi.emplace_back(ds, blockId, activeField, solverType, vecFieldType, resultType);
}
return dsi;
auto dataset = input.GetPartition(i);
// Build the field for the current dataset
FieldType field = this->GetField(dataset);
// Build the termination for the current dataset
TerminationType termination = this->GetTermination(dataset);
// Build the analysis for the current dataset
AnalysisType analysis = this->GetAnalysis(dataset);
dsi.emplace_back(blockId, field, dataset, this->SolverType, termination, analysis);
}
} //anonymous namespace
vtkm::filter::flow::internal::ParticleAdvector<DSIType> pav(
boundsMap, dsi, this->UseThreadedAlgorithm, this->UseAsynchronousCommunication);
VTKM_CONT vtkm::cont::PartitionedDataSet FilterParticleAdvectionSteadyState::DoExecutePartitions(
const vtkm::cont::PartitionedDataSet& input)
{
using DSIType = vtkm::filter::flow::internal::DataSetIntegratorSteadyState;
this->ValidateOptions();
using VariantType = vtkm::cont::Variant<std::string, std::pair<std::string, std::string>>;
VariantType variant;
if (this->VecFieldType == vtkm::filter::flow::VectorFieldType::VELOCITY_FIELD_TYPE)
{
const auto& field = this->GetActiveFieldName();
variant.Emplace<DSIType::VelocityFieldNameType>(field);
}
else if (this->VecFieldType == vtkm::filter::flow::VectorFieldType::ELECTRO_MAGNETIC_FIELD_TYPE)
{
const auto& electric = this->GetEField();
const auto& magnetic = this->GetBField();
variant.Emplace<DSIType::ElectroMagneticFieldNameType>(electric, magnetic);
}
vtkm::filter::flow::internal::BoundsMap boundsMap;
if (this->BlockIdsSet)
boundsMap = vtkm::filter::flow::internal::BoundsMap(input, this->BlockIds);
else
boundsMap = vtkm::filter::flow::internal::BoundsMap(input);
auto dsi = CreateDataSetIntegrators(
input, variant, boundsMap, this->SolverType, this->VecFieldType, this->GetResultType());
vtkm::filter::flow::internal::ParticleAdvector<DSIType> pav(boundsMap,
dsi,
this->UseThreadedAlgorithm,
this->UseAsynchronousCommunication,
this->GetResultType());
return pav.Execute(this->NumberOfSteps, this->StepSize, this->Seeds);
vtkm::cont::ArrayHandle<ParticleType> particles;
this->Seeds.AsArrayHandle(particles);
return pav.Execute(particles, this->StepSize);
}
}
}
} // namespace vtkm::filter::flow
#include <vtkm/filter/flow/ParticleAdvection.h>
#include <vtkm/filter/flow/Streamline.h>
#include <vtkm/filter/flow/WarpXStreamline.h>
namespace vtkm
{
namespace filter
{
namespace flow
{
template class FilterParticleAdvectionSteadyState<vtkm::filter::flow::ParticleAdvection>;
template class FilterParticleAdvectionSteadyState<vtkm::filter::flow::Streamline>;
template class FilterParticleAdvectionSteadyState<vtkm::filter::flow::WarpXStreamline>;
} // namespace flow
} // namespace filter
} // namespace vtkm

24
vtkm/filter/flow/FilterParticleAdvectionSteadyState.h
View File

@ -21,13 +21,29 @@ namespace filter
{
namespace flow
{
template <typename Derived>
struct FlowTraits;
template <typename Derived>
class VTKM_FILTER_FLOW_EXPORT FilterParticleAdvectionSteadyState : public FilterParticleAdvection
{
private:
VTKM_CONT vtkm::cont::PartitionedDataSet DoExecutePartitions(
const vtkm::cont::PartitionedDataSet& inData) override;
};
public:
using ParticleType = typename FlowTraits<Derived>::ParticleType;
using FieldType = typename FlowTraits<Derived>::FieldType;
using TerminationType = typename FlowTraits<Derived>::TerminationType;
using AnalysisType = typename FlowTraits<Derived>::AnalysisType;
private:
VTKM_CONT FieldType GetField(const vtkm::cont::DataSet& data) const;
VTKM_CONT TerminationType GetTermination(const vtkm::cont::DataSet& data) const;
VTKM_CONT AnalysisType GetAnalysis(const vtkm::cont::DataSet& data) const;
VTKM_CONT vtkm::cont::PartitionedDataSet DoExecutePartitions(
const vtkm::cont::PartitionedDataSet& input) override;
};
}
}
} // namespace vtkm::filter::flow

121
vtkm/filter/flow/FilterParticleAdvectionUnsteadyState.cxx
View File

@ -9,6 +9,8 @@
//============================================================================
#include <vtkm/filter/flow/FilterParticleAdvectionUnsteadyState.h>
#include <vtkm/filter/flow/internal/BoundsMap.h>
#include <vtkm/filter/flow/internal/DataSetIntegratorUnsteadyState.h>
#include <vtkm/filter/flow/internal/ParticleAdvector.h>
@ -18,72 +20,95 @@ namespace filter
{
namespace flow
{
namespace
template <typename Derived>
VTKM_CONT typename FilterParticleAdvectionUnsteadyState<Derived>::FieldType
FilterParticleAdvectionUnsteadyState<Derived>::GetField(const vtkm::cont::DataSet& data) const
{
VTKM_CONT std::vector<vtkm::filter::flow::internal::DataSetIntegratorUnsteadyState>
CreateDataSetIntegrators(const vtkm::cont::PartitionedDataSet& input,
const vtkm::cont::PartitionedDataSet& input2,
const std::string& activeField,
const vtkm::FloatDefault timer1,
const vtkm::FloatDefault timer2,
const vtkm::filter::flow::internal::BoundsMap& boundsMap,
const vtkm::filter::flow::IntegrationSolverType solverType,
const vtkm::filter::flow::VectorFieldType vecFieldType,
const vtkm::filter::flow::FlowResultType resultType)
const Derived* inst = static_cast<const Derived*>(this);
return inst->GetField(data);
}
template <typename Derived>
VTKM_CONT typename FilterParticleAdvectionUnsteadyState<Derived>::TerminationType
FilterParticleAdvectionUnsteadyState<Derived>::GetTermination(const vtkm::cont::DataSet& data) const
{
using DSIType = vtkm::filter::flow::internal::DataSetIntegratorUnsteadyState;
const Derived* inst = static_cast<const Derived*>(this);
return inst->GetTermination(data);
}
template <typename Derived>
VTKM_CONT typename FilterParticleAdvectionUnsteadyState<Derived>::AnalysisType
FilterParticleAdvectionUnsteadyState<Derived>::GetAnalysis(const vtkm::cont::DataSet& data) const
{
const Derived* inst = static_cast<const Derived*>(this);
return inst->GetAnalysis(data);
}
template <typename Derived>
VTKM_CONT vtkm::cont::PartitionedDataSet
FilterParticleAdvectionUnsteadyState<Derived>::DoExecutePartitions(
const vtkm::cont::PartitionedDataSet& input)
{
this->ValidateOptions();
using DSIType = vtkm::filter::flow::internal::
DataSetIntegratorUnsteadyState<ParticleType, FieldType, TerminationType, AnalysisType>;
vtkm::filter::flow::internal::BoundsMap boundsMap(input);
std::vector<DSIType> dsi;
for (vtkm::Id i = 0; i < input.GetNumberOfPartitions(); i++)
{
vtkm::Id blockId = boundsMap.GetLocalBlockId(i);
auto ds1 = input.GetPartition(i);
auto ds2 = input2.GetPartition(i);
if ((!ds1.HasPointField(activeField) && !ds1.HasCellField(activeField)) ||
(!ds2.HasPointField(activeField) && !ds2.HasCellField(activeField)))
throw vtkm::cont::ErrorFilterExecution("Unsupported field assocation");
auto ds2 = this->Input2.GetPartition(i);
dsi.emplace_back(
ds1, ds2, timer1, timer2, blockId, activeField, solverType, vecFieldType, resultType);
}
// Build the field for the current dataset
FieldType field1 = this->GetField(ds1);
FieldType field2 = this->GetField(ds2);
return dsi;
}
} // anonymous namespace
// Build the termination for the current dataset
TerminationType termination = this->GetTermination(ds1);
VTKM_CONT void FilterParticleAdvectionUnsteadyState::ValidateOptions() const
{
this->FilterParticleAdvection::ValidateOptions();
if (this->Time1 >= this->Time2)
throw vtkm::cont::ErrorFilterExecution("PreviousTime must be less than NextTime");
}
AnalysisType analysis = this->GetAnalysis(ds1);
VTKM_CONT vtkm::cont::PartitionedDataSet FilterParticleAdvectionUnsteadyState::DoExecutePartitions(
const vtkm::cont::PartitionedDataSet& input)
{
using DSIType = vtkm::filter::flow::internal::DataSetIntegratorUnsteadyState;
this->ValidateOptions();
vtkm::filter::flow::internal::BoundsMap boundsMap(input);
auto dsi = CreateDataSetIntegrators(input,
this->Input2,
this->GetActiveFieldName(),
dsi.emplace_back(blockId,
field1,
field2,
ds1,
ds2,
this->Time1,
this->Time2,
boundsMap,
this->SolverType,
this->VecFieldType,
this->GetResultType());
termination,
analysis);
}
vtkm::filter::flow::internal::ParticleAdvector<DSIType> pav(
boundsMap, dsi, this->UseThreadedAlgorithm, this->UseAsynchronousCommunication);
vtkm::filter::flow::internal::ParticleAdvector<DSIType> pav(boundsMap,
dsi,
this->UseThreadedAlgorithm,
this->UseAsynchronousCommunication,
this->GetResultType());
return pav.Execute(this->NumberOfSteps, this->StepSize, this->Seeds);
vtkm::cont::ArrayHandle<ParticleType> particles;
this->Seeds.AsArrayHandle(particles);
return pav.Execute(particles, this->StepSize);
}
}
}
} // namespace vtkm::filter::flow
#include <vtkm/filter/flow/PathParticle.h>
#include <vtkm/filter/flow/Pathline.h>
namespace vtkm
{
namespace filter
{
namespace flow
{
template class FilterParticleAdvectionUnsteadyState<vtkm::filter::flow::PathParticle>;
template class FilterParticleAdvectionUnsteadyState<vtkm::filter::flow::Pathline>;
} // namespace flow
} // namespace filter
} // namespace vtkm

20
vtkm/filter/flow/FilterParticleAdvectionUnsteadyState.h
View File

@ -21,9 +21,18 @@ namespace filter
namespace flow
{
template <typename Derived>
struct FlowTraits;
template <typename Derived>
class VTKM_FILTER_FLOW_EXPORT FilterParticleAdvectionUnsteadyState : public FilterParticleAdvection
{
public:
using ParticleType = typename FlowTraits<Derived>::ParticleType;
using FieldType = typename FlowTraits<Derived>::FieldType;
using TerminationType = typename FlowTraits<Derived>::TerminationType;
using AnalysisType = typename FlowTraits<Derived>::AnalysisType;
VTKM_CONT void SetPreviousTime(vtkm::FloatDefault t1) { this->Time1 = t1; }
VTKM_CONT void SetNextTime(vtkm::FloatDefault t2) { this->Time2 = t2; }
@ -32,12 +41,15 @@ public:
VTKM_CONT void SetNextDataSet(const vtkm::cont::PartitionedDataSet& pds) { this->Input2 = pds; }
protected:
VTKM_CONT virtual void ValidateOptions() const override;
private:
VTKM_CONT FieldType GetField(const vtkm::cont::DataSet& data) const;
VTKM_CONT TerminationType GetTermination(const vtkm::cont::DataSet& data) const;
VTKM_CONT AnalysisType GetAnalysis(const vtkm::cont::DataSet& data) const;
VTKM_CONT vtkm::cont::PartitionedDataSet DoExecutePartitions(
const vtkm::cont::PartitionedDataSet& inData) override;
const vtkm::cont::PartitionedDataSet& input);
vtkm::cont::PartitionedDataSet Input2;
vtkm::FloatDefault Time1 = -1;

16
vtkm/filter/flow/Lagrangian.cxx
View File

@ -21,6 +21,14 @@
#include <vtkm/filter/flow/worklet/RK4Integrator.h>
#include <vtkm/filter/flow/worklet/Stepper.h>
#include <vtkm/filter/flow/worklet/Analysis.h>
#include <vtkm/filter/flow/worklet/Termination.h>
//#include <vtkm/worklet/WorkletMapField.h>
//#include <cstring>
//#include <sstream>
//#include <string.h>
namespace vtkm
{
namespace filter
@ -220,7 +228,6 @@ VTKM_CONT vtkm::cont::DataSet Lagrangian::DoExecute(const vtkm::cont::DataSet& i
using Stepper = vtkm::worklet::flow::Stepper<RK4Type, GridEvalType>;
vtkm::worklet::flow::ParticleAdvection particleadvection;
vtkm::worklet::flow::ParticleAdvectionResult<vtkm::Particle> res;
const auto field = input.GetField(this->GetActiveFieldName());
FieldType velocities(field.GetData().AsArrayHandle<vtkm::cont::ArrayHandle<vtkm::Vec3f>>(),
@ -228,9 +235,10 @@ VTKM_CONT vtkm::cont::DataSet Lagrangian::DoExecute(const vtkm::cont::DataSet& i
GridEvalType gridEval(coords, cells, velocities);
Stepper rk4(gridEval, static_cast<vtkm::Float32>(this->StepSize));
res = particleadvection.Run(rk4, basisParticleArray, 1); // Taking a single step
auto particles = res.Particles;
vtkm::worklet::flow::NormalTermination termination(1);
vtkm::worklet::flow::NoAnalysis<vtkm::Particle> analysis;
particleadvection.Run(rk4, basisParticleArray, termination, analysis); // Taking a single step
auto particles = analysis.Particles;
vtkm::cont::DataSet outputData;

14
vtkm/filter/flow/LagrangianStructures.cxx
View File

@ -13,14 +13,17 @@
#include <vtkm/cont/ArrayCopy.h>
#include <vtkm/cont/DataSetBuilderUniform.h>
#include <vtkm/cont/ErrorFilterExecution.h>
#include <vtkm/cont/Invoker.h>
#include <vtkm/filter/flow/LagrangianStructures.h>
#include <vtkm/filter/flow/worklet/Analysis.h>
#include <vtkm/filter/flow/worklet/Field.h>
#include <vtkm/filter/flow/worklet/GridEvaluators.h>
#include <vtkm/filter/flow/worklet/LagrangianStructures.h>
#include <vtkm/filter/flow/worklet/ParticleAdvection.h>
#include <vtkm/filter/flow/worklet/RK4Integrator.h>
#include <vtkm/filter/flow/worklet/Stepper.h>
#include <vtkm/filter/flow/worklet/Termination.h>
namespace
{
@ -144,11 +147,14 @@ VTKM_CONT vtkm::cont::DataSet LagrangianStructures::DoExecute(const vtkm::cont::
GridEvaluator evaluator(input.GetCoordinateSystem(), input.GetCellSet(), velocities);
Stepper integrator(evaluator, stepSize);
vtkm::worklet::flow::ParticleAdvection particles;
vtkm::worklet::flow::ParticleAdvectionResult<vtkm::Particle> advectionResult;
vtkm::worklet::flow::NormalTermination termination(numberOfSteps);
vtkm::worklet::flow::NoAnalysis<vtkm::Particle> analysis;
vtkm::cont::ArrayHandle<vtkm::Particle> advectionPoints;
this->Invoke(detail::MakeParticles{}, lcsInputPoints, advectionPoints);
advectionResult = particles.Run(integrator, advectionPoints, numberOfSteps);
this->Invoke(detail::ExtractParticlePosition{}, advectionResult.Particles, lcsOutputPoints);
vtkm::cont::Invoker invoke;
invoke(detail::MakeParticles{}, lcsInputPoints, advectionPoints);
particles.Run(integrator, advectionPoints, termination, analysis);
invoke(detail::ExtractParticlePosition{}, analysis.Particles, lcsOutputPoints);
}
// FTLE output field
vtkm::cont::ArrayHandle<vtkm::FloatDefault> outputField;

33
vtkm/filter/flow/ParticleAdvection.cxx
View File

@ -17,9 +17,38 @@ namespace filter
namespace flow
{
VTKM_CONT vtkm::filter::flow::FlowResultType ParticleAdvection::GetResultType() const
// using ParticleType = vtkm::Particle;
// using TerminationType = vtkm::worklet::flow::NormalTermination;
// using AnalysisType = vtkm::worklet::flow::Streamline;
// using ArrayType = vtkm::cont::ArrayHandle<vtkm::Vec3f>;
// using FieldType = vtkm::worklet::flow::VelocityField<ArrayType>;
VTKM_CONT ParticleAdvection::FieldType ParticleAdvection::GetField(
const vtkm::cont::DataSet& dataset) const
{
return vtkm::filter::flow::FlowResultType::PARTICLE_ADVECT_TYPE;
const auto& fieldNm = this->GetActiveFieldName();
if (!dataset.HasPointField(fieldNm) && !dataset.HasCellField(fieldNm))
throw vtkm::cont::ErrorFilterExecution("Unsupported field assocation");
auto assoc = dataset.GetField(fieldNm).GetAssociation();
ArrayType arr;
vtkm::cont::ArrayCopyShallowIfPossible(dataset.GetField(fieldNm).GetData(), arr);
return ParticleAdvection::FieldType(arr, assoc);
}
VTKM_CONT ParticleAdvection::TerminationType ParticleAdvection::GetTermination(
const vtkm::cont::DataSet& dataset) const
{
// dataset not used
(void)dataset;
return ParticleAdvection::TerminationType(this->NumberOfSteps);
}
VTKM_CONT ParticleAdvection::AnalysisType ParticleAdvection::GetAnalysis(
const vtkm::cont::DataSet& dataset) const
{
// dataset not used
(void)dataset;
return ParticleAdvection::AnalysisType{};
}
}

31
vtkm/filter/flow/ParticleAdvection.h
View File

@ -15,6 +15,10 @@
#include <vtkm/filter/flow/FlowTypes.h>
#include <vtkm/filter/flow/vtkm_filter_flow_export.h>
#include <vtkm/filter/flow/worklet/Analysis.h>
#include <vtkm/filter/flow/worklet/Field.h>
#include <vtkm/filter/flow/worklet/Termination.h>
namespace vtkm
{
namespace filter
@ -22,17 +26,38 @@ namespace filter
namespace flow
{
class ParticleAdvection;
template <>
struct FlowTraits<ParticleAdvection>
{
using ParticleType = vtkm::Particle;
using TerminationType = vtkm::worklet::flow::NormalTermination;
using AnalysisType = vtkm::worklet::flow::NoAnalysis<ParticleType>;
using ArrayType = vtkm::cont::ArrayHandle<vtkm::Vec3f>;
using FieldType = vtkm::worklet::flow::VelocityField<ArrayType>;
};
/// \brief Advect particles in a vector field.
/// Takes as input a vector field and seed locations and generates the
/// end points for each seed through the vector field.
class VTKM_FILTER_FLOW_EXPORT ParticleAdvection
: public vtkm::filter::flow::FilterParticleAdvectionSteadyState
: public vtkm::filter::flow::FilterParticleAdvectionSteadyState<ParticleAdvection>
{
public:
private:
VTKM_CONT vtkm::filter::flow::FlowResultType GetResultType() const override;
using ParticleType = typename FlowTraits<ParticleAdvection>::ParticleType;
using TerminationType = typename FlowTraits<ParticleAdvection>::TerminationType;
using AnalysisType = typename FlowTraits<ParticleAdvection>::AnalysisType;
using ArrayType = typename FlowTraits<ParticleAdvection>::ArrayType;
using FieldType = typename FlowTraits<ParticleAdvection>::FieldType;
VTKM_CONT FieldType GetField(const vtkm::cont::DataSet& data) const;
VTKM_CONT TerminationType GetTermination(const vtkm::cont::DataSet& data) const;
VTKM_CONT AnalysisType GetAnalysis(const vtkm::cont::DataSet& data) const;
};
}

31
vtkm/filter/flow/PathParticle.cxx
View File

@ -17,11 +17,38 @@ namespace filter
namespace flow
{
VTKM_CONT vtkm::filter::flow::FlowResultType PathParticle::GetResultType() const
VTKM_CONT PathParticle::FieldType PathParticle::GetField(const vtkm::cont::DataSet& dataset) const
{
return vtkm::filter::flow::FlowResultType::PARTICLE_ADVECT_TYPE;
const auto& fieldNm = this->GetActiveFieldName();
if (!dataset.HasPointField(fieldNm) && !dataset.HasCellField(fieldNm))
throw vtkm::cont::ErrorFilterExecution("Unsupported field assocation");
auto assoc = dataset.GetField(fieldNm).GetAssociation();
ArrayType arr;
vtkm::cont::ArrayCopyShallowIfPossible(dataset.GetField(fieldNm).GetData(), arr);
return PathParticle::FieldType(arr, assoc);
}
VTKM_CONT PathParticle::TerminationType PathParticle::GetTermination(
const vtkm::cont::DataSet& dataset) const
{
// dataset not used
(void)dataset;
return PathParticle::TerminationType(this->NumberOfSteps);
}
VTKM_CONT PathParticle::AnalysisType PathParticle::GetAnalysis(
const vtkm::cont::DataSet& dataset) const
{
// dataset not used
(void)dataset;
return PathParticle::AnalysisType();
}
//VTKM_CONT vtkm::filter::flow::FlowResultType PathParticle::GetResultType() const
//{
// return vtkm::filter::flow::FlowResultType::PARTICLE_ADVECT_TYPE;
//}
}
}
} // namespace vtkm::filter::flow

32
vtkm/filter/flow/PathParticle.h
View File

@ -15,6 +15,10 @@
#include <vtkm/filter/flow/FlowTypes.h>
#include <vtkm/filter/flow/vtkm_filter_flow_export.h>
#include <vtkm/filter/flow/worklet/Analysis.h>
#include <vtkm/filter/flow/worklet/Field.h>
#include <vtkm/filter/flow/worklet/Termination.h>
namespace vtkm
{
namespace filter
@ -22,16 +26,38 @@ namespace filter
namespace flow
{
class PathParticle;
template <>
struct FlowTraits<PathParticle>
{
using ParticleType = vtkm::Particle;
using TerminationType = vtkm::worklet::flow::NormalTermination;
using AnalysisType = vtkm::worklet::flow::NoAnalysis<ParticleType>;
using ArrayType = vtkm::cont::ArrayHandle<vtkm::Vec3f>;
using FieldType = vtkm::worklet::flow::VelocityField<ArrayType>;
};
/// \brief Advect particles in a vector field.
/// Takes as input a vector field and seed locations and generates the
/// end points for each seed through the vector field.
class VTKM_FILTER_FLOW_EXPORT PathParticle
: public vtkm::filter::flow::FilterParticleAdvectionUnsteadyState
: public vtkm::filter::flow::FilterParticleAdvectionUnsteadyState<PathParticle>
{
private:
VTKM_CONT vtkm::filter::flow::FlowResultType GetResultType() const override;
public:
using ParticleType = typename FlowTraits<PathParticle>::ParticleType;
using TerminationType = typename FlowTraits<PathParticle>::TerminationType;
using AnalysisType = typename FlowTraits<PathParticle>::AnalysisType;
using ArrayType = typename FlowTraits<PathParticle>::ArrayType;
using FieldType = typename FlowTraits<PathParticle>::FieldType;
VTKM_CONT FieldType GetField(const vtkm::cont::DataSet& data) const;
VTKM_CONT TerminationType GetTermination(const vtkm::cont::DataSet& data) const;
VTKM_CONT AnalysisType GetAnalysis(const vtkm::cont::DataSet& data) const;
};
}

29
vtkm/filter/flow/Pathline.cxx
View File

@ -17,11 +17,36 @@ namespace filter
namespace flow
{
VTKM_CONT vtkm::filter::flow::FlowResultType Pathline::GetResultType() const
VTKM_CONT Pathline::FieldType Pathline::GetField(const vtkm::cont::DataSet& dataset) const
{
return vtkm::filter::flow::FlowResultType::STREAMLINE_TYPE;
const auto& fieldNm = this->GetActiveFieldName();
if (!dataset.HasPointField(fieldNm) && !dataset.HasCellField(fieldNm))
throw vtkm::cont::ErrorFilterExecution("Unsupported field assocation");
auto assoc = dataset.GetField(fieldNm).GetAssociation();
ArrayType arr;
vtkm::cont::ArrayCopyShallowIfPossible(dataset.GetField(fieldNm).GetData(), arr);
return vtkm::worklet::flow::VelocityField<ArrayType>(arr, assoc);
}
VTKM_CONT Pathline::TerminationType Pathline::GetTermination(
const vtkm::cont::DataSet& dataset) const
{
// dataset not used
(void)dataset;
return Pathline::TerminationType(this->NumberOfSteps);
}
VTKM_CONT Pathline::AnalysisType Pathline::GetAnalysis(const vtkm::cont::DataSet& dataset) const
{
// dataset not used
(void)dataset;
return Pathline::AnalysisType(this->NumberOfSteps);
}
//VTKM_CONT vtkm::filter::flow::FlowResultType Pathline::GetResultType() const
//{
// return vtkm::filter::flow::FlowResultType::STREAMLINE_TYPE;
//}
}
}
} // namespace vtkm::filter::flow

32
vtkm/filter/flow/Pathline.h
View File

@ -15,6 +15,10 @@
#include <vtkm/filter/flow/FlowTypes.h>
#include <vtkm/filter/flow/vtkm_filter_flow_export.h>
#include <vtkm/filter/flow/worklet/Analysis.h>
#include <vtkm/filter/flow/worklet/Field.h>
#include <vtkm/filter/flow/worklet/Termination.h>
namespace vtkm
{
namespace filter
@ -22,16 +26,38 @@ namespace filter
namespace flow
{
class Pathline;
template <>
struct FlowTraits<Pathline>
{
using ParticleType = vtkm::Particle;
using TerminationType = vtkm::worklet::flow::NormalTermination;
using AnalysisType = vtkm::worklet::flow::StreamlineAnalysis<ParticleType>;
using ArrayType = vtkm::cont::ArrayHandle<vtkm::Vec3f>;
using FieldType = vtkm::worklet::flow::VelocityField<ArrayType>;
};
/// \brief Advect particles in a vector field.
/// Takes as input a vector field and seed locations and generates the
/// end points for each seed through the vector field.
class VTKM_FILTER_FLOW_EXPORT Pathline
: public vtkm::filter::flow::FilterParticleAdvectionUnsteadyState
: public vtkm::filter::flow::FilterParticleAdvectionUnsteadyState<Pathline>
{
private:
VTKM_CONT vtkm::filter::flow::FlowResultType GetResultType() const override;
public:
using ParticleType = typename FlowTraits<Pathline>::ParticleType;
using TerminationType = typename FlowTraits<Pathline>::TerminationType;
using AnalysisType = typename FlowTraits<Pathline>::AnalysisType;
using ArrayType = typename FlowTraits<Pathline>::ArrayType;
using FieldType = typename FlowTraits<Pathline>::FieldType;
VTKM_CONT FieldType GetField(const vtkm::cont::DataSet& data) const;
VTKM_CONT TerminationType GetTermination(const vtkm::cont::DataSet& data) const;
VTKM_CONT AnalysisType GetAnalysis(const vtkm::cont::DataSet& data) const;
};
}

19
vtkm/filter/flow/StreamSurface.cxx
View File

@ -1,4 +1,5 @@
//============================================================================
//=============================================================================
//
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.txt for details.
@ -6,7 +7,8 @@
// 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.
//============================================================================
//
//=============================================================================
#include <vtkm/cont/ArrayCopy.h>
#include <vtkm/cont/ErrorFilterExecution.h>
@ -65,19 +67,22 @@ VTKM_CONT vtkm::cont::DataSet StreamSurface::DoExecute(const vtkm::cont::DataSet
GridEvalType eval(coords, cells, velocities);
Stepper rk4(eval, this->StepSize);
vtkm::worklet::flow::Streamline streamline;
using ParticleArray = vtkm::cont::ArrayHandle<vtkm::Particle>;
vtkm::cont::ArrayHandle<vtkm::Particle> seedArray;
vtkm::cont::ArrayCopy(this->Seeds.AsArrayHandle<ParticleArray>(), seedArray);
auto res = streamline.Run(rk4, seedArray, this->NumberOfSteps);
vtkm::worklet::flow::ParticleAdvection worklet;
vtkm::worklet::flow::NormalTermination termination(this->NumberOfSteps);
vtkm::worklet::flow::StreamlineAnalysis<vtkm::Particle> analysis(this->NumberOfSteps);
worklet.Run(rk4, seedArray, termination, analysis);
//compute surface from streamlines
vtkm::worklet::flow::StreamSurface streamSurface;
vtkm::cont::ArrayHandle<vtkm::Vec3f> srfPoints;
vtkm::cont::CellSetSingleType<> srfCells;
vtkm::cont::CoordinateSystem slCoords("coordinates", res.Positions);
streamSurface.Run(slCoords, res.PolyLines, srfPoints, srfCells);
vtkm::cont::CoordinateSystem slCoords("coordinates", analysis.Streams);
streamSurface.Run(slCoords, analysis.PolyLines, srfPoints, srfCells);
vtkm::cont::DataSet outData;
vtkm::cont::CoordinateSystem outputCoords("coordinates", srfPoints);

31
vtkm/filter/flow/Streamline.cxx
View File

@ -17,9 +17,36 @@ namespace filter
namespace flow
{
VTKM_CONT vtkm::filter::flow::FlowResultType Streamline::GetResultType() const
// using ParticleType = vtkm::Particle;
// using TerminationType = vtkm::worklet::flow::NormalTermination;
// using AnalysisType = vtkm::worklet::flow::Streamline;
// using ArrayType = vtkm::cont::ArrayHandle<vtkm::Vec3f>;
// using FieldType = vtkm::worklet::flow::VelocityField<ArrayType>;
VTKM_CONT Streamline::FieldType Streamline::GetField(const vtkm::cont::DataSet& dataset) const
{
return vtkm::filter::flow::FlowResultType::STREAMLINE_TYPE;
const auto& fieldNm = this->GetActiveFieldName();
if (!dataset.HasPointField(fieldNm) && !dataset.HasCellField(fieldNm))
throw vtkm::cont::ErrorFilterExecution("Unsupported field assocation");
auto assoc = dataset.GetField(fieldNm).GetAssociation();
ArrayType arr;
vtkm::cont::ArrayCopyShallowIfPossible(dataset.GetField(fieldNm).GetData(), arr);
return vtkm::worklet::flow::VelocityField<ArrayType>(arr, assoc);
}
VTKM_CONT Streamline::TerminationType Streamline::GetTermination(
const vtkm::cont::DataSet& dataset) const
{
// dataset not used
(void)dataset;
return Streamline::TerminationType(this->NumberOfSteps);
}
VTKM_CONT Streamline::AnalysisType Streamline::GetAnalysis(const vtkm::cont::DataSet& dataset) const
{
// dataset not used
(void)dataset;
return Streamline::AnalysisType(this->NumberOfSteps);
}
}

32
vtkm/filter/flow/Streamline.h
View File

@ -15,6 +15,10 @@
#include <vtkm/filter/flow/FlowTypes.h>
#include <vtkm/filter/flow/vtkm_filter_flow_export.h>
#include <vtkm/filter/flow/worklet/Analysis.h>
#include <vtkm/filter/flow/worklet/Field.h>
#include <vtkm/filter/flow/worklet/Termination.h>
namespace vtkm
{
namespace filter
@ -22,16 +26,38 @@ namespace filter
namespace flow
{
class Streamline;
template <>
struct FlowTraits<Streamline>
{
using ParticleType = vtkm::Particle;
using TerminationType = vtkm::worklet::flow::NormalTermination;
using AnalysisType = vtkm::worklet::flow::StreamlineAnalysis<ParticleType>;
using ArrayType = vtkm::cont::ArrayHandle<vtkm::Vec3f>;
using FieldType = vtkm::worklet::flow::VelocityField<ArrayType>;
};
/// \brief Advect particles in a vector field.
/// Takes as input a vector field and seed locations and generates the
/// end points for each seed through the vector field.
class VTKM_FILTER_FLOW_EXPORT Streamline
: public vtkm::filter::flow::FilterParticleAdvectionSteadyState
: public vtkm::filter::flow::FilterParticleAdvectionSteadyState<Streamline>
{
private:
VTKM_CONT vtkm::filter::flow::FlowResultType GetResultType() const override;
public:
using ParticleType = typename FlowTraits<Streamline>::ParticleType;
using TerminationType = typename FlowTraits<Streamline>::TerminationType;
using AnalysisType = typename FlowTraits<Streamline>::AnalysisType;
using ArrayType = typename FlowTraits<Streamline>::ArrayType;
using FieldType = typename FlowTraits<Streamline>::FieldType;
VTKM_CONT FieldType GetField(const vtkm::cont::DataSet& data) const;
VTKM_CONT TerminationType GetTermination(const vtkm::cont::DataSet& data) const;
VTKM_CONT AnalysisType GetAnalysis(const vtkm::cont::DataSet& data) const;
};
}

59
vtkm/filter/flow/WarpXStreamline.cxx Normal file
View File

@ -0,0 +1,59 @@
//============================================================================
// 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.
//============================================================================
#include <vtkm/filter/flow/WarpXStreamline.h>
namespace vtkm
{
namespace filter
{
namespace flow
{
VTKM_CONT WarpXStreamline::FieldType WarpXStreamline::GetField(
const vtkm::cont::DataSet& dataset) const
{
const auto& electric = this->GetEField();
const auto& magnetic = this->GetBField();
if (!dataset.HasPointField(electric) && !dataset.HasCellField(electric))
throw vtkm::cont::ErrorFilterExecution("Unsupported field assocation");
if (!dataset.HasPointField(magnetic) && !dataset.HasCellField(magnetic))
throw vtkm::cont::ErrorFilterExecution("Unsupported field assocation");
auto eAssoc = dataset.GetField(electric).GetAssociation();
auto bAssoc = dataset.GetField(magnetic).GetAssociation();
if (eAssoc != bAssoc)
{
throw vtkm::cont::ErrorFilterExecution("E and B field need to have same association");
}
ArrayType eField, bField;
vtkm::cont::ArrayCopyShallowIfPossible(dataset.GetField(electric).GetData(), eField);
vtkm::cont::ArrayCopyShallowIfPossible(dataset.GetField(magnetic).GetData(), bField);
return vtkm::worklet::flow::ElectroMagneticField<ArrayType>(eField, bField, eAssoc);
}
VTKM_CONT WarpXStreamline::TerminationType WarpXStreamline::GetTermination(
const vtkm::cont::DataSet& dataset) const
{
// dataset not used
(void)dataset;
return WarpXStreamline::TerminationType(this->NumberOfSteps);
}
VTKM_CONT WarpXStreamline::AnalysisType WarpXStreamline::GetAnalysis(
const vtkm::cont::DataSet& dataset) const
{
// dataset not used
(void)dataset;
return WarpXStreamline::AnalysisType(this->NumberOfSteps);
}
}
}
} // namespace vtkm::filter::flow

76
vtkm/filter/flow/WarpXStreamline.h Normal file
View File

@ -0,0 +1,76 @@
//============================================================================
// 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_filter_flow_WarpXStreamline_h
#define vtk_m_filter_flow_WarpXStreamline_h
#include <vtkm/filter/flow/FilterParticleAdvectionSteadyState.h>
#include <vtkm/filter/flow/FlowTypes.h>
#include <vtkm/filter/flow/vtkm_filter_flow_export.h>
#include <vtkm/filter/flow/worklet/Analysis.h>
#include <vtkm/filter/flow/worklet/Field.h>
#include <vtkm/filter/flow/worklet/Termination.h>
namespace vtkm
{
namespace filter
{
namespace flow
{
class WarpXStreamline;
template <>
struct FlowTraits<WarpXStreamline>
{
using ParticleType = vtkm::ChargedParticle;
using TerminationType = vtkm::worklet::flow::NormalTermination;
using AnalysisType = vtkm::worklet::flow::StreamlineAnalysis<ParticleType>;
using ArrayType = vtkm::cont::ArrayHandle<vtkm::Vec3f>;
using FieldType = vtkm::worklet::flow::ElectroMagneticField<ArrayType>;
};
/// \brief Advect particles in a vector field.
/// Takes as input a vector field and seed locations and generates the
/// end points for each seed through the vector field.
class VTKM_FILTER_FLOW_EXPORT WarpXStreamline
: public vtkm::filter::flow::FilterParticleAdvectionSteadyState<WarpXStreamline>
{
public:
using ParticleType = typename FlowTraits<WarpXStreamline>::ParticleType;
using TerminationType = typename FlowTraits<WarpXStreamline>::TerminationType;
using AnalysisType = typename FlowTraits<WarpXStreamline>::AnalysisType;
using ArrayType = typename FlowTraits<WarpXStreamline>::ArrayType;
using FieldType = typename FlowTraits<WarpXStreamline>::FieldType;
VTKM_CONT WarpXStreamline() { this->SetSolverEuler(); }
VTKM_CONT FieldType GetField(const vtkm::cont::DataSet& data) const;
VTKM_CONT TerminationType GetTermination(const vtkm::cont::DataSet& data) const;
VTKM_CONT AnalysisType GetAnalysis(const vtkm::cont::DataSet& data) const;
VTKM_CONT void SetEField(const std::string& name) { this->SetActiveField(0, name); }
VTKM_CONT void SetBField(const std::string& name) { this->SetActiveField(1, name); }
VTKM_CONT std::string GetEField() const { return this->GetActiveFieldName(0); }
VTKM_CONT std::string GetBField() const { return this->GetActiveFieldName(1); }
};
}
}
} // namespace vtkm::filter::flow
#endif // vtk_m_filter_flow_WarpXStreamline_h

24
vtkm/filter/flow/internal/AdvectAlgorithm.h
View File

@ -25,10 +25,12 @@ namespace flow
namespace internal
{
template <typename DSIType, template <typename> class ResultType, typename ParticleType>
template <typename DSIType>
class AdvectAlgorithm
{
public:
using ParticleType = typename DSIType::PType;
AdvectAlgorithm(const vtkm::filter::flow::internal::BoundsMap& bm,
std::vector<DSIType>& blocks,
bool useAsyncComm)
@ -40,11 +42,8 @@ public:
{
}
void Execute(vtkm::Id maxNumSteps,
vtkm::FloatDefault stepSize,
const vtkm::cont::ArrayHandle<ParticleType>& seeds)
void Execute(const vtkm::cont::ArrayHandle<ParticleType>& seeds, vtkm::FloatDefault stepSize)
{
this->SetMaxNumberOfSteps(maxNumSteps);
this->SetStepSize(stepSize);
this->SetSeeds(seeds);
this->Go();
@ -53,19 +52,17 @@ public:
vtkm::cont::PartitionedDataSet GetOutput() const
{
vtkm::cont::PartitionedDataSet output;
for (const auto& b : this->Blocks)
{
vtkm::cont::DataSet ds;
if (b.template GetOutput<ParticleType>(ds))
if (b.GetOutput(ds))
output.AppendPartition(ds);
}
return output;
}
void SetStepSize(vtkm::FloatDefault stepSize) { this->StepSize = stepSize; }
void SetMaxNumberOfSteps(vtkm::Id numSteps) { this->MaxNumberOfSteps = numSteps; }
void SetSeeds(const vtkm::cont::ArrayHandle<ParticleType>& seeds)
{
this->ClearParticles();
@ -110,10 +107,9 @@ public:
{
//make this a pointer to avoid the copy?
auto& block = this->GetDataSet(blockId);
DSIHelperInfoType bb =
DSIHelperInfo<ParticleType>(v, this->BoundsMap, this->ParticleBlockIDsMap);
block.Advect(bb, this->StepSize, this->MaxNumberOfSteps);
numTerm = this->UpdateResult(bb.Get<DSIHelperInfo<ParticleType>>());
DSIHelperInfo<ParticleType> bb(v, this->BoundsMap, this->ParticleBlockIDsMap);
block.Advect(bb, this->StepSize);
numTerm = this->UpdateResult(bb);
}
vtkm::Id numTermMessages = 0;
@ -145,7 +141,7 @@ public:
#endif
}
DataSetIntegrator<DSIType>& GetDataSet(vtkm::Id id)
DataSetIntegrator<DSIType, ParticleType>& GetDataSet(vtkm::Id id)
{
for (auto& it : this->Blocks)
if (it.GetID() == id)

35
vtkm/filter/flow/internal/AdvectAlgorithmThreaded.h
View File

@ -28,14 +28,16 @@ namespace flow
namespace internal
{
template <typename DSIType, template <typename> class ResultType, typename ParticleType>
class AdvectAlgorithmThreaded : public AdvectAlgorithm<DSIType, ResultType, ParticleType>
template <typename DSIType>
class AdvectAlgorithmThreaded : public AdvectAlgorithm<DSIType>
{
public:
using ParticleType = typename DSIType::PType;
AdvectAlgorithmThreaded(const vtkm::filter::flow::internal::BoundsMap& bm,
std::vector<DSIType>& blocks,
bool useAsyncComm)
: AdvectAlgorithm<DSIType, ResultType, ParticleType>(bm, blocks, useAsyncComm)
: AdvectAlgorithm<DSIType>(bm, blocks, useAsyncComm)
, Done(false)
, WorkerActivate(false)
{
@ -64,8 +66,7 @@ protected:
bool GetActiveParticles(std::vector<ParticleType>& particles, vtkm::Id& blockId) override
{
std::lock_guard<std::mutex> lock(this->Mutex);
bool val = this->AdvectAlgorithm<DSIType, ResultType, ParticleType>::GetActiveParticles(
particles, blockId);
bool val = this->AdvectAlgorithm<DSIType>::GetActiveParticles(particles, blockId);
this->WorkerActivate = val;
return val;
}
@ -76,7 +77,7 @@ protected:
if (!particles.empty())
{
std::lock_guard<std::mutex> lock(this->Mutex);
this->AdvectAlgorithm<DSIType, ResultType, ParticleType>::UpdateActive(particles, idsMap);
this->AdvectAlgorithm<DSIType>::UpdateActive(particles, idsMap);
//Let workers know there is new work
this->WorkerActivateCondition.notify_all();
@ -105,7 +106,7 @@ protected:
this->WorkerActivateCondition.wait(lock, [this] { return WorkerActivate || Done; });
}
void UpdateWorkerResult(vtkm::Id blockId, DSIHelperInfoType& b)
void UpdateWorkerResult(vtkm::Id blockId, DSIHelperInfo<ParticleType>& b)
{
std::lock_guard<std::mutex> lock(this->Mutex);
auto& it = this->WorkerResults[blockId];
@ -121,9 +122,8 @@ protected:
if (this->GetActiveParticles(v, blockId))
{
auto& block = this->GetDataSet(blockId);
DSIHelperInfoType bb =
DSIHelperInfo<ParticleType>(v, this->BoundsMap, this->ParticleBlockIDsMap);
block.Advect(bb, this->StepSize, this->MaxNumberOfSteps);
DSIHelperInfo<ParticleType> bb(v, this->BoundsMap, this->ParticleBlockIDsMap);
block.Advect(bb, this->StepSize);
this->UpdateWorkerResult(blockId, bb);
}
else
@ -146,17 +146,14 @@ protected:
while (this->TotalNumTerminatedParticles < this->TotalNumParticles)
{
std::unordered_map<vtkm::Id, std::vector<DSIHelperInfoType>> workerResults;
std::unordered_map<vtkm::Id, std::vector<DSIHelperInfo<ParticleType>>> workerResults;
this->GetWorkerResults(workerResults);
vtkm::Id numTerm = 0;
for (auto& it : workerResults)
{
for (auto& r : it.second)
{
numTerm += this->UpdateResult(r.Get<DSIHelperInfo<ParticleType>>());
}
numTerm += this->UpdateResult(r);
}
vtkm::Id numTermMessages = 0;
@ -177,12 +174,12 @@ protected:
if (this->Done)
return true;
return (this->AdvectAlgorithm<DSIType, ResultType, ParticleType>::GetBlockAndWait(
syncComm, numLocalTerm) &&
return (this->AdvectAlgorithm<DSIType>::GetBlockAndWait(syncComm, numLocalTerm) &&
!this->WorkerActivate && this->WorkerResults.empty());
}
void GetWorkerResults(std::unordered_map<vtkm::Id, std::vector<DSIHelperInfoType>>& results)
void GetWorkerResults(
std::unordered_map<vtkm::Id, std::vector<DSIHelperInfo<ParticleType>>>& results)
{
results.clear();
@ -198,7 +195,7 @@ protected:
std::mutex Mutex;
bool WorkerActivate;
std::condition_variable WorkerActivateCondition;
std::unordered_map<vtkm::Id, std::vector<DSIHelperInfoType>> WorkerResults;
std::unordered_map<vtkm::Id, std::vector<DSIHelperInfo<ParticleType>>> WorkerResults;
};
}

217
vtkm/filter/flow/internal/DataSetIntegrator.h
View File

@ -92,8 +92,8 @@ public:
this->TermID.emplace_back(pID);
}
const vtkm::filter::flow::internal::BoundsMap BoundsMap;
const std::unordered_map<vtkm::Id, std::vector<vtkm::Id>> ParticleBlockIDsMap;
vtkm::filter::flow::internal::BoundsMap BoundsMap;
std::unordered_map<vtkm::Id, std::vector<vtkm::Id>> ParticleBlockIDsMap;
ParticleBlockIds InBounds;
ParticleBlockIds OutOfBounds;
@ -102,36 +102,13 @@ public:
std::vector<vtkm::Id> TermIdx;
};
using DSIHelperInfoType =
vtkm::cont::Variant<DSIHelperInfo<vtkm::Particle>, DSIHelperInfo<vtkm::ChargedParticle>>;
template <typename Derived>
template <typename Derived, typename ParticleType>
class DataSetIntegrator
{
public:
using VelocityFieldNameType = std::string;
using ElectroMagneticFieldNameType = std::pair<std::string, std::string>;
protected:
using FieldNameType = vtkm::cont::Variant<VelocityFieldNameType, ElectroMagneticFieldNameType>;
using RType =
vtkm::cont::Variant<vtkm::worklet::flow::ParticleAdvectionResult<vtkm::Particle>,
vtkm::worklet::flow::ParticleAdvectionResult<vtkm::ChargedParticle>,
vtkm::worklet::flow::StreamlineResult<vtkm::Particle>,
vtkm::worklet::flow::StreamlineResult<vtkm::ChargedParticle>>;
public:
DataSetIntegrator(vtkm::Id id,
const FieldNameType& fieldName,
vtkm::filter::flow::IntegrationSolverType solverType,
vtkm::filter::flow::VectorFieldType vecFieldType,
vtkm::filter::flow::FlowResultType resultType)
: FieldName(fieldName)
, Id(id)
DataSetIntegrator(vtkm::Id id, vtkm::filter::flow::IntegrationSolverType solverType)
: Id(id)
, SolverType(solverType)
, VecFieldType(vecFieldType)
, AdvectionResType(resultType)
, Rank(this->Comm.rank())
{
//check that things are valid.
@ -141,57 +118,33 @@ public:
VTKM_CONT void SetCopySeedFlag(bool val) { this->CopySeedArray = val; }
VTKM_CONT
void Advect(DSIHelperInfoType& b,
vtkm::FloatDefault stepSize, //move these to member data(?)
vtkm::Id maxSteps)
void Advect(DSIHelperInfo<ParticleType>& b,
vtkm::FloatDefault stepSize) //move these to member data(?)
{
Derived* inst = static_cast<Derived*>(this);
//Cast the DSIHelperInfo<ParticleType> to the concrete type and call DoAdvect.
b.CastAndCall([&](auto& concrete) { inst->DoAdvect(concrete, stepSize, maxSteps); });
inst->DoAdvect(b, stepSize);
}
template <typename ParticleType>
VTKM_CONT bool GetOutput(vtkm::cont::DataSet& ds) const;
VTKM_CONT bool GetOutput(vtkm::cont::DataSet& dataset) const
{
Derived* inst = static_cast<Derived*>(this);
return inst->GetOutput(dataset);
}
protected:
template <typename ParticleType, template <typename> class ResultType>
VTKM_CONT void UpdateResult(const ResultType<ParticleType>& result,
DSIHelperInfo<ParticleType>& dsiInfo);
VTKM_CONT bool IsParticleAdvectionResult() const
{
return this->AdvectionResType == FlowResultType::PARTICLE_ADVECT_TYPE;
}
VTKM_CONT bool IsStreamlineResult() const
{
return this->AdvectionResType == FlowResultType::STREAMLINE_TYPE;
}
template <typename ParticleType>
VTKM_CONT inline void ClassifyParticles(const vtkm::cont::ArrayHandle<ParticleType>& particles,
DSIHelperInfo<ParticleType>& dsiInfo) const;
//Data members.
vtkm::cont::Variant<VelocityFieldNameType, ElectroMagneticFieldNameType> FieldName;
vtkm::Id Id;
vtkm::filter::flow::IntegrationSolverType SolverType;
vtkm::filter::flow::VectorFieldType VecFieldType;
vtkm::filter::flow::FlowResultType AdvectionResType =
vtkm::filter::flow::FlowResultType::UNKNOWN_TYPE;
vtkmdiy::mpi::communicator Comm = vtkm::cont::EnvironmentTracker::GetCommunicator();
vtkm::Id Rank;
bool CopySeedArray = false;
std::vector<RType> Results;
};
template <typename Derived>
template <typename ParticleType>
VTKM_CONT inline void DataSetIntegrator<Derived>::ClassifyParticles(
template <typename Derived, typename ParticleType>
VTKM_CONT inline void DataSetIntegrator<Derived, ParticleType>::ClassifyParticles(
const vtkm::cont::ArrayHandle<ParticleType>& particles,
DSIHelperInfo<ParticleType>& dsiInfo) const
{
@ -275,146 +228,6 @@ VTKM_CONT inline void DataSetIntegrator<Derived>::ClassifyParticles(
dsiInfo.Validate(n);
}
template <typename Derived>
template <typename ParticleType, template <typename> class ResultType>
VTKM_CONT inline void DataSetIntegrator<Derived>::UpdateResult(
const ResultType<ParticleType>& result,
DSIHelperInfo<ParticleType>& dsiInfo)
{
this->ClassifyParticles(result.Particles, dsiInfo);
if (this->IsParticleAdvectionResult())
{
if (dsiInfo.TermIdx.empty())
return;
using ResType = vtkm::worklet::flow::ParticleAdvectionResult<ParticleType>;
auto indicesAH = vtkm::cont::make_ArrayHandle(dsiInfo.TermIdx, vtkm::CopyFlag::Off);
auto termPerm = vtkm::cont::make_ArrayHandlePermutation(indicesAH, result.Particles);
vtkm::cont::ArrayHandle<ParticleType> termParticles;
vtkm::cont::Algorithm::Copy(termPerm, termParticles);
ResType termRes(termParticles);
this->Results.emplace_back(termRes);
}
else if (this->IsStreamlineResult())
this->Results.emplace_back(result);
}
template <typename Derived>
template <typename ParticleType>
VTKM_CONT inline bool DataSetIntegrator<Derived>::GetOutput(vtkm::cont::DataSet& ds) const
{
std::size_t nResults = this->Results.size();
if (nResults == 0)
return false;
if (this->IsParticleAdvectionResult())
{
using ResType = vtkm::worklet::flow::ParticleAdvectionResult<ParticleType>;
std::vector<vtkm::cont::ArrayHandle<ParticleType>> allParticles;
allParticles.reserve(nResults);
for (const auto& vres : this->Results)
allParticles.emplace_back(vres.template Get<ResType>().Particles);
vtkm::cont::ArrayHandle<vtkm::Vec3f> pts;
vtkm::cont::ParticleArrayCopy(allParticles, pts);
vtkm::Id numPoints = pts.GetNumberOfValues();
if (numPoints > 0)
{
//Create coordinate system and vertex cell set.
ds.AddCoordinateSystem(vtkm::cont::CoordinateSystem("coordinates", pts));
vtkm::cont::CellSetSingleType<> cells;
vtkm::cont::ArrayHandleIndex conn(numPoints);
vtkm::cont::ArrayHandle<vtkm::Id> connectivity;
vtkm::cont::ArrayCopy(conn, connectivity);
cells.Fill(numPoints, vtkm::CELL_SHAPE_VERTEX, 1, connectivity);
ds.SetCellSet(cells);
}
}
else if (this->IsStreamlineResult())
{
using ResType = vtkm::worklet::flow::StreamlineResult<ParticleType>;
//Easy case with one result.
if (nResults == 1)
{
const auto& res = this->Results[0].template Get<ResType>();
ds.AddCoordinateSystem(vtkm::cont::CoordinateSystem("coordinates", res.Positions));
ds.SetCellSet(res.PolyLines);
}
else
{
std::vector<vtkm::Id> posOffsets(nResults, 0);
vtkm::Id totalNumCells = 0, totalNumPts = 0;
for (std::size_t i = 0; i < nResults; i++)
{
const auto& res = this->Results[i].template Get<ResType>();
if (i == 0)
posOffsets[i] = 0;
else
posOffsets[i] = totalNumPts;
totalNumPts += res.Positions.GetNumberOfValues();
totalNumCells += res.PolyLines.GetNumberOfCells();
}
//Append all the points together.
vtkm::cont::ArrayHandle<vtkm::Vec3f> appendPts;
appendPts.Allocate(totalNumPts);
for (std::size_t i = 0; i < nResults; i++)
{
const auto& res = this->Results[i].template Get<ResType>();
// copy all values into appendPts starting at offset.
vtkm::cont::Algorithm::CopySubRange(
res.Positions, 0, res.Positions.GetNumberOfValues(), appendPts, posOffsets[i]);
}
ds.AddCoordinateSystem(vtkm::cont::CoordinateSystem("coordinates", appendPts));
//Create polylines.
std::vector<vtkm::Id> numPtsPerCell(static_cast<std::size_t>(totalNumCells));
std::size_t off = 0;
for (std::size_t i = 0; i < nResults; i++)
{
const auto& res = this->Results[i].template Get<ResType>();
vtkm::Id nCells = res.PolyLines.GetNumberOfCells();
for (vtkm::Id j = 0; j < nCells; j++)
numPtsPerCell[off++] = static_cast<vtkm::Id>(res.PolyLines.GetNumberOfPointsInCell(j));
}
auto numPointsPerCellArray = vtkm::cont::make_ArrayHandle(numPtsPerCell, vtkm::CopyFlag::Off);
vtkm::cont::ArrayHandle<vtkm::Id> cellIndex;
vtkm::Id connectivityLen =
vtkm::cont::Algorithm::ScanExclusive(numPointsPerCellArray, cellIndex);
vtkm::cont::ArrayHandleIndex connCount(connectivityLen);
vtkm::cont::ArrayHandle<vtkm::Id> connectivity;
vtkm::cont::ArrayCopy(connCount, connectivity);
vtkm::cont::ArrayHandle<vtkm::UInt8> cellTypes;
auto polyLineShape = vtkm::cont::make_ArrayHandleConstant<vtkm::UInt8>(
vtkm::CELL_SHAPE_POLY_LINE, totalNumCells);
vtkm::cont::ArrayCopy(polyLineShape, cellTypes);
auto offsets = vtkm::cont::ConvertNumComponentsToOffsets(numPointsPerCellArray);
vtkm::cont::CellSetExplicit<> polyLines;
polyLines.Fill(totalNumPts, cellTypes, connectivity, offsets);
ds.SetCellSet(polyLines);
}
}
else
{
throw vtkm::cont::ErrorFilterExecution("Unsupported ParticleAdvectionResultType");
}
return true;
}
}
}
}

336
vtkm/filter/flow/internal/DataSetIntegratorSteadyState.h
View File

@ -22,252 +22,150 @@ namespace flow
{
namespace internal
{
class DataSetIntegratorSteadyState
: public vtkm::filter::flow::internal::DataSetIntegrator<DataSetIntegratorSteadyState>
{
public:
DataSetIntegratorSteadyState(const vtkm::cont::DataSet& ds,
vtkm::Id id,
const FieldNameType& fieldName,
vtkm::filter::flow::IntegrationSolverType solverType,
vtkm::filter::flow::VectorFieldType vecFieldType,
vtkm::filter::flow::FlowResultType resultType)
: vtkm::filter::flow::internal::DataSetIntegrator<DataSetIntegratorSteadyState>(id,
fieldName,
solverType,
vecFieldType,
resultType)
, DataSet(ds)
{
}
VTKM_CONT inline void DoAdvect(DSIHelperInfo<vtkm::Particle>& b,
vtkm::FloatDefault stepSize,
vtkm::Id maxSteps);
VTKM_CONT inline void DoAdvect(DSIHelperInfo<vtkm::ChargedParticle>& b,
vtkm::FloatDefault stepSize,
vtkm::Id maxSteps);
protected:
template <typename ArrayType>
VTKM_CONT void GetVelocityField(
vtkm::worklet::flow::VelocityField<ArrayType>& velocityField) const
{
if (this->FieldName.IsType<VelocityFieldNameType>())
{
const auto& fieldNm = this->FieldName.Get<VelocityFieldNameType>();
auto assoc = this->DataSet.GetField(fieldNm).GetAssociation();
ArrayType arr;
vtkm::cont::ArrayCopyShallowIfPossible(this->DataSet.GetField(fieldNm).GetData(), arr);
velocityField = vtkm::worklet::flow::VelocityField<ArrayType>(arr, assoc);
}
else
throw vtkm::cont::ErrorFilterExecution("Velocity field vector type not available");
}
template <typename ArrayType>
VTKM_CONT void GetElectroMagneticField(
vtkm::worklet::flow::ElectroMagneticField<ArrayType>& ebField) const
{
if (this->FieldName.IsType<ElectroMagneticFieldNameType>())
{
const auto& fieldNm = this->FieldName.Get<ElectroMagneticFieldNameType>();
const auto& electric = fieldNm.first;
const auto& magnetic = fieldNm.second;
auto eAssoc = this->DataSet.GetField(electric).GetAssociation();
auto bAssoc = this->DataSet.GetField(magnetic).GetAssociation();
if (eAssoc != bAssoc)
{
throw vtkm::cont::ErrorFilterExecution("E and B field need to have same association");
}
ArrayType eField, bField;
vtkm::cont::ArrayCopyShallowIfPossible(this->DataSet.GetField(electric).GetData(), eField);
vtkm::cont::ArrayCopyShallowIfPossible(this->DataSet.GetField(magnetic).GetData(), bField);
ebField = vtkm::worklet::flow::ElectroMagneticField<ArrayType>(eField, bField, eAssoc);
}
else
throw vtkm::cont::ErrorFilterExecution("Electromagnetic field vector type not available");
}
private:
vtkm::cont::DataSet DataSet;
};
namespace internal
{
using ArrayType = vtkm::cont::ArrayHandle<vtkm::Vec3f>;
using VelocityFieldType = vtkm::worklet::flow::VelocityField<ArrayType>;
using VelocityEvalType = vtkm::worklet::flow::GridEvaluator<VelocityFieldType>;
using EBFieldType = vtkm::worklet::flow::ElectroMagneticField<ArrayType>;
using EBEvalType = vtkm::worklet::flow::GridEvaluator<EBFieldType>;
//template <typename GridEvalType, typename ParticleType>
//class AdvectHelper;
template <typename FieldType, typename ParticleType>
class AdvectHelper //<FieldType, ParticleType>
namespace detail
{
template <typename ParticleType,
typename FieldType,
typename TerminationType,
typename AnalysisType>
class AdvectHelperSteadyState
{
public:
using WorkletType = vtkm::worklet::flow::ParticleAdvection;
using SteadyStateGridEvalType = vtkm::worklet::flow::GridEvaluator<FieldType>;
static void Advect(const FieldType& vecField,
const vtkm::cont::DataSet& ds,
vtkm::cont::ArrayHandle<ParticleType>& seedArray,
template <template <typename> class SolverType>
static void DoAdvect(vtkm::cont::ArrayHandle<ParticleType>& seedArray,
const FieldType& field,
const vtkm::cont::DataSet& dataset,
const TerminationType& termination,
vtkm::FloatDefault stepSize,
vtkm::Id maxSteps,
const IntegrationSolverType& solverType,
vtkm::worklet::flow::ParticleAdvectionResult<ParticleType>& result)
{
if (solverType == IntegrationSolverType::RK4_TYPE)
{
DoAdvect<vtkm::worklet::flow::ParticleAdvection,
vtkm::worklet::flow::ParticleAdvectionResult,
vtkm::worklet::flow::RK4Integrator>(
vecField, ds, seedArray, stepSize, maxSteps, result);
}
else if (solverType == IntegrationSolverType::EULER_TYPE)
{
DoAdvect<vtkm::worklet::flow::ParticleAdvection,
vtkm::worklet::flow::ParticleAdvectionResult,
vtkm::worklet::flow::EulerIntegrator>(
vecField, ds, seedArray, stepSize, maxSteps, result);
}
else
throw vtkm::cont::ErrorFilterExecution("Unsupported Integrator type");
}
static void Advect(const FieldType& vecField,
const vtkm::cont::DataSet& ds,
vtkm::cont::ArrayHandle<ParticleType>& seedArray,
vtkm::FloatDefault stepSize,
vtkm::Id maxSteps,
const IntegrationSolverType& solverType,
vtkm::worklet::flow::StreamlineResult<ParticleType>& result)
{
if (solverType == IntegrationSolverType::RK4_TYPE)
{
DoAdvect<vtkm::worklet::flow::Streamline,
vtkm::worklet::flow::StreamlineResult,
vtkm::worklet::flow::RK4Integrator>(
vecField, ds, seedArray, stepSize, maxSteps, result);
}
else if (solverType == IntegrationSolverType::EULER_TYPE)
{
DoAdvect<vtkm::worklet::flow::Streamline,
vtkm::worklet::flow::StreamlineResult,
vtkm::worklet::flow::EulerIntegrator>(
vecField, ds, seedArray, stepSize, maxSteps, result);
}
else
throw vtkm::cont::ErrorFilterExecution("Unsupported Integrator type");
}
template <typename WorkletType,
template <typename>
class ResultType,
template <typename>
class SolverType>
static void DoAdvect(const FieldType& vecField,
const vtkm::cont::DataSet& ds,
vtkm::cont::ArrayHandle<ParticleType>& seedArray,
vtkm::FloatDefault stepSize,
vtkm::Id maxSteps,
ResultType<ParticleType>& result)
AnalysisType& analysis)
{
using StepperType =
vtkm::worklet::flow::Stepper<SolverType<SteadyStateGridEvalType>, SteadyStateGridEvalType>;
SteadyStateGridEvalType eval(dataset, field);
StepperType stepper(eval, stepSize);
WorkletType worklet;
SteadyStateGridEvalType eval(ds, vecField);
StepperType stepper(eval, stepSize);
result = worklet.Run(stepper, seedArray, maxSteps);
worklet.Run(stepper, seedArray, termination, analysis);
}
static void Advect(vtkm::cont::ArrayHandle<ParticleType>& seedArray,
const FieldType& field,
const vtkm::cont::DataSet& dataset,
const TerminationType& termination,
const IntegrationSolverType& solverType,
vtkm::FloatDefault stepSize,
AnalysisType& analysis)
{
if (solverType == IntegrationSolverType::RK4_TYPE)
{
DoAdvect<vtkm::worklet::flow::RK4Integrator>(
seedArray, field, dataset, termination, stepSize, analysis);
}
else if (solverType == IntegrationSolverType::EULER_TYPE)
{
DoAdvect<vtkm::worklet::flow::EulerIntegrator>(
seedArray, field, dataset, termination, stepSize, analysis);
}
else
throw vtkm::cont::ErrorFilterExecution("Unsupported Integrator type");
}
};
} // namespace detail
template <typename ParticleType,
typename FieldType,
typename TerminationType,
typename AnalysisType>
class DataSetIntegratorSteadyState
: public vtkm::filter::flow::internal::DataSetIntegrator<
DataSetIntegratorSteadyState<ParticleType, FieldType, TerminationType, AnalysisType>,
ParticleType>
{
public:
using BaseType = vtkm::filter::flow::internal::DataSetIntegrator<
DataSetIntegratorSteadyState<ParticleType, FieldType, TerminationType, AnalysisType>,
ParticleType>;
using PType = ParticleType;
using FType = FieldType;
using TType = TerminationType;
using AType = AnalysisType;
DataSetIntegratorSteadyState(vtkm::Id id,
const FieldType& field,
const vtkm::cont::DataSet& dataset,
vtkm::filter::flow::IntegrationSolverType solverType,
const TerminationType& termination,
const AnalysisType& analysis)
: BaseType(id, solverType)
, Field(field)
, Dataset(dataset)
, Termination(termination)
, Analysis(analysis)
{
}
VTKM_CONT inline void DataSetIntegratorSteadyState::DoAdvect(DSIHelperInfo<vtkm::Particle>& b,
vtkm::FloatDefault stepSize,
vtkm::Id maxSteps)
VTKM_CONT inline void DoAdvect(vtkm::filter::flow::internal::DSIHelperInfo<ParticleType>& block,
vtkm::FloatDefault stepSize)
{
using ArrayType = vtkm::cont::ArrayHandle<vtkm::Vec3f>;
auto copyFlag = (this->CopySeedArray ? vtkm::CopyFlag::On : vtkm::CopyFlag::Off);
auto seedArray = vtkm::cont::make_ArrayHandle(b.Particles, copyFlag);
auto seedArray = vtkm::cont::make_ArrayHandle(block.Particles, copyFlag);
if (this->VecFieldType == VectorFieldType::VELOCITY_FIELD_TYPE)
{
using FieldType = vtkm::worklet::flow::VelocityField<ArrayType>;
FieldType vecField;
this->GetVelocityField(vecField);
using AdvectionHelper =
detail::AdvectHelperSteadyState<ParticleType, FieldType, TerminationType, AnalysisType>;
AnalysisType analysis;
analysis.UseAsTemplate(this->Analysis);
using AHType = internal::AdvectHelper<internal::VelocityFieldType, vtkm::Particle>;
AdvectionHelper::Advect(seedArray,
this->Field,
this->Dataset,
this->Termination,
this->SolverType,
stepSize,
analysis);
if (this->IsParticleAdvectionResult())
{
vtkm::worklet::flow::ParticleAdvectionResult<vtkm::Particle> result;
AHType::Advect(
vecField, this->DataSet, seedArray, stepSize, maxSteps, this->SolverType, result);
this->UpdateResult(result, b);
this->UpdateResult(analysis, block);
}
else if (this->IsStreamlineResult())
VTKM_CONT void UpdateResult(AnalysisType& analysis,
vtkm::filter::flow::internal::DSIHelperInfo<ParticleType>& dsiInfo)
{
vtkm::worklet::flow::StreamlineResult<vtkm::Particle> result;
AHType::Advect(
vecField, this->DataSet, seedArray, stepSize, maxSteps, this->SolverType, result);
this->UpdateResult(result, b);
this->ClassifyParticles(analysis.Particles, dsiInfo);
if (std::is_same<AnalysisType, vtkm::worklet::flow::NoAnalysis<ParticleType>>::value)
{
if (dsiInfo.TermIdx.empty())
return;
auto indicesAH = vtkm::cont::make_ArrayHandle(dsiInfo.TermIdx, vtkm::CopyFlag::Off);
auto termPerm = vtkm::cont::make_ArrayHandlePermutation(indicesAH, analysis.Particles);
vtkm::cont::ArrayHandle<ParticleType> termParticles;
vtkm::cont::Algorithm::Copy(termPerm, termParticles);
analysis.FinalizeAnalysis(termParticles);
this->Analyses.emplace_back(analysis);
}
else
throw vtkm::cont::ErrorFilterExecution("Unsupported result type");
}
else
throw vtkm::cont::ErrorFilterExecution("Unsupported vector field type");
}
VTKM_CONT inline void DataSetIntegratorSteadyState::DoAdvect(
DSIHelperInfo<vtkm::ChargedParticle>& b,
vtkm::FloatDefault stepSize,
vtkm::Id maxSteps)
{
using ArrayType = vtkm::cont::ArrayHandle<vtkm::Vec3f>;
this->Analyses.emplace_back(analysis);
}
}
auto copyFlag = (this->CopySeedArray ? vtkm::CopyFlag::On : vtkm::CopyFlag::Off);
auto seedArray = vtkm::cont::make_ArrayHandle(b.Particles, copyFlag);
if (this->VecFieldType == VectorFieldType::ELECTRO_MAGNETIC_FIELD_TYPE)
VTKM_CONT bool GetOutput(vtkm::cont::DataSet& ds) const
{
using FieldType = vtkm::worklet::flow::ElectroMagneticField<ArrayType>;
FieldType ebField;
this->GetElectroMagneticField(ebField);
std::size_t nAnalyses = this->Analyses.size();
if (nAnalyses == 0)
return false;
return AnalysisType::MakeDataSet(ds, this->Analyses);
}
using AHType = internal::AdvectHelper<internal::EBFieldType, vtkm::ChargedParticle>;
if (this->IsParticleAdvectionResult())
{
vtkm::worklet::flow::ParticleAdvectionResult<vtkm::ChargedParticle> result;
AHType::Advect(
ebField, this->DataSet, seedArray, stepSize, maxSteps, this->SolverType, result);
this->UpdateResult(result, b);
}
else if (this->IsStreamlineResult())
{
vtkm::worklet::flow::StreamlineResult<vtkm::ChargedParticle> result;
AHType::Advect(
ebField, this->DataSet, seedArray, stepSize, maxSteps, this->SolverType, result);
this->UpdateResult(result, b);
}
else
throw vtkm::cont::ErrorFilterExecution("Unsupported result type");
}
else
throw vtkm::cont::ErrorFilterExecution("Unsupported vector field type");
}
private:
FieldType Field;
vtkm::cont::DataSet Dataset;
TerminationType Termination;
// Used as a template to initialize successive analysis objects.
AnalysisType Analysis;
std::vector<AnalysisType> Analyses;
};
}
}

351
vtkm/filter/flow/internal/DataSetIntegratorUnsteadyState.h
View File

@ -23,242 +23,173 @@ namespace flow
{
namespace internal
{
class DataSetIntegratorUnsteadyState
: public vtkm::filter::flow::internal::DataSetIntegrator<DataSetIntegratorUnsteadyState>
namespace detail
{
template <typename ParticleType,
typename FieldType,
typename TerminationType,
typename AnalysisType>
class AdvectHelperUnsteadyState
{
public:
DataSetIntegratorUnsteadyState(const vtkm::cont::DataSet& ds1,
using WorkletType = vtkm::worklet::flow::ParticleAdvection;
using UnsteadyStateGridEvalType = vtkm::worklet::flow::TemporalGridEvaluator<FieldType>;
template <template <typename> class SolverType>
static void DoAdvect(vtkm::cont::ArrayHandle<ParticleType>& seedArray,
const FieldType& field1,
const vtkm::cont::DataSet& ds1,
vtkm::FloatDefault t1,
const FieldType& field2,
const vtkm::cont::DataSet& ds2,
vtkm::FloatDefault t2,
const TerminationType& termination,
vtkm::FloatDefault stepSize,
AnalysisType& analysis)
{
using StepperType = vtkm::worklet::flow::Stepper<SolverType<UnsteadyStateGridEvalType>,
UnsteadyStateGridEvalType>;
WorkletType worklet;
UnsteadyStateGridEvalType eval(ds1, t1, field1, ds2, t2, field2);
StepperType stepper(eval, stepSize);
worklet.Run(stepper, seedArray, termination, analysis);
}
static void Advect(vtkm::cont::ArrayHandle<ParticleType>& seedArray,
const FieldType& field1,
const vtkm::cont::DataSet& ds1,
vtkm::FloatDefault t1,
const FieldType& field2,
const vtkm::cont::DataSet& ds2,
vtkm::FloatDefault t2,
const TerminationType& termination,
const IntegrationSolverType& solverType,
vtkm::FloatDefault stepSize,
AnalysisType& analysis)
{
if (solverType == IntegrationSolverType::RK4_TYPE)
{
DoAdvect<vtkm::worklet::flow::RK4Integrator>(
seedArray, field1, ds1, t1, field2, ds2, t2, termination, stepSize, analysis);
}
else if (solverType == IntegrationSolverType::EULER_TYPE)
{
DoAdvect<vtkm::worklet::flow::EulerIntegrator>(
seedArray, field1, ds1, t1, field2, ds2, t2, termination, stepSize, analysis);
}
else
throw vtkm::cont::ErrorFilterExecution("Unsupported Integrator type");
}
};
} //namespace detail
template <typename ParticleType,
typename FieldType,
typename TerminationType,
typename AnalysisType>
class DataSetIntegratorUnsteadyState
: public vtkm::filter::flow::internal::DataSetIntegrator<
DataSetIntegratorUnsteadyState<ParticleType, FieldType, TerminationType, AnalysisType>,
ParticleType>
{
public:
using BaseType = vtkm::filter::flow::internal::DataSetIntegrator<
DataSetIntegratorUnsteadyState<ParticleType, FieldType, TerminationType, AnalysisType>,
ParticleType>;
using PType = ParticleType;
using FType = FieldType;
using TType = TerminationType;
using AType = AnalysisType;
DataSetIntegratorUnsteadyState(vtkm::Id id,
const FieldType& field1,
const FieldType& field2,
const vtkm::cont::DataSet& ds1,
const vtkm::cont::DataSet& ds2,
vtkm::FloatDefault t1,
vtkm::FloatDefault t2,
vtkm::Id id,
const vtkm::filter::flow::internal::DataSetIntegrator<
DataSetIntegratorUnsteadyState>::FieldNameType& fieldName,
vtkm::filter::flow::IntegrationSolverType solverType,
vtkm::filter::flow::VectorFieldType vecFieldType,
vtkm::filter::flow::FlowResultType resultType)
: vtkm::filter::flow::internal::DataSetIntegrator<DataSetIntegratorUnsteadyState>(id,
fieldName,
solverType,
vecFieldType,
resultType)
const TerminationType& termination,
const AnalysisType& analysis)
: BaseType(id, solverType)
, Field1(field1)
, Field2(field2)
, DataSet1(ds1)
, DataSet2(ds2)
, Time1(t1)
, Time2(t2)
, Termination(termination)
, Analysis(analysis)
{
}
VTKM_CONT inline void DoAdvect(DSIHelperInfo<vtkm::Particle>& b,
vtkm::FloatDefault stepSize,
vtkm::Id maxSteps);
VTKM_CONT inline void DoAdvect(DSIHelperInfo<vtkm::ChargedParticle>& b,
vtkm::FloatDefault stepSize,
vtkm::Id maxSteps);
protected:
template <typename ArrayType>
VTKM_CONT void GetVelocityFields(
vtkm::worklet::flow::VelocityField<ArrayType>& velocityField1,
vtkm::worklet::flow::VelocityField<ArrayType>& velocityField2) const
VTKM_CONT inline void DoAdvect(vtkm::filter::flow::internal::DSIHelperInfo<ParticleType>& block,
vtkm::FloatDefault stepSize)
{
if (this->FieldName.GetIndex() == this->FieldName.GetIndexOf<VelocityFieldNameType>())
auto copyFlag = (this->CopySeedArray ? vtkm::CopyFlag::On : vtkm::CopyFlag::Off);
auto seedArray = vtkm::cont::make_ArrayHandle(block.Particles, copyFlag);
using AdvectionHelper =
detail::AdvectHelperUnsteadyState<ParticleType, FieldType, TerminationType, AnalysisType>;
AnalysisType analysis;
analysis.UseAsTemplate(this->Analysis);
AdvectionHelper::Advect(seedArray,
this->Field1,
this->DataSet1,
this->Time1,
this->Field2,
this->DataSet2,
this->Time2,
this->Termination,
this->SolverType,
stepSize,
analysis);
this->UpdateResult(analysis, block);
}
VTKM_CONT void UpdateResult(AnalysisType& analysis,
vtkm::filter::flow::internal::DSIHelperInfo<ParticleType>& dsiInfo)
{
const auto& fieldNm = this->FieldName.Get<VelocityFieldNameType>();
auto assoc = this->DataSet1.GetField(fieldNm).GetAssociation();
if (assoc != this->DataSet2.GetField(fieldNm).GetAssociation())
throw vtkm::cont::ErrorFilterExecution(
"Unsteady state velocity fields have differnt associations");
ArrayType arr1, arr2;
vtkm::cont::ArrayCopyShallowIfPossible(this->DataSet1.GetField(fieldNm).GetData(), arr1);
vtkm::cont::ArrayCopyShallowIfPossible(this->DataSet2.GetField(fieldNm).GetData(), arr2);
velocityField1 = vtkm::worklet::flow::VelocityField<ArrayType>(arr1, assoc);
velocityField2 = vtkm::worklet::flow::VelocityField<ArrayType>(arr2, assoc);
this->ClassifyParticles(analysis.Particles, dsiInfo);
if (std::is_same<AnalysisType, vtkm::worklet::flow::NoAnalysis<ParticleType>>::value)
{
if (dsiInfo.TermIdx.empty())
return;
auto indicesAH = vtkm::cont::make_ArrayHandle(dsiInfo.TermIdx, vtkm::CopyFlag::Off);
auto termPerm = vtkm::cont::make_ArrayHandlePermutation(indicesAH, analysis.Particles);
vtkm::cont::ArrayHandle<ParticleType> termParticles;
vtkm::cont::Algorithm::Copy(termPerm, termParticles);
analysis.FinalizeAnalysis(termParticles);
this->Analyses.emplace_back(analysis);
}
else
throw vtkm::cont::ErrorFilterExecution("Velocity field vector type not available");
{
this->Analyses.emplace_back(analysis);
}
}
VTKM_CONT bool GetOutput(vtkm::cont::DataSet& ds) const
{
std::size_t nAnalyses = this->Analyses.size();
if (nAnalyses == 0)
return false;
return AnalysisType::MakeDataSet(ds, this->Analyses);
}
private:
FieldType Field1;
FieldType Field2;
vtkm::cont::DataSet DataSet1;
vtkm::cont::DataSet DataSet2;
vtkm::FloatDefault Time1;
vtkm::FloatDefault Time2;
TerminationType Termination;
AnalysisType Analysis;
std::vector<AnalysisType> Analyses;
};
namespace internal
{
using ArrayType = vtkm::cont::ArrayHandle<vtkm::Vec3f>;
using VelocityFieldType = vtkm::worklet::flow::VelocityField<ArrayType>;
using UnsteadyStateGridEvalType = vtkm::worklet::flow::TemporalGridEvaluator<VelocityFieldType>;
template <typename GridEvalType, typename ParticleType>
class AdvectHelper;
template <typename ParticleType>
class AdvectHelper<UnsteadyStateGridEvalType, ParticleType>
{
public:
static void Advect(const VelocityFieldType& velField1,
const vtkm::cont::DataSet& ds1,
vtkm::FloatDefault t1,
const VelocityFieldType& velField2,
const vtkm::cont::DataSet& ds2,
vtkm::FloatDefault t2,
vtkm::cont::ArrayHandle<ParticleType>& seedArray,
vtkm::FloatDefault stepSize,
vtkm::Id maxSteps,
const IntegrationSolverType& solverType,
vtkm::worklet::flow::ParticleAdvectionResult<ParticleType>& result)
{
if (solverType == IntegrationSolverType::RK4_TYPE)
{
DoAdvect<vtkm::worklet::flow::ParticleAdvection,
vtkm::worklet::flow::ParticleAdvectionResult,
vtkm::worklet::flow::RK4Integrator>(
velField1, ds1, t1, velField2, ds2, t2, seedArray, stepSize, maxSteps, result);
}
else if (solverType == IntegrationSolverType::EULER_TYPE)
{
DoAdvect<vtkm::worklet::flow::ParticleAdvection,
vtkm::worklet::flow::ParticleAdvectionResult,
vtkm::worklet::flow::EulerIntegrator>(
velField1, ds1, t1, velField2, ds2, t2, seedArray, stepSize, maxSteps, result);
}
else
throw vtkm::cont::ErrorFilterExecution("Unsupported Integrator type");
}
static void Advect(const VelocityFieldType& velField1,
const vtkm::cont::DataSet& ds1,
vtkm::FloatDefault t1,
const VelocityFieldType& velField2,
const vtkm::cont::DataSet& ds2,
vtkm::FloatDefault t2,
vtkm::cont::ArrayHandle<ParticleType>& seedArray,
vtkm::FloatDefault stepSize,
vtkm::Id maxSteps,
const IntegrationSolverType& solverType,
vtkm::worklet::flow::StreamlineResult<ParticleType>& result)
{
if (solverType == IntegrationSolverType::RK4_TYPE)
{
DoAdvect<vtkm::worklet::flow::Streamline,
vtkm::worklet::flow::StreamlineResult,
vtkm::worklet::flow::RK4Integrator>(
velField1, ds1, t1, velField2, ds2, t2, seedArray, stepSize, maxSteps, result);
}
else if (solverType == IntegrationSolverType::EULER_TYPE)
{
DoAdvect<vtkm::worklet::flow::Streamline,
vtkm::worklet::flow::StreamlineResult,
vtkm::worklet::flow::EulerIntegrator>(
velField1, ds1, t1, velField2, ds2, t2, seedArray, stepSize, maxSteps, result);
}
else
throw vtkm::cont::ErrorFilterExecution("Unsupported Integrator type");
}
template <typename WorkletType,
template <typename>
class ResultType,
template <typename>
class SolverType>
static void DoAdvect(const VelocityFieldType& velField1,
const vtkm::cont::DataSet& ds1,
vtkm::FloatDefault t1,
const VelocityFieldType& velField2,
const vtkm::cont::DataSet& ds2,
vtkm::FloatDefault t2,
vtkm::cont::ArrayHandle<ParticleType>& seedArray,
vtkm::FloatDefault stepSize,
vtkm::Id maxSteps,
ResultType<ParticleType>& result)
{
using StepperType = vtkm::worklet::flow::Stepper<SolverType<UnsteadyStateGridEvalType>,
UnsteadyStateGridEvalType>;
WorkletType worklet;
UnsteadyStateGridEvalType eval(ds1, t1, velField1, ds2, t2, velField2);
StepperType stepper(eval, stepSize);
result = worklet.Run(stepper, seedArray, maxSteps);
}
};
}
VTKM_CONT inline void DataSetIntegratorUnsteadyState::DoAdvect(DSIHelperInfo<vtkm::Particle>& b,
vtkm::FloatDefault stepSize,
vtkm::Id maxSteps)
{
using ArrayType = vtkm::cont::ArrayHandle<vtkm::Vec3f>;
auto copyFlag = (this->CopySeedArray ? vtkm::CopyFlag::On : vtkm::CopyFlag::Off);
auto seedArray = vtkm::cont::make_ArrayHandle(b.Particles, copyFlag);
using AHType = internal::AdvectHelper<internal::UnsteadyStateGridEvalType, vtkm::Particle>;
if (this->VecFieldType == VectorFieldType::VELOCITY_FIELD_TYPE)
{
using FieldType = vtkm::worklet::flow::VelocityField<ArrayType>;
FieldType velField1, velField2;
this->GetVelocityFields(velField1, velField2);
if (this->IsParticleAdvectionResult())
{
vtkm::worklet::flow::ParticleAdvectionResult<vtkm::Particle> result;
AHType::Advect(velField1,
this->DataSet1,
this->Time1,
velField2,
this->DataSet2,
this->Time2,
seedArray,
stepSize,
maxSteps,
this->SolverType,
result);
this->UpdateResult(result, b);
}
else if (this->IsStreamlineResult())
{
vtkm::worklet::flow::StreamlineResult<vtkm::Particle> result;
AHType::Advect(velField1,
this->DataSet1,
this->Time1,
velField2,
this->DataSet2,
this->Time2,
seedArray,
stepSize,
maxSteps,
this->SolverType,
result);
this->UpdateResult(result, b);
}
else
throw vtkm::cont::ErrorFilterExecution("Unsupported result type");
}
else
throw vtkm::cont::ErrorFilterExecution("Unsupported vector field type");
}
VTKM_CONT inline void DataSetIntegratorUnsteadyState::DoAdvect(
DSIHelperInfo<vtkm::ChargedParticle>& vtkmNotUsed(b),
vtkm::FloatDefault vtkmNotUsed(stepSize),
vtkm::Id vtkmNotUsed(maxSteps))
{
throw vtkm::cont::ErrorFilterExecution(
"Unsupported operation : charged particles and electromagnetic fielfs currently only supported "
"for steady state");
}
}
}
}

86
vtkm/filter/flow/internal/ParticleAdvector.h
View File

@ -29,94 +29,48 @@ template <typename DSIType>
class ParticleAdvector
{
public:
using ParticleType = typename DSIType::PType;
ParticleAdvector(const vtkm::filter::flow::internal::BoundsMap& bm,
const std::vector<DSIType>& blocks,
const bool& useThreaded,
const bool& useAsyncComm,
const vtkm::filter::flow::FlowResultType& parType)
const bool& useAsyncComm)
: Blocks(blocks)
, BoundsMap(bm)
, ResultType(parType)
, UseAsynchronousCommunication(useAsyncComm)
, UseThreadedAlgorithm(useThreaded)
, UseAsynchronousCommunication(useAsyncComm)
{
}
vtkm::cont::PartitionedDataSet Execute(vtkm::Id numSteps,
vtkm::FloatDefault stepSize,
const vtkm::cont::UnknownArrayHandle& seeds)
{
using ParticleTypeList = vtkm::List<vtkm::Particle, vtkm::ChargedParticle>;
vtkm::cont::PartitionedDataSet result;
seeds.CastAndCallForTypes<ParticleTypeList, VTKM_DEFAULT_STORAGE_LIST>(
[&](const auto& concreteSeeds) {
result = this->Execute(numSteps, stepSize, concreteSeeds);
});
return result;
}
private:
template <typename AlgorithmType, typename ParticleType>
vtkm::cont::PartitionedDataSet RunAlgo(vtkm::Id numSteps,
vtkm::FloatDefault stepSize,
const vtkm::cont::ArrayHandle<ParticleType>& seeds)
{
AlgorithmType algo(this->BoundsMap, this->Blocks, this->UseAsynchronousCommunication);
algo.Execute(numSteps, stepSize, seeds);
return algo.GetOutput();
}
template <typename ParticleType>
vtkm::cont::PartitionedDataSet Execute(vtkm::Id numSteps,
vtkm::FloatDefault stepSize,
const vtkm::cont::ArrayHandle<ParticleType>& seeds)
vtkm::cont::PartitionedDataSet Execute(const vtkm::cont::ArrayHandle<ParticleType>& seeds,
vtkm::FloatDefault stepSize)
{
if (!this->UseThreadedAlgorithm)
{
if (this->ResultType == vtkm::filter::flow::FlowResultType::PARTICLE_ADVECT_TYPE)
{
using AlgorithmType = vtkm::filter::flow::internal::
AdvectAlgorithm<DSIType, vtkm::worklet::flow::ParticleAdvectionResult, ParticleType>;
return this->RunAlgo<AlgorithmType, ParticleType>(numSteps, stepSize, seeds);
using AlgorithmType = vtkm::filter::flow::internal::AdvectAlgorithm<DSIType>;
return this->RunAlgo<AlgorithmType>(seeds, stepSize);
}
else
{
using AlgorithmType = vtkm::filter::flow::internal::
AdvectAlgorithm<DSIType, vtkm::worklet::flow::StreamlineResult, ParticleType>;
using AlgorithmType = vtkm::filter::flow::internal::AdvectAlgorithmThreaded<DSIType>;
return this->RunAlgo<AlgorithmType>(seeds, stepSize);
}
}
return this->RunAlgo<AlgorithmType, ParticleType>(numSteps, stepSize, seeds);
}
}
else
private:
template <typename AlgorithmType>
vtkm::cont::PartitionedDataSet RunAlgo(const vtkm::cont::ArrayHandle<ParticleType>& seeds,
vtkm::FloatDefault stepSize)
{
if (this->ResultType == vtkm::filter::flow::FlowResultType::PARTICLE_ADVECT_TYPE)
{
using AlgorithmType = vtkm::filter::flow::internal::AdvectAlgorithmThreaded<
DSIType,
vtkm::worklet::flow::ParticleAdvectionResult,
ParticleType>;
return this->RunAlgo<AlgorithmType, ParticleType>(numSteps, stepSize, seeds);
AlgorithmType algo(this->BoundsMap, this->Blocks, this->UseAsynchronousCommunication);
algo.Execute(seeds, stepSize);
return algo.GetOutput();
}
else
{
using AlgorithmType = vtkm::filter::flow::internal::
AdvectAlgorithmThreaded<DSIType, vtkm::worklet::flow::StreamlineResult, ParticleType>;
return this->RunAlgo<AlgorithmType, ParticleType>(numSteps, stepSize, seeds);
}
}
}
std::vector<DSIType> Blocks;
vtkm::filter::flow::internal::BoundsMap BoundsMap;
FlowResultType ResultType;
bool UseAsynchronousCommunication = true;
bool UseThreadedAlgorithm;
bool UseAsynchronousCommunication = true;
};
}

4
vtkm/filter/flow/testing/UnitTestStreamlineFilter.cxx
View File

@ -372,7 +372,6 @@ void TestPartitionedDataSet(vtkm::Id num, bool useGhost, FilterType fType)
seedArray = vtkm::cont::make_ArrayHandle({ vtkm::Particle(vtkm::Vec3f(.2f, 1.0f, .2f), 0),
vtkm::Particle(vtkm::Vec3f(.2f, 2.0f, .2f), 1) });
vtkm::Id numSeeds = seedArray.GetNumberOfValues();
if (fType == FilterType::STREAMLINE || fType == FilterType::PATHLINE)
{
vtkm::cont::PartitionedDataSet out;
@ -382,7 +381,6 @@ void TestPartitionedDataSet(vtkm::Id num, bool useGhost, FilterType fType)
streamline.SetStepSize(0.1f);
streamline.SetNumberOfSteps(100000);
streamline.SetSeeds(seedArray);
streamline.SetActiveField(fieldName);
out = streamline.Execute(pds);
}
@ -590,8 +588,6 @@ void TestStreamlineFilters()
FilterType::PATHLINE,
FilterType::PATH_PARTICLE };
//fTypes = {FilterType::PARTICLE_ADVECTION,FilterType::STREAMLINE};
//fTypes = {FilterType::PATHLINE};
for (int n = 1; n < 3; n++)
{
for (auto useGhost : flags)

10
vtkm/filter/flow/testing/UnitTestStreamlineFilterWarpX.cxx
View File

@ -11,7 +11,7 @@
#include <vtkm/Particle.h>
#include <vtkm/cont/testing/Testing.h>
#include <vtkm/filter/flow/Streamline.h>
#include <vtkm/filter/flow/WarpXStreamline.h>
#include <vtkm/io/VTKDataSetReader.h>
namespace
@ -43,7 +43,7 @@ void GenerateChargedParticles(const vtkm::cont::ArrayHandle<vtkm::Vec3f>& pos,
}
}
void TestStreamlineFilters()
void TestFilters()
{
std::string particleFile = vtkm::cont::testing::Testing::DataPath("misc/warpXparticles.vtk");
std::string fieldFile = vtkm::cont::testing::Testing::DataPath("misc/warpXfields.vtk");
@ -88,12 +88,10 @@ void TestStreamlineFilters()
1.0 / (SPEED_OF_LIGHT * vtkm::Sqrt(1. / spacing[0] + 1. / spacing[1] + 1. / spacing[2])));
std::cout << "CFL length : " << length << std::endl;
vtkm::filter::flow::Streamline streamline;
vtkm::filter::flow::WarpXStreamline streamline;
streamline.SetStepSize(length);
streamline.SetNumberOfSteps(steps);
streamline.SetSeeds(seeds);
streamline.SetVectorFieldType(vtkm::filter::flow::VectorFieldType::ELECTRO_MAGNETIC_FIELD_TYPE);
streamline.SetEField("E");
streamline.SetBField("B");
@ -109,5 +107,5 @@ void TestStreamlineFilters()
int UnitTestStreamlineFilterWarpX(int argc, char* argv[])
{
return vtkm::cont::testing::Testing::Run(TestStreamlineFilters, argc, argv);
return vtkm::cont::testing::Testing::Run(TestFilters, argc, argv);
}

121
vtkm/filter/flow/testing/UnitTestWorkletParticleAdvection.cxx
View File

@ -14,6 +14,7 @@
#include <vtkm/cont/ArrayHandle.h>
#include <vtkm/cont/DataSet.h>
#include <vtkm/cont/testing/Testing.h>
#include <vtkm/filter/flow/worklet/Analysis.h>
#include <vtkm/filter/flow/worklet/EulerIntegrator.h>
#include <vtkm/filter/flow/worklet/Field.h>
#include <vtkm/filter/flow/worklet/GridEvaluators.h>
@ -21,6 +22,7 @@
#include <vtkm/filter/flow/worklet/Particles.h>
#include <vtkm/filter/flow/worklet/RK4Integrator.h>
#include <vtkm/filter/flow/worklet/Stepper.h>
#include <vtkm/filter/flow/worklet/Termination.h>
#include <vtkm/filter/mesh_info/GhostCellClassify.h>
#include <vtkm/io/VTKDataSetReader.h>
#include <vtkm/worklet/testing/GenerateTestDataSets.h>
@ -369,6 +371,8 @@ void TestGhostCellEvaluators()
using GridEvalType = vtkm::worklet::flow::GridEvaluator<FieldType>;
using RK4Type = vtkm::worklet::flow::RK4Integrator<GridEvalType>;
using Stepper = vtkm::worklet::flow::Stepper<RK4Type, GridEvalType>;
using Termination = vtkm::worklet::flow::NormalTermination;
using Analysis = vtkm::worklet::flow::NoAnalysis<vtkm::Particle>;
constexpr vtkm::Id nX = 6;
constexpr vtkm::Id nY = 6;
@ -407,9 +411,11 @@ void TestGhostCellEvaluators()
seeds.push_back(vtkm::Particle(vtkm::Vec3f(3, 3, 3), 3));
auto seedArray = vtkm::cont::make_ArrayHandle(seeds, vtkm::CopyFlag::Off);
auto res = pa.Run(rk4, seedArray, 10000);
Termination termination(10000);
Analysis analysis;
pa.Run(rk4, seedArray, termination, analysis);
auto posPortal = res.Particles.ReadPortal();
auto posPortal = analysis.Particles.ReadPortal();
vtkm::Id numSeeds = seedArray.GetNumberOfValues();
for (vtkm::Id i = 0; i < numSeeds; i++)
{
@ -428,8 +434,7 @@ void TestGhostCellEvaluators()
}
}
void ValidateParticleAdvectionResult(
const vtkm::worklet::flow::ParticleAdvectionResult<vtkm::Particle>& res,
void ValidateParticleAdvectionResult(const vtkm::worklet::flow::NoAnalysis<vtkm::Particle>& res,
vtkm::Id nSeeds,
vtkm::Id maxSteps)
{
@ -449,7 +454,7 @@ void ValidateParticleAdvectionResult(
}
}
void ValidateStreamlineResult(const vtkm::worklet::flow::StreamlineResult<vtkm::Particle>& res,
void ValidateStreamlineResult(const vtkm::worklet::flow::StreamlineAnalysis<vtkm::Particle>& res,
vtkm::Id nSeeds,
vtkm::Id maxSteps)
{
@ -471,6 +476,9 @@ void TestIntegrators()
using FieldHandle = vtkm::cont::ArrayHandle<vtkm::Vec3f>;
using FieldType = vtkm::worklet::flow::VelocityField<FieldHandle>;
using GridEvalType = vtkm::worklet::flow::GridEvaluator<FieldType>;
using Termination = vtkm::worklet::flow::NormalTermination;
using Analysis = vtkm::worklet::flow::NoAnalysis<vtkm::Particle>;
const vtkm::Id3 dims(5, 5, 5);
const vtkm::Bounds bounds(0., 1., 0., 1., .0, .1);
@ -496,22 +504,25 @@ void TestIntegrators()
GenerateRandomParticles(points, 3, bounds);
vtkm::worklet::flow::ParticleAdvection pa;
vtkm::worklet::flow::ParticleAdvectionResult<vtkm::Particle> res;
Termination termination(maxSteps);
Analysis analysis;
{
auto seeds = vtkm::cont::make_ArrayHandle(points, vtkm::CopyFlag::On);
using IntegratorType = vtkm::worklet::flow::RK4Integrator<GridEvalType>;
using Stepper = vtkm::worklet::flow::Stepper<IntegratorType, GridEvalType>;
Stepper rk4(eval, stepSize);
res = pa.Run(rk4, seeds, maxSteps);
ValidateParticleAdvectionResult(res, nSeeds, maxSteps);
pa.Run(rk4, seeds, termination, analysis);
ValidateParticleAdvectionResult(analysis, nSeeds, maxSteps);
}
{
auto seeds = vtkm::cont::make_ArrayHandle(points, vtkm::CopyFlag::On);
using IntegratorType = vtkm::worklet::flow::EulerIntegrator<GridEvalType>;
using Stepper = vtkm::worklet::flow::Stepper<IntegratorType, GridEvalType>;
Stepper euler(eval, stepSize);
res = pa.Run(euler, seeds, maxSteps);
ValidateParticleAdvectionResult(res, nSeeds, maxSteps);
pa.Run(euler, seeds, termination, analysis);
ValidateParticleAdvectionResult(analysis, nSeeds, maxSteps);
//res = pa.Run(euler, seeds, maxSteps);
//ValidateParticleAdvectionResult(res, nSeeds, maxSteps);
}
}
}
@ -523,6 +534,9 @@ void TestParticleWorkletsWithDataSetTypes()
using GridEvalType = vtkm::worklet::flow::GridEvaluator<FieldType>;
using RK4Type = vtkm::worklet::flow::RK4Integrator<GridEvalType>;
using Stepper = vtkm::worklet::flow::Stepper<RK4Type, GridEvalType>;
using Termination = vtkm::worklet::flow::NormalTermination;
using PAnalysis = vtkm::worklet::flow::NoAnalysis<vtkm::Particle>;
using SAnalysis = vtkm::worklet::flow::StreamlineAnalysis<vtkm::Particle>;
vtkm::FloatDefault stepSize = 0.01f;
const vtkm::Id3 dims(5, 5, 5);
@ -574,34 +588,36 @@ void TestParticleWorkletsWithDataSetTypes()
if (i < 2)
{
vtkm::worklet::flow::ParticleAdvection pa;
vtkm::worklet::flow::ParticleAdvectionResult<vtkm::Particle> res;
Termination termination(maxSteps);
PAnalysis analysis;
if (i == 0)
{
auto seeds = vtkm::cont::make_ArrayHandle(pts, vtkm::CopyFlag::On);
res = pa.Run(rk4, seeds, maxSteps);
pa.Run(rk4, seeds, termination, analysis);
}
else
{
auto seeds = vtkm::cont::make_ArrayHandle(pts2, vtkm::CopyFlag::On);
res = pa.Run(rk4, seeds, maxSteps);
pa.Run(rk4, seeds, termination, analysis);
}
ValidateParticleAdvectionResult(res, nSeeds, maxSteps);
ValidateParticleAdvectionResult(analysis, nSeeds, maxSteps);
}
else
{
vtkm::worklet::flow::Streamline s;
vtkm::worklet::flow::StreamlineResult<vtkm::Particle> res;
vtkm::worklet::flow::ParticleAdvection pa;
Termination termination(maxSteps);
SAnalysis analysis(maxSteps);
if (i == 2)
{
auto seeds = vtkm::cont::make_ArrayHandle(pts, vtkm::CopyFlag::On);
res = s.Run(rk4, seeds, maxSteps);
pa.Run(rk4, seeds, termination, analysis);
}
else
{
auto seeds = vtkm::cont::make_ArrayHandle(pts2, vtkm::CopyFlag::On);
res = s.Run(rk4, seeds, maxSteps);
pa.Run(rk4, seeds, termination, analysis);
}
ValidateStreamlineResult(res, nSeeds, maxSteps);
ValidateStreamlineResult(analysis, nSeeds, maxSteps);
}
}
}
@ -629,6 +645,8 @@ void TestParticleStatus()
using GridEvalType = vtkm::worklet::flow::GridEvaluator<FieldType>;
using RK4Type = vtkm::worklet::flow::RK4Integrator<GridEvalType>;
using Stepper = vtkm::worklet::flow::Stepper<RK4Type, GridEvalType>;
using Termination = vtkm::worklet::flow::NormalTermination;
using Analysis = vtkm::worklet::flow::NoAnalysis<vtkm::Particle>;
vtkm::Id maxSteps = 1000;
vtkm::FloatDefault stepSize = 0.01f;
@ -644,7 +662,10 @@ void TestParticleStatus()
pts.push_back(vtkm::Particle(vtkm::Vec3f(-1, -1, -1), 1));
auto seedsArray = vtkm::cont::make_ArrayHandle(pts, vtkm::CopyFlag::On);
pa.Run(rk4, seedsArray, maxSteps);
Termination termination(maxSteps);
Analysis analysis;
pa.Run(rk4, seedsArray, termination, analysis);
auto portal = seedsArray.ReadPortal();
bool tookStep0 = portal.Get(0).GetStatus().CheckTookAnySteps();
@ -679,6 +700,9 @@ void TestWorkletsBasic()
using GridEvalType = vtkm::worklet::flow::GridEvaluator<FieldType>;
using RK4Type = vtkm::worklet::flow::RK4Integrator<GridEvalType>;
using Stepper = vtkm::worklet::flow::Stepper<RK4Type, GridEvalType>;
using Termination = vtkm::worklet::flow::NormalTermination;
using PAnalysis = vtkm::worklet::flow::NoAnalysis<vtkm::Particle>;
using SAnalysis = vtkm::worklet::flow::StreamlineAnalysis<vtkm::Particle>;
vtkm::FloatDefault stepSize = 0.01f;
const vtkm::Id3 dims(5, 5, 5);
@ -738,15 +762,15 @@ void TestWorkletsBasic()
if (w == "particleAdvection")
{
vtkm::worklet::flow::ParticleAdvection pa;
vtkm::worklet::flow::ParticleAdvectionResult<vtkm::Particle> res;
res = pa.Run(rk4, seedsArray, maxSteps);
Termination termination(maxSteps);
PAnalysis analysis;
pa.Run(rk4, seedsArray, termination, analysis);
vtkm::Id numRequiredPoints = static_cast<vtkm::Id>(endPts.size());
VTKM_TEST_ASSERT(res.Particles.GetNumberOfValues() == numRequiredPoints,
VTKM_TEST_ASSERT(analysis.Particles.GetNumberOfValues() == numRequiredPoints,
"Wrong number of points in particle advection result.");
auto portal = res.Particles.ReadPortal();
for (vtkm::Id i = 0; i < res.Particles.GetNumberOfValues(); i++)
auto portal = analysis.Particles.ReadPortal();
for (vtkm::Id i = 0; i < analysis.Particles.GetNumberOfValues(); i++)
{
VTKM_TEST_ASSERT(portal.Get(i).GetPosition() == endPts[static_cast<std::size_t>(i)],
"Particle advection point is wrong");
@ -762,18 +786,19 @@ void TestWorkletsBasic()
}
else if (w == "streamline")
{
vtkm::worklet::flow::Streamline s;
vtkm::worklet::flow::StreamlineResult<vtkm::Particle> res;
res = s.Run(rk4, seedsArray, maxSteps);
vtkm::worklet::flow::ParticleAdvection pa;
Termination termination(maxSteps);
SAnalysis analysis(maxSteps);
pa.Run(rk4, seedsArray, termination, analysis);
vtkm::Id numRequiredPoints = static_cast<vtkm::Id>(samplePts.size());
VTKM_TEST_ASSERT(res.Positions.GetNumberOfValues() == numRequiredPoints,
VTKM_TEST_ASSERT(analysis.Streams.GetNumberOfValues() == numRequiredPoints,
"Wrong number of points in streamline result.");
//Make sure all the points match.
auto parPortal = res.Particles.ReadPortal();
for (vtkm::Id i = 0; i < res.Particles.GetNumberOfValues(); i++)
auto parPortal = analysis.Particles.ReadPortal();
for (vtkm::Id i = 0; i < analysis.Particles.GetNumberOfValues(); i++)
{
VTKM_TEST_ASSERT(parPortal.Get(i).GetPosition() == endPts[static_cast<std::size_t>(i)],
"Streamline end point is wrong");
@ -786,19 +811,19 @@ void TestWorkletsBasic()
"Streamline particle did not terminate");
}
auto posPortal = res.Positions.ReadPortal();
for (vtkm::Id i = 0; i < res.Positions.GetNumberOfValues(); i++)
auto posPortal = analysis.Streams.ReadPortal();
for (vtkm::Id i = 0; i < analysis.Streams.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::Id numCells = analysis.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,
VTKM_TEST_ASSERT(analysis.PolyLines.GetCellShape(i) == vtkm::CELL_SHAPE_POLY_LINE,
"Wrong cell type in streamline.");
VTKM_TEST_ASSERT(res.PolyLines.GetNumberOfPointsInCell(i) ==
VTKM_TEST_ASSERT(analysis.PolyLines.GetNumberOfPointsInCell(i) ==
static_cast<vtkm::Id>(maxSteps + 1),
"Wrong number of points in streamline cell");
}
@ -869,6 +894,9 @@ void TestParticleAdvectionFile(const std::string& fileName,
using GridEvalType = vtkm::worklet::flow::GridEvaluator<FieldType>;
using RK4Type = vtkm::worklet::flow::RK4Integrator<GridEvalType>;
using Stepper = vtkm::worklet::flow::Stepper<RK4Type, GridEvalType>;
using Termination = vtkm::worklet::flow::NormalTermination;
using PAnalysis = vtkm::worklet::flow::NoAnalysis<vtkm::Particle>;
using SAnalysis = vtkm::worklet::flow::StreamlineAnalysis<vtkm::Particle>;
VTKM_TEST_ASSERT(ds.HasField(fieldName), "Data set missing a field named ", fieldName);
vtkm::cont::Field& field = ds.GetField(fieldName);
@ -887,6 +915,7 @@ void TestParticleAdvectionFile(const std::string& fileName,
FieldType velocities(fieldArray);
GridEvalType eval(ds.GetCoordinateSystem(), ds.GetCellSet(), velocities);
Stepper rk4(eval, stepSize);
Termination termination(maxSteps);
for (int i = 0; i < 2; i++)
{
@ -898,18 +927,16 @@ void TestParticleAdvectionFile(const std::string& fileName,
if (i == 0)
{
vtkm::worklet::flow::ParticleAdvection pa;
vtkm::worklet::flow::ParticleAdvectionResult<vtkm::Particle> res;
res = pa.Run(rk4, seedArray, maxSteps);
ValidateResult(res, maxSteps, endPts);
PAnalysis analysis;
pa.Run(rk4, seedArray, termination, analysis);
ValidateResult(analysis, maxSteps, endPts);
}
else if (i == 1)
{
vtkm::worklet::flow::Streamline s;
vtkm::worklet::flow::StreamlineResult<vtkm::Particle> res;
res = s.Run(rk4, seedArray, maxSteps);
ValidateResult(res, maxSteps, endPts);
vtkm::worklet::flow::ParticleAdvection pa;
SAnalysis analysis(maxSteps);
pa.Run(rk4, seedArray, termination, analysis);
ValidateResult(analysis, maxSteps, endPts);
}
}
}

1
vtkm/filter/flow/vtkm.module
View File

@ -11,6 +11,7 @@ OPTIONAL_DEPENDS
TEST_DEPENDS
vtkm_filter_geometry_refinement
vtkm_filter_mesh_info
vtkm_filter_flow
vtkm_io
TEST_OPTIONAL_DEPENDS
vtkm_rendering_testing

231
vtkm/filter/flow/worklet/Analysis.cxx Normal file
View File

@ -0,0 +1,231 @@
//=============================================================================
//
// 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.
//
//=============================================================================
#define vtk_m_filter_flow_worklet_Analysis_cxx
#include <vtkm/filter/flow/worklet/Analysis.h>
#include <vtkm/cont/Algorithm.h>
#include <vtkm/cont/ArrayCopy.h>
#include <vtkm/cont/ConvertNumComponentsToOffsets.h>
#include <vtkm/cont/DataSet.h>
#include <vtkm/cont/Invoker.h>
#include <vtkm/cont/ParticleArrayCopy.h>
#include <vtkm/worklet/WorkletMapField.h>
namespace vtkm
{
namespace worklet
{
namespace flow
{
template <typename ParticleType>
VTKM_CONT bool NoAnalysis<ParticleType>::MakeDataSet(
vtkm::cont::DataSet& dataset,
const std::vector<NoAnalysis<ParticleType>>& results)
{
size_t nResults = results.size();
std::vector<vtkm::cont::ArrayHandle<ParticleType>> allParticles;
allParticles.reserve(nResults);
for (const auto& vres : results)
allParticles.emplace_back(vres.Particles);
vtkm::cont::ArrayHandle<vtkm::Vec3f> pts;
vtkm::cont::ParticleArrayCopy(allParticles, pts);
vtkm::Id numPoints = pts.GetNumberOfValues();
if (numPoints > 0)
{
//Create coordinate system and vertex cell set.
dataset.AddCoordinateSystem(vtkm::cont::CoordinateSystem("coordinates", pts));
vtkm::cont::CellSetSingleType<> cells;
vtkm::cont::ArrayHandleIndex conn(numPoints);
vtkm::cont::ArrayHandle<vtkm::Id> connectivity;
vtkm::cont::ArrayCopy(conn, connectivity);
cells.Fill(numPoints, vtkm::CELL_SHAPE_VERTEX, 1, connectivity);
dataset.SetCellSet(cells);
}
return true;
}
namespace detail
{
class GetSteps : public vtkm::worklet::WorkletMapField
{
public:
VTKM_CONT
GetSteps() {}
using ControlSignature = void(FieldIn, FieldOut);
using ExecutionSignature = void(_1, _2);
template <typename ParticleType>
VTKM_EXEC void operator()(const ParticleType& p, vtkm::Id& numSteps) const
{
numSteps = p.GetNumberOfSteps();
}
};
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)
template <typename ParticleType>
VTKM_EXEC void operator()(const ParticleType& p,
const vtkm::Id& initialNumSteps,
vtkm::Id& diff) const
{
diff = 1 + p.GetNumberOfSteps() - initialNumSteps;
}
};
} // namespace detail
template <typename ParticleType>
VTKM_CONT void StreamlineAnalysis<ParticleType>::InitializeAnalysis(
const vtkm::cont::ArrayHandle<ParticleType>& particles)
{
this->NumParticles = particles.GetNumberOfValues();
//Create ValidPointArray initialized to zero.
vtkm::cont::ArrayHandleConstant<vtkm::Id> validity(0, this->NumParticles * (this->MaxSteps + 1));
vtkm::cont::ArrayCopy(validity, this->Validity);
//Create StepCountArray initialized to zero.
vtkm::cont::ArrayHandleConstant<vtkm::Id> streamLengths(0, this->NumParticles);
vtkm::cont::ArrayCopy(streamLengths, this->StreamLengths);
// Initialize InitLengths
vtkm::Id numSeeds = static_cast<vtkm::Id>(particles.GetNumberOfValues());
vtkm::cont::ArrayHandleIndex idxArray(numSeeds);
vtkm::cont::Invoker invoker;
invoker(detail::GetSteps{}, particles, this->InitialLengths);
}
template <typename ParticleType>
VTKM_CONT void StreamlineAnalysis<ParticleType>::FinalizeAnalysis(
vtkm::cont::ArrayHandle<ParticleType>& particles)
{
vtkm::Id numSeeds = particles.GetNumberOfValues();
vtkm::cont::ArrayHandle<vtkm::Vec3f> positions;
vtkm::cont::Algorithm::CopyIf(this->Streams, this->Validity, positions, IsOne());
vtkm::cont::Algorithm::Copy(positions, this->Streams);
// Create the cells
vtkm::cont::ArrayHandle<vtkm::Id> numPoints;
vtkm::cont::Invoker invoker;
invoker(detail::ComputeNumPoints{}, particles, this->InitialLengths, numPoints);
vtkm::cont::ArrayHandle<vtkm::Id> cellIndex;
vtkm::Id connectivityLen = vtkm::cont::Algorithm::ScanExclusive(numPoints, cellIndex);
vtkm::cont::ArrayHandleIndex connCount(connectivityLen);
vtkm::cont::ArrayHandle<vtkm::Id> connectivity;
vtkm::cont::ArrayCopy(connCount, connectivity);
vtkm::cont::ArrayHandle<vtkm::UInt8> cellTypes;
auto polyLineShape =
vtkm::cont::make_ArrayHandleConstant<vtkm::UInt8>(vtkm::CELL_SHAPE_POLY_LINE, numSeeds);
vtkm::cont::ArrayCopy(polyLineShape, cellTypes);
auto offsets = vtkm::cont::ConvertNumComponentsToOffsets(numPoints);
this->PolyLines.Fill(this->Streams.GetNumberOfValues(), cellTypes, connectivity, offsets);
this->Particles = particles;
}
template <typename ParticleType>
VTKM_CONT bool StreamlineAnalysis<ParticleType>::MakeDataSet(
vtkm::cont::DataSet& dataset,
const std::vector<StreamlineAnalysis<ParticleType>>& results)
{
size_t nResults = results.size();
if (nResults == 1)
{
const auto& res = results[0];
dataset.AddCoordinateSystem(vtkm::cont::CoordinateSystem("coordinates", res.Streams));
dataset.SetCellSet(res.PolyLines);
}
else
{
std::vector<vtkm::Id> posOffsets(nResults, 0);
vtkm::Id totalNumCells = 0, totalNumPts = 0;
for (std::size_t i = 0; i < nResults; i++)
{
const auto& res = results[i];
if (i == 0)
posOffsets[i] = 0;
else
posOffsets[i] = totalNumPts;
totalNumPts += res.Streams.GetNumberOfValues();
totalNumCells += res.PolyLines.GetNumberOfCells();
}
//Append all the points together.
vtkm::cont::ArrayHandle<vtkm::Vec3f> appendPts;
appendPts.Allocate(totalNumPts);
for (std::size_t i = 0; i < nResults; i++)
{
const auto& res = results[i];
// copy all values into appendPts starting at offset.
vtkm::cont::Algorithm::CopySubRange(
res.Streams, 0, res.Streams.GetNumberOfValues(), appendPts, posOffsets[i]);
}
dataset.AddCoordinateSystem(vtkm::cont::CoordinateSystem("coordinates", appendPts));
//Create polylines.
std::vector<vtkm::Id> numPtsPerCell(static_cast<std::size_t>(totalNumCells));
std::size_t off = 0;
for (std::size_t i = 0; i < nResults; i++)
{
const auto& res = results[i];
vtkm::Id nCells = res.PolyLines.GetNumberOfCells();
for (vtkm::Id j = 0; j < nCells; j++)
numPtsPerCell[off++] = static_cast<vtkm::Id>(res.PolyLines.GetNumberOfPointsInCell(j));
}
auto numPointsPerCellArray = vtkm::cont::make_ArrayHandle(numPtsPerCell, vtkm::CopyFlag::Off);
vtkm::cont::ArrayHandle<vtkm::Id> cellIndex;
vtkm::Id connectivityLen =
vtkm::cont::Algorithm::ScanExclusive(numPointsPerCellArray, cellIndex);
vtkm::cont::ArrayHandleIndex connCount(connectivityLen);
vtkm::cont::ArrayHandle<vtkm::Id> connectivity;
vtkm::cont::ArrayCopy(connCount, connectivity);
vtkm::cont::ArrayHandle<vtkm::UInt8> cellTypes;
auto polyLineShape =
vtkm::cont::make_ArrayHandleConstant<vtkm::UInt8>(vtkm::CELL_SHAPE_POLY_LINE, totalNumCells);
vtkm::cont::ArrayCopy(polyLineShape, cellTypes);
auto offsets = vtkm::cont::ConvertNumComponentsToOffsets(numPointsPerCellArray);
vtkm::cont::CellSetExplicit<> polyLines;
polyLines.Fill(totalNumPts, cellTypes, connectivity, offsets);
dataset.SetCellSet(polyLines);
}
return true;
}
template class VTKM_FILTER_FLOW_EXPORT NoAnalysis<vtkm::Particle>;
template class VTKM_FILTER_FLOW_EXPORT NoAnalysis<vtkm::ChargedParticle>;
template class VTKM_FILTER_FLOW_EXPORT StreamlineAnalysis<vtkm::Particle>;
template class VTKM_FILTER_FLOW_EXPORT StreamlineAnalysis<vtkm::ChargedParticle>;
} // namespace flow
} // namespace worklet
} // namespace vtkm

247
vtkm/filter/flow/worklet/Analysis.h Normal file
View File

@ -0,0 +1,247 @@
//=============================================================================
//
// 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 vtkm_worklet_particleadvection_analysis
#define vtkm_worklet_particleadvection_analysis
#include <vtkm/Particle.h>
#include <vtkm/Types.h>
#include <vtkm/cont/ArrayHandle.h>
#include <vtkm/cont/DataSet.h>
#include <vtkm/cont/ExecutionObjectBase.h>
#include <vtkm/filter/flow/vtkm_filter_flow_export.h>
namespace vtkm
{
namespace worklet
{
namespace flow
{
template <typename ParticleType>
class VTKM_FILTER_FLOW_EXPORT NoAnalysisExec
{
public:
VTKM_EXEC_CONT
NoAnalysisExec() {}
VTKM_EXEC void PreStepAnalyze(const vtkm::Id index, const ParticleType& particle)
{
(void)index;
(void)particle;
}
//template <typename ParticleType>
VTKM_EXEC void Analyze(const vtkm::Id index,
const ParticleType& oldParticle,
const ParticleType& newParticle)
{
// Do nothing
(void)index;
(void)oldParticle;
(void)newParticle;
}
};
template <typename ParticleType>
class NoAnalysis : public vtkm::cont::ExecutionObjectBase
{
public:
// Intended to store advected particles after Finalize
vtkm::cont::ArrayHandle<ParticleType> Particles;
VTKM_CONT
NoAnalysis()
: Particles()
{
}
VTKM_CONT
void UseAsTemplate(const NoAnalysis& other) { (void)other; }
VTKM_CONT
//template <typename ParticleType>
void InitializeAnalysis(const vtkm::cont::ArrayHandle<ParticleType>& particles)
{
(void)particles;
}
VTKM_CONT
//template <typename ParticleType>
void FinalizeAnalysis(vtkm::cont::ArrayHandle<ParticleType>& particles)
{
this->Particles = particles; //, vtkm::CopyFlag::Off);
}
VTKM_CONT NoAnalysisExec<ParticleType> PrepareForExecution(vtkm::cont::DeviceAdapterId device,
vtkm::cont::Token& token) const
{
(void)device;
(void)token;
return NoAnalysisExec<ParticleType>();
}
VTKM_CONT bool SupportPushOutOfBounds() const { return true; }
VTKM_CONT static bool MakeDataSet(vtkm::cont::DataSet& dataset,
const std::vector<NoAnalysis>& results);
};
template <typename ParticleType>
class VTKM_FILTER_FLOW_EXPORT StreamlineAnalysisExec
{
public:
VTKM_EXEC_CONT
StreamlineAnalysisExec()
: NumParticles(0)
, MaxSteps(0)
, Streams()
, StreamLengths()
, Validity()
{
}
VTKM_CONT
StreamlineAnalysisExec(vtkm::Id numParticles,
vtkm::Id maxSteps,
const vtkm::cont::ArrayHandle<vtkm::Vec3f>& streams,
const vtkm::cont::ArrayHandle<vtkm::Id>& streamLengths,
const vtkm::cont::ArrayHandle<vtkm::Id>& validity,
vtkm::cont::DeviceAdapterId device,
vtkm::cont::Token& token)
: NumParticles(numParticles)
, MaxSteps(maxSteps + 1)
{
Streams = streams.PrepareForOutput(this->NumParticles * this->MaxSteps, device, token);
StreamLengths = streamLengths.PrepareForInPlace(device, token);
Validity = validity.PrepareForInPlace(device, token);
}
VTKM_EXEC void PreStepAnalyze(const vtkm::Id index, const ParticleType& particle)
{
vtkm::Id streamLength = this->StreamLengths.Get(index);
if (streamLength == 0)
{
this->StreamLengths.Set(index, 1);
vtkm::Id loc = index * MaxSteps;
this->Streams.Set(loc, particle.GetPosition());
this->Validity.Set(loc, 1);
}
}
//template <typename ParticleType>
VTKM_EXEC void Analyze(const vtkm::Id index,
const ParticleType& oldParticle,
const ParticleType& newParticle)
{
(void)oldParticle;
vtkm::Id streamLength = this->StreamLengths.Get(index);
vtkm::Id loc = index * MaxSteps + streamLength;
this->StreamLengths.Set(index, ++streamLength);
this->Streams.Set(loc, newParticle.GetPosition());
this->Validity.Set(loc, 1);
}
private:
using IdPortal = typename vtkm::cont::ArrayHandle<vtkm::Id>::WritePortalType;
using VecPortal = typename vtkm::cont::ArrayHandle<vtkm::Vec3f>::WritePortalType;
vtkm::Id NumParticles;
vtkm::Id MaxSteps;
VecPortal Streams;
IdPortal StreamLengths;
IdPortal Validity;
};
template <typename ParticleType>
class StreamlineAnalysis : public vtkm::cont::ExecutionObjectBase
{
public:
// Intended to store advected particles after Finalize
vtkm::cont::ArrayHandle<ParticleType> Particles;
vtkm::cont::ArrayHandle<vtkm::Vec3f> Streams;
vtkm::cont::CellSetExplicit<> PolyLines;
//Helper functor for compacting history
struct IsOne
{
template <typename T>
VTKM_EXEC_CONT bool operator()(const T& x) const
{
return x == T(1);
}
};
VTKM_CONT
StreamlineAnalysis()
: Particles()
, MaxSteps(0)
{
}
VTKM_CONT
StreamlineAnalysis(vtkm::Id maxSteps)
: Particles()
, MaxSteps(maxSteps)
{
}
VTKM_CONT
void UseAsTemplate(const StreamlineAnalysis& other) { this->MaxSteps = other.MaxSteps; }
VTKM_CONT StreamlineAnalysisExec<ParticleType> PrepareForExecution(
vtkm::cont::DeviceAdapterId device,
vtkm::cont::Token& token) const
{
return StreamlineAnalysisExec<ParticleType>(this->NumParticles,
this->MaxSteps,
this->Streams,
this->StreamLengths,
this->Validity,
device,
token);
}
VTKM_CONT bool SupportPushOutOfBounds() const { return true; }
VTKM_CONT
void InitializeAnalysis(const vtkm::cont::ArrayHandle<ParticleType>& particles);
VTKM_CONT
//template <typename ParticleType>
void FinalizeAnalysis(vtkm::cont::ArrayHandle<ParticleType>& particles);
VTKM_CONT static bool MakeDataSet(vtkm::cont::DataSet& dataset,
const std::vector<StreamlineAnalysis>& results);
private:
vtkm::Id NumParticles;
vtkm::Id MaxSteps;
vtkm::cont::ArrayHandle<vtkm::Id> StreamLengths;
vtkm::cont::ArrayHandle<vtkm::Id> InitialLengths;
vtkm::cont::ArrayHandle<vtkm::Id> Validity;
};
#ifndef vtk_m_filter_flow_worklet_Analysis_cxx
extern template class VTKM_FILTER_FLOW_TEMPLATE_EXPORT NoAnalysis<vtkm::Particle>;
extern template class VTKM_FILTER_FLOW_TEMPLATE_EXPORT NoAnalysis<vtkm::ChargedParticle>;
extern template class VTKM_FILTER_FLOW_TEMPLATE_EXPORT StreamlineAnalysis<vtkm::Particle>;
extern template class VTKM_FILTER_FLOW_TEMPLATE_EXPORT StreamlineAnalysis<vtkm::ChargedParticle>;
#endif //!vtk_m_filter_flow_worklet_Analysis_cxx
} // namespace flow
} // namespace worklet
} // namespace vtkm
#endif //vtkm_worklet_particleadvection_analysis

5
vtkm/filter/flow/worklet/CMakeLists.txt
View File

@ -9,7 +9,7 @@
##============================================================================
set(headers
ParticleAdvection.h
Analysis.h
CellInterpolationHelper.h
EulerIntegrator.h
Field.h
@ -18,11 +18,14 @@ set(headers
IntegratorStatus.h
LagrangianStructures.h
Particles.h
ParticleAdvection.h
ParticleAdvectionWorklets.h
RK4Integrator.h
TemporalGridEvaluators.h
Stepper.h
StreamSurface.h
TemporalGridEvaluators.h
Termination.h
)
#-----------------------------------------------------------------------------

2
vtkm/filter/flow/worklet/EulerIntegrator.h
View File

@ -34,7 +34,7 @@ public:
}
template <typename Particle>
VTKM_EXEC IntegratorStatus CheckStep(Particle& particle,
VTKM_EXEC IntegratorStatus CheckStep(const Particle& particle,
vtkm::FloatDefault stepLength,
vtkm::Vec3f& velocity) const
{

100
vtkm/filter/flow/worklet/Field.h
View File

@ -31,6 +31,7 @@ class ExecutionVelocityField
public:
using FieldPortalType = typename FieldArrayType::ReadPortalType;
using Association = vtkm::cont::Field::Association;
using DelegateToField = std::false_type;
VTKM_CONT
ExecutionVelocityField(const FieldArrayType& velocityValues,
@ -68,17 +69,67 @@ public:
value = vtkm::make_Vec(velocityInterp);
}
template <typename Point, typename Locator, typename Helper>
VTKM_EXEC bool GetValue(const Point& vtkmNotUsed(point),
const vtkm::FloatDefault& vtkmNotUsed(time),
vtkm::VecVariable<Point, 2>& vtkmNotUsed(out),
const Locator& vtkmNotUsed(locator),
const Helper& vtkmNotUsed(helper)) const
{
//TODO Raise Error : Velocity Field should not allow this path
return false;
}
private:
FieldPortalType VelocityValues;
Association Assoc;
};
template <typename FieldArrayType>
class VelocityField : public vtkm::cont::ExecutionObjectBase
{
public:
using ExecutionType = ExecutionVelocityField<FieldArrayType>;
using Association = vtkm::cont::Field::Association;
VTKM_CONT
VelocityField() = default;
VTKM_CONT
VelocityField(const FieldArrayType& fieldValues)
: FieldValues(fieldValues)
, Assoc(vtkm::cont::Field::Association::Points)
{
}
VTKM_CONT
VelocityField(const FieldArrayType& fieldValues, const Association assoc)
: FieldValues(fieldValues)
, Assoc(assoc)
{
if (assoc != Association::Points && assoc != Association::Cells)
throw("Unsupported field association");
}
VTKM_CONT
const ExecutionType PrepareForExecution(vtkm::cont::DeviceAdapterId device,
vtkm::cont::Token& token) const
{
return ExecutionType(this->FieldValues, this->Assoc, device, token);
}
private:
FieldArrayType FieldValues;
Association Assoc;
};
template <typename FieldArrayType>
class ExecutionElectroMagneticField
{
public:
using FieldPortalType = typename FieldArrayType::ReadPortalType;
using Association = vtkm::cont::Field::Association;
using DelegateToField = std::false_type;
VTKM_CONT
ExecutionElectroMagneticField(const FieldArrayType& electricValues,
@ -124,50 +175,23 @@ public:
value = vtkm::make_Vec(electricInterp, magneticInterp);
}
template <typename Point, typename Locator, typename Helper>
VTKM_EXEC bool GetValue(const Point& vtkmNotUsed(point),
const vtkm::FloatDefault& vtkmNotUsed(time),
vtkm::VecVariable<Point, 2>& vtkmNotUsed(out),
const Locator& vtkmNotUsed(locator),
const Helper& vtkmNotUsed(helper)) const
{
//TODO : Raise Error : Velocity Field should not allow this path
return false;
}
private:
FieldPortalType ElectricValues;
FieldPortalType MagneticValues;
Association Assoc;
};
template <typename FieldArrayType>
class VelocityField : public vtkm::cont::ExecutionObjectBase
{
public:
using ExecutionType = ExecutionVelocityField<FieldArrayType>;
using Association = vtkm::cont::Field::Association;
VTKM_CONT
VelocityField() = default;
VTKM_CONT
VelocityField(const FieldArrayType& fieldValues)
: FieldValues(fieldValues)
, Assoc(vtkm::cont::Field::Association::Points)
{
}
VTKM_CONT
VelocityField(const FieldArrayType& fieldValues, const Association assoc)
: FieldValues(fieldValues)
, Assoc(assoc)
{
if (assoc != Association::Points && assoc != Association::Cells)
throw("Unsupported field association");
}
VTKM_CONT
const ExecutionType PrepareForExecution(vtkm::cont::DeviceAdapterId device,
vtkm::cont::Token& token) const
{
return ExecutionType(this->FieldValues, this->Assoc, device, token);
}
private:
FieldArrayType FieldValues;
Association Assoc;
};
template <typename FieldArrayType>
class ElectroMagneticField : public vtkm::cont::ExecutionObjectBase
{

37
vtkm/filter/flow/worklet/GridEvaluators.h
View File

@ -35,6 +35,7 @@ template <typename FieldType>
class ExecutionGridEvaluator
{
using GhostCellArrayType = vtkm::cont::ArrayHandle<vtkm::UInt8>;
using ExecFieldType = typename FieldType::ExecutionType;
public:
VTKM_CONT
@ -86,7 +87,7 @@ public:
}
template <typename Point>
VTKM_EXEC GridEvaluatorStatus Evaluate(const Point& point,
VTKM_EXEC GridEvaluatorStatus HelpEvaluate(const Point& point,
const vtkm::FloatDefault& time,
vtkm::VecVariable<Point, 2>& out) const
{
@ -134,10 +135,40 @@ public:
status.SetOk();
}
return status;
}
template <typename Point>
VTKM_EXEC GridEvaluatorStatus DeligateEvaluateToField(const Point& point,
const vtkm::FloatDefault& time,
vtkm::VecVariable<Point, 2>& out) const
{
GridEvaluatorStatus status;
status.SetOk();
// TODO: Allow for getting status from deligated work from Field
if (!this->Field.GetValue(point, time, out, this->Locator, this->InterpolationHelper))
{
status.SetFail();
status.SetSpatialBounds();
}
return status;
}
template <typename Point>
VTKM_EXEC GridEvaluatorStatus Evaluate(const Point& point,
const vtkm::FloatDefault& time,
vtkm::VecVariable<Point, 2>& out) const
{
if (!ExecFieldType::DelegateToField::value)
{
return this->HelpEvaluate(point, time, out);
}
else
{
return this->DeligateEvaluateToField(point, time, out);
}
}
private:
VTKM_EXEC bool InGhostCell(const vtkm::Id& cellId) const
{
@ -150,7 +181,7 @@ private:
using GhostCellPortal = typename vtkm::cont::ArrayHandle<vtkm::UInt8>::ReadPortalType;
vtkm::Bounds Bounds;
typename FieldType::ExecutionType Field;
ExecFieldType Field;
GhostCellPortal GhostCells;
bool HaveGhostCells;
vtkm::exec::CellInterpolationHelper InterpolationHelper;

113
vtkm/filter/flow/worklet/ParticleAdvection.h
View File

@ -48,58 +48,40 @@ public:
} //detail
template <typename ParticleType>
struct ParticleAdvectionResult
{
ParticleAdvectionResult()
: Particles()
{
}
ParticleAdvectionResult(const vtkm::cont::ArrayHandle<ParticleType>& p)
: Particles(p)
{
}
vtkm::cont::ArrayHandle<ParticleType> Particles;
};
class ParticleAdvection
{
public:
ParticleAdvection() {}
template <typename IntegratorType, typename ParticleType, typename ParticleStorage>
void Run2(const IntegratorType& it,
template <typename IntegratorType,
typename ParticleType,
typename ParticleStorage,
typename TerminationType,
typename AnalysisType>
void Run(const IntegratorType& it,
vtkm::cont::ArrayHandle<ParticleType, ParticleStorage>& particles,
vtkm::Id MaxSteps,
ParticleAdvectionResult<ParticleType>& result)
const TerminationType& termination,
AnalysisType& analysis)
{
vtkm::worklet::flow::ParticleAdvectionWorklet<IntegratorType, ParticleType> worklet;
worklet.Run(it, particles, MaxSteps);
result = ParticleAdvectionResult<ParticleType>(particles);
vtkm::worklet::flow::
ParticleAdvectionWorklet<IntegratorType, ParticleType, TerminationType, AnalysisType>
worklet;
worklet.Run(it, particles, termination, analysis);
}
template <typename IntegratorType, typename ParticleType, typename ParticleStorage>
ParticleAdvectionResult<ParticleType> Run(
const IntegratorType& it,
vtkm::cont::ArrayHandle<ParticleType, ParticleStorage>& particles,
vtkm::Id MaxSteps)
{
vtkm::worklet::flow::ParticleAdvectionWorklet<IntegratorType, ParticleType> worklet;
worklet.Run(it, particles, MaxSteps);
return ParticleAdvectionResult<ParticleType>(particles);
}
template <typename IntegratorType, typename ParticleType, typename PointStorage>
ParticleAdvectionResult<ParticleType> Run(
const IntegratorType& it,
template <typename IntegratorType,
typename ParticleType,
typename PointStorage,
typename TerminationType,
typename AnalysisType>
void Run(const IntegratorType& it,
const vtkm::cont::ArrayHandle<vtkm::Vec3f, PointStorage>& points,
vtkm::Id MaxSteps)
const TerminationType& termination,
AnalysisType& analysis)
{
vtkm::worklet::flow::ParticleAdvectionWorklet<IntegratorType, ParticleType> worklet;
vtkm::worklet::flow::
ParticleAdvectionWorklet<IntegratorType, ParticleType, TerminationType, AnalysisType>
worklet;
vtkm::cont::ArrayHandle<ParticleType> particles;
vtkm::cont::ArrayHandle<vtkm::Id> step, ids;
@ -117,54 +99,7 @@ public:
vtkm::cont::ArrayCopy(id, ids);
invoke(detail::CopyToParticle{}, points, ids, time, step, particles);
worklet.Run(it, particles, MaxSteps);
return ParticleAdvectionResult<ParticleType>(particles);
}
};
template <typename ParticleType>
struct StreamlineResult
{
StreamlineResult()
: Particles()
, Positions()
, PolyLines()
{
}
StreamlineResult(const vtkm::cont::ArrayHandle<ParticleType>& part,
const vtkm::cont::ArrayHandle<vtkm::Vec3f>& pos,
const vtkm::cont::CellSetExplicit<>& lines)
: Particles(part)
, Positions(pos)
, PolyLines(lines)
{
}
vtkm::cont::ArrayHandle<ParticleType> Particles;
vtkm::cont::ArrayHandle<vtkm::Vec3f> Positions;
vtkm::cont::CellSetExplicit<> PolyLines;
};
class Streamline
{
public:
Streamline() {}
template <typename IntegratorType, typename ParticleType, typename ParticleStorage>
StreamlineResult<ParticleType> Run(
const IntegratorType& it,
vtkm::cont::ArrayHandle<ParticleType, ParticleStorage>& particles,
vtkm::Id MaxSteps)
{
vtkm::worklet::flow::StreamlineWorklet<IntegratorType, ParticleType> worklet;
vtkm::cont::ArrayHandle<vtkm::Vec3f> positions;
vtkm::cont::CellSetExplicit<> polyLines;
worklet.Run(it, particles, MaxSteps, positions, polyLines);
return StreamlineResult<ParticleType>(particles, positions, polyLines);
worklet.Run(it, particles, termination, analysis);
}
};

163
vtkm/filter/flow/worklet/ParticleAdvectionWorklets.h
View File

@ -15,6 +15,7 @@
#include <vtkm/cont/CellSetExplicit.h>
#include <vtkm/cont/ConvertNumComponentsToOffsets.h>
#include <vtkm/cont/ExecutionObjectBase.h>
#include <vtkm/cont/Invoker.h>
#include <vtkm/Particle.h>
#include <vtkm/filter/flow/worklet/Particles.h>
@ -34,18 +35,26 @@ namespace flow
class ParticleAdvectWorklet : public vtkm::worklet::WorkletMapField
{
public:
using ControlSignature = void(FieldIn idx,
ExecObject integrator,
ExecObject integralCurve,
FieldIn maxSteps);
using ExecutionSignature = void(_1 idx, _2 integrator, _3 integralCurve, _4 maxSteps);
VTKM_EXEC_CONT
ParticleAdvectWorklet()
: PushOutOfBounds(true)
{
}
VTKM_EXEC_CONT
ParticleAdvectWorklet(bool pushOutOfBounds)
: PushOutOfBounds(pushOutOfBounds)
{
}
using ControlSignature = void(FieldIn idx, ExecObject integrator, ExecObject integralCurve);
using ExecutionSignature = void(_1 idx, _2 integrator, _3 integralCurve);
using InputDomain = _1;
template <typename IntegratorType, typename IntegralCurveType>
VTKM_EXEC void operator()(const vtkm::Id& idx,
const IntegratorType& integrator,
IntegralCurveType& integralCurve,
const vtkm::Id& maxSteps) const
IntegralCurveType& integralCurve) const
{
auto particle = integralCurve.GetParticle(idx);
vtkm::FloatDefault time = particle.GetTime();
@ -55,7 +64,7 @@ public:
// 1. you could have success AND at temporal boundary.
// 2. could you have success AND at spatial?
// 3. all three?
integralCurve.PreStepUpdate(idx);
integralCurve.PreStepUpdate(idx, particle);
do
{
particle = integralCurve.GetParticle(idx);
@ -69,7 +78,7 @@ public:
//We can't take a step inside spatial boundary.
//Try and take a step just past the boundary.
else if (status.CheckSpatialBounds())
else if (status.CheckSpatialBounds() && this->PushOutOfBounds)
{
status = integrator.SmallStep(particle, time, outpos);
if (status.CheckOk())
@ -78,16 +87,22 @@ public:
tookAnySteps = true;
}
}
integralCurve.StatusUpdate(idx, status, maxSteps);
integralCurve.StatusUpdate(idx, status);
} while (integralCurve.CanContinue(idx));
//Mark if any steps taken
integralCurve.UpdateTookSteps(idx, tookAnySteps);
}
private:
bool PushOutOfBounds;
};
template <typename IntegratorType, typename ParticleType>
template <typename IntegratorType,
typename ParticleType,
typename TerminationType,
typename AnalysisType>
class ParticleAdvectionWorklet
{
public:
@ -97,19 +112,17 @@ public:
void Run(const IntegratorType& integrator,
vtkm::cont::ArrayHandle<ParticleType>& particles,
vtkm::Id& MaxSteps)
const TerminationType& termination,
AnalysisType& analysis)
{
using ParticleAdvectWorkletType = vtkm::worklet::flow::ParticleAdvectWorklet;
using ParticleWorkletDispatchType =
typename vtkm::worklet::DispatcherMapField<ParticleAdvectWorkletType>;
using ParticleArrayType = vtkm::worklet::flow::Particles<ParticleType>;
using ParticleArrayType =
vtkm::worklet::flow::Particles<ParticleType, TerminationType, AnalysisType>;
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::ArrayHandleConstant<vtkm::Id> maxSteps(MaxSteps, numSeeds);
vtkm::cont::ArrayHandleIndex idxArray(numSeeds);
// TODO: The particle advection sometimes behaves incorrectly on CUDA if the stack size
// is not changed thusly. This is concerning as the compiler should be able to determine
// staticly the required stack depth. What is even more concerning is that the runtime
@ -121,114 +134,20 @@ public:
(void)stack;
#endif // VTKM_CUDA
ParticleArrayType particlesObj(particles, MaxSteps);
// Initialize all the pre-requisites needed to start analysis
// It's based on the existing properties of the particles,
// for e.g. the number of steps they've already taken
analysis.InitializeAnalysis(particles);
//Invoke particle advection worklet
ParticleWorkletDispatchType particleWorkletDispatch;
ParticleArrayType particlesObj(particles, termination, analysis);
particleWorkletDispatch.Invoke(idxArray, integrator, particlesObj, maxSteps);
}
};
vtkm::worklet::flow::ParticleAdvectWorklet worklet(analysis.SupportPushOutOfBounds());
namespace detail
{
class GetSteps : public vtkm::worklet::WorkletMapField
{
public:
VTKM_CONT
GetSteps() {}
using ControlSignature = void(FieldIn, FieldOut);
using ExecutionSignature = void(_1, _2);
vtkm::cont::Invoker invoker;
invoker(worklet, idxArray, integrator, particlesObj);
template <typename ParticleType>
VTKM_EXEC void operator()(const ParticleType& p, vtkm::Id& numSteps) const
{
numSteps = p.GetNumberOfSteps();
}
};
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.GetNumberOfSteps() - initalNumSteps)
template <typename ParticleType>
VTKM_EXEC void operator()(const ParticleType& p,
const vtkm::Id& initialNumSteps,
vtkm::Id& diff) const
{
diff = 1 + p.GetNumberOfSteps() - initialNumSteps;
}
};
} // namespace detail
template <typename IntegratorType, typename ParticleType>
class StreamlineWorklet
{
public:
template <typename PointStorage, typename PointStorage2>
void Run(const IntegratorType& it,
vtkm::cont::ArrayHandle<ParticleType, PointStorage>& particles,
vtkm::Id& MaxSteps,
vtkm::cont::ArrayHandle<vtkm::Vec3f, PointStorage2>& positions,
vtkm::cont::CellSetExplicit<>& polyLines)
{
using ParticleWorkletDispatchType =
typename vtkm::worklet::DispatcherMapField<vtkm::worklet::flow::ParticleAdvectWorklet>;
using StreamlineArrayType = vtkm::worklet::flow::StateRecordingParticles<ParticleType>;
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);
// This method uses the same workklet as ParticleAdvectionWorklet::Run (and more). Yet for
// some reason ParticleAdvectionWorklet::Run needs this adjustment while this method does
// not.
#ifdef VTKM_CUDA
// // This worklet needs some extra space on CUDA.
// vtkm::cont::cuda::internal::ScopedCudaStackSize stack(16 * 1024);
// (void)stack;
#endif // VTKM_CUDA
//Run streamline worklet
StreamlineArrayType streamlines(particles, MaxSteps);
ParticleWorkletDispatchType particleWorkletDispatch;
vtkm::cont::ArrayHandleConstant<vtkm::Id> maxSteps(MaxSteps, numSeeds);
particleWorkletDispatch.Invoke(idxArray, it, streamlines, maxSteps);
//Get the positions
streamlines.GetCompactedHistory(positions);
//Create the cells
vtkm::cont::ArrayHandle<vtkm::Id> numPoints;
vtkm::worklet::DispatcherMapField<detail::ComputeNumPoints> computeNumPointsDispatcher{ (
detail::ComputeNumPoints{}) };
computeNumPointsDispatcher.Invoke(particles, initialStepsTaken, numPoints);
vtkm::cont::ArrayHandle<vtkm::Id> cellIndex;
vtkm::Id connectivityLen = vtkm::cont::Algorithm::ScanExclusive(numPoints, cellIndex);
vtkm::cont::ArrayHandleIndex connCount(connectivityLen);
vtkm::cont::ArrayHandle<vtkm::Id> connectivity;
vtkm::cont::ArrayCopy(connCount, connectivity);
vtkm::cont::ArrayHandle<vtkm::UInt8> cellTypes;
auto polyLineShape =
vtkm::cont::make_ArrayHandleConstant<vtkm::UInt8>(vtkm::CELL_SHAPE_POLY_LINE, numSeeds);
vtkm::cont::ArrayCopy(polyLineShape, cellTypes);
auto offsets = vtkm::cont::ConvertNumComponentsToOffsets(numPoints);
polyLines.Fill(positions.GetNumberOfValues(), cellTypes, connectivity, offsets);
// Finalize the analysis and clear intermittant arrays.
analysis.FinalizeAnalysis(particles);
}
};

215
vtkm/filter/flow/worklet/Particles.h
View File

@ -14,6 +14,9 @@
#include <vtkm/cont/ExecutionObjectBase.h>
#include <vtkm/filter/flow/worklet/IntegratorStatus.h>
#include <vtkm/filter/flow/worklet/Analysis.h>
#include <vtkm/filter/flow/worklet/Termination.h>
namespace vtkm
{
namespace worklet
@ -21,31 +24,37 @@ namespace worklet
namespace flow
{
template <typename ParticleType>
template <typename ParticleType, typename TerminationType, typename AnalysisType>
class ParticleExecutionObject
{
public:
VTKM_EXEC_CONT
ParticleExecutionObject()
: Particles()
, MaxSteps(0)
, Termination()
, Analysis()
{
}
ParticleExecutionObject(vtkm::cont::ArrayHandle<ParticleType> particleArray,
vtkm::Id maxSteps,
const TerminationType& termination,
const AnalysisType& analysis,
vtkm::cont::DeviceAdapterId device,
vtkm::cont::Token& token)
: Termination(termination)
, Analysis(analysis)
{
Particles = particleArray.PrepareForInPlace(device, token);
MaxSteps = maxSteps;
}
VTKM_EXEC
ParticleType GetParticle(const vtkm::Id& idx) { return this->Particles.Get(idx); }
VTKM_EXEC
void PreStepUpdate(const vtkm::Id& vtkmNotUsed(idx)) {}
void PreStepUpdate(const vtkm::Id& idx, const ParticleType& particle)
{
this->Analysis.PreStepAnalyze(idx, particle);
}
VTKM_EXEC
void StepUpdate(const vtkm::Id& idx,
@ -57,19 +66,15 @@ public:
newParticle.SetPosition(pt);
newParticle.SetTime(time);
newParticle.SetNumberOfSteps(particle.GetNumberOfSteps() + 1);
this->Analysis.Analyze(idx, particle, newParticle);
this->Particles.Set(idx, newParticle);
}
VTKM_EXEC
void StatusUpdate(const vtkm::Id& idx,
const vtkm::worklet::flow::IntegratorStatus& status,
vtkm::Id maxSteps)
void StatusUpdate(const vtkm::Id& idx, const vtkm::worklet::flow::IntegratorStatus& status)
{
ParticleType p(this->GetParticle(idx));
if (p.GetNumberOfSteps() == maxSteps)
p.GetStatus().SetTerminate();
if (status.CheckFail())
p.GetStatus().SetFail();
if (status.CheckSpatialBounds())
@ -88,7 +93,13 @@ public:
}
VTKM_EXEC
bool CanContinue(const vtkm::Id& idx) { return this->GetParticle(idx).GetStatus().CanContinue(); }
bool CanContinue(const vtkm::Id& idx)
{
ParticleType particle(this->GetParticle(idx));
auto terminate = this->Termination.CheckTermination(particle);
this->Particles.Set(idx, particle);
return terminate;
}
VTKM_EXEC
void UpdateTookSteps(const vtkm::Id& idx, bool val)
@ -103,176 +114,48 @@ public:
protected:
using ParticlePortal = typename vtkm::cont::ArrayHandle<ParticleType>::WritePortalType;
ParticlePortal Particles;
vtkm::Id MaxSteps;
TerminationType Termination;
AnalysisType Analysis;
};
template <typename ParticleType>
template <typename ParticleType, typename TerminationType, typename AnalysisType>
class Particles : public vtkm::cont::ExecutionObjectBase
{
protected:
vtkm::cont::ArrayHandle<ParticleType> ParticleArray;
TerminationType Termination;
AnalysisType Analysis;
public:
VTKM_CONT vtkm::worklet::flow::ParticleExecutionObject<ParticleType> PrepareForExecution(
vtkm::cont::DeviceAdapterId device,
vtkm::cont::Token& token) const
{
return vtkm::worklet::flow::ParticleExecutionObject<ParticleType>(
this->ParticleArray, this->MaxSteps, device, token);
}
VTKM_CONT
Particles() = default;
VTKM_CONT
Particles(vtkm::cont::ArrayHandle<ParticleType>& pArray, vtkm::Id& maxSteps)
Particles(vtkm::cont::ArrayHandle<ParticleType>& pArray,
TerminationType termination,
AnalysisType analysis)
: ParticleArray(pArray)
, MaxSteps(maxSteps)
, Termination(termination)
, Analysis(analysis)
{
}
Particles() {}
protected:
vtkm::cont::ArrayHandle<ParticleType> ParticleArray;
vtkm::Id MaxSteps;
};
template <typename ParticleType>
class StateRecordingParticleExecutionObject : public ParticleExecutionObject<ParticleType>
VTKM_CONT auto PrepareForExecution(vtkm::cont::DeviceAdapterId device,
vtkm::cont::Token& token) const
-> vtkm::worklet::flow::ParticleExecutionObject<
ParticleType,
decltype(this->Termination.PrepareForExecution(device, token)),
decltype(this->Analysis.PrepareForExecution(device, token))>
{
public:
VTKM_EXEC_CONT
StateRecordingParticleExecutionObject()
: ParticleExecutionObject<ParticleType>()
, History()
, Length(0)
, StepCount()
, ValidPoint()
{
}
StateRecordingParticleExecutionObject(vtkm::cont::ArrayHandle<ParticleType> pArray,
vtkm::cont::ArrayHandle<vtkm::Vec3f> historyArray,
vtkm::cont::ArrayHandle<vtkm::Id> validPointArray,
vtkm::cont::ArrayHandle<vtkm::Id> stepCountArray,
vtkm::Id maxSteps,
vtkm::cont::DeviceAdapterId device,
vtkm::cont::Token& token)
: ParticleExecutionObject<ParticleType>(pArray, maxSteps, device, token)
, Length(maxSteps + 1)
{
vtkm::Id numPos = pArray.GetNumberOfValues();
History = historyArray.PrepareForOutput(numPos * Length, device, token);
ValidPoint = validPointArray.PrepareForInPlace(device, token);
StepCount = stepCountArray.PrepareForInPlace(device, token);
}
VTKM_EXEC
void PreStepUpdate(const vtkm::Id& idx)
{
ParticleType p = this->ParticleExecutionObject<ParticleType>::GetParticle(idx);
if (this->StepCount.Get(idx) == 0)
{
vtkm::Id loc = idx * Length;
this->History.Set(loc, p.GetPosition());
this->ValidPoint.Set(loc, 1);
this->StepCount.Set(idx, 1);
}
}
VTKM_EXEC
void StepUpdate(const vtkm::Id& idx,
const ParticleType& particle,
vtkm::FloatDefault time,
const vtkm::Vec3f& pt)
{
this->ParticleExecutionObject<ParticleType>::StepUpdate(idx, particle, time, pt);
//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);
this->StepCount.Set(idx, stepCount + 1);
}
protected:
using IdPortal = typename vtkm::cont::ArrayHandle<vtkm::Id>::WritePortalType;
using HistoryPortal = typename vtkm::cont::ArrayHandle<vtkm::Vec3f>::WritePortalType;
HistoryPortal History;
vtkm::Id Length;
IdPortal StepCount;
IdPortal ValidPoint;
};
template <typename ParticleType>
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);
auto termination = this->Termination.PrepareForExecution(device, token);
auto analysis = this->Analysis.PrepareForExecution(device, token);
return vtkm::worklet::flow::
ParticleExecutionObject<ParticleType, decltype(termination), decltype(analysis)>(
this->ParticleArray, termination, analysis, device, token);
}
};
VTKM_CONT vtkm::worklet::flow::StateRecordingParticleExecutionObject<ParticleType>
PrepareForExecution(vtkm::cont::DeviceAdapterId device, vtkm::cont::Token& token) const
{
return vtkm::worklet::flow::StateRecordingParticleExecutionObject<ParticleType>(
this->ParticleArray,
this->HistoryArray,
this->ValidPointArray,
this->StepCountArray,
this->MaxSteps,
device,
token);
}
VTKM_CONT
StateRecordingParticles(vtkm::cont::ArrayHandle<ParticleType>& 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<ParticleType>& pArray,
vtkm::cont::ArrayHandle<vtkm::Vec3f>& historyArray,
vtkm::cont::ArrayHandle<vtkm::Id>& validPointArray,
vtkm::Id& maxSteps)
{
ParticleArray = pArray;
HistoryArray = historyArray;
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::Vec3f> HistoryArray;
vtkm::Id MaxSteps;
vtkm::cont::ArrayHandle<ParticleType> ParticleArray;
vtkm::cont::ArrayHandle<vtkm::Id> StepCountArray;
vtkm::cont::ArrayHandle<vtkm::Id> ValidPointArray;
};
}
}
} //vtkm::worklet::flow

2
vtkm/filter/flow/worklet/RK4Integrator.h
View File

@ -34,7 +34,7 @@ public:
}
template <typename Particle>
VTKM_EXEC IntegratorStatus CheckStep(Particle& particle,
VTKM_EXEC IntegratorStatus CheckStep(const Particle& particle,
vtkm::FloatDefault stepLength,
vtkm::Vec3f& velocity) const
{

173
vtkm/filter/flow/worklet/Termination.h Normal file
View File

@ -0,0 +1,173 @@
//=============================================================================
//
// 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 vtkm_worklet_particleadvection_termination
#define vtkm_worklet_particleadvection_termination
#include <vtkm/Types.h>
#include <vtkm/cont/ExecutionObjectBase.h>
namespace vtkm
{
namespace worklet
{
namespace flow
{
class NormalTerminationExec
{
public:
VTKM_EXEC_CONT
NormalTerminationExec()
: MaxSteps(0)
{
}
VTKM_EXEC_CONT
NormalTerminationExec(vtkm::Id maxSteps)
: MaxSteps(maxSteps)
{
}
template <typename ParticleType>
VTKM_EXEC bool CheckTermination(ParticleType& particle) const
{
/// Checks particle properties to make a decision for termination
/// -- Check if the particle is out of spatial boundaries
/// -- Check if the particle has reached the maximum number of steps
/// -- Check if the particle is in a zero velocity region
auto& status = particle.GetStatus();
if (particle.GetNumberOfSteps() == this->MaxSteps)
{
status.SetTerminate();
particle.SetStatus(status);
}
bool terminate = status.CheckOk() && !status.CheckTerminate() && !status.CheckSpatialBounds() &&
!status.CheckTemporalBounds() && !status.CheckInGhostCell() && !status.CheckZeroVelocity();
return terminate;
}
private:
vtkm::Id MaxSteps;
};
class NormalTermination : public vtkm::cont::ExecutionObjectBase
{
public:
VTKM_CONT
NormalTermination()
: MaxSteps(0)
{
}
VTKM_CONT
NormalTermination(vtkm::Id maxSteps)
: MaxSteps(maxSteps)
{
}
VTKM_CONT
vtkm::Id AllocationSize() const { return this->MaxSteps; }
VTKM_CONT vtkm::worklet::flow::NormalTerminationExec PrepareForExecution(
vtkm::cont::DeviceAdapterId device,
vtkm::cont::Token& token) const
{
(void)device;
(void)token;
return vtkm::worklet::flow::NormalTerminationExec(MaxSteps);
}
private:
vtkm::Id MaxSteps;
};
class PoincareTerminationExec
{
public:
VTKM_EXEC_CONT
PoincareTerminationExec()
: MaxSteps(0)
, MaxPunctures(0)
{
}
VTKM_EXEC_CONT
PoincareTerminationExec(vtkm::Id maxSteps, vtkm::Id maxPunctures)
: MaxSteps(maxSteps)
, MaxPunctures(maxPunctures)
{
}
template <typename ParticleType>
VTKM_EXEC bool CheckTermination(ParticleType& particle) const
{
/// Checks particle properties to make a decision for termination
/// -- Check if the particle is out of spatial boundaries
/// -- Check if the particle has reached the maximum number of steps
/// -- Check if the particle is in a zero velocity region
auto& status = particle.GetStatus();
if (particle.GetNumberOfSteps() >= this->MaxSteps ||
particle.GetNumberOfPunctures() >= this->MaxPunctures)
{
status.SetTerminate();
particle.SetStatus(status);
}
bool terminate = status.CheckOk() && !status.CheckTerminate() && !status.CheckSpatialBounds() &&
!status.CheckTemporalBounds() && !status.CheckInGhostCell() && !status.CheckZeroVelocity();
return terminate;
}
private:
vtkm::Id MaxSteps;
vtkm::Id MaxPunctures;
};
class PoincareTermination : public vtkm::cont::ExecutionObjectBase
{
public:
VTKM_CONT
PoincareTermination()
: MaxSteps(0)
, MaxPunctures(0)
{
}
VTKM_CONT
PoincareTermination(vtkm::Id maxSteps, vtkm::Id maxPunctures)
: MaxSteps(maxSteps)
, MaxPunctures(maxPunctures)
{
}
VTKM_CONT
vtkm::Id AllocationSize() const { return this->MaxPunctures; }
VTKM_CONT vtkm::worklet::flow::PoincareTerminationExec PrepareForExecution(
vtkm::cont::DeviceAdapterId device,
vtkm::cont::Token& token) const
{
(void)device;
(void)token;
return vtkm::worklet::flow::PoincareTerminationExec(MaxSteps, MaxPunctures);
}
private:
vtkm::Id MaxSteps;
vtkm::Id MaxPunctures;
};
} // namespace particleadvection
} // namespace worklet
} // namespace vtkm
#endif