diff --git a/vtkm/Particle.h b/vtkm/Particle.h index 32ae7b214..6e545446c 100644 --- a/vtkm/Particle.h +++ b/vtkm/Particle.h @@ -228,6 +228,14 @@ public: return this->Pos + translation; } + inline VTKM_CONT friend std::ostream& operator<<(std::ostream& out, + const vtkm::ChargedParticle& p) + { + out << "v(" << p.Time << ") = " << p.Pos << ", ID: " << p.ID << ", NumSteps: " << p.NumSteps + << ", Status: " << p.Status; + return out; + } + vtkm::Vec3f Pos; vtkm::Id ID = -1; vtkm::Id NumSteps = 0; diff --git a/vtkm/filter/NewFilterParticleAdvection.h b/vtkm/filter/NewFilterParticleAdvection.h new file mode 100644 index 000000000..a9e41f0c9 --- /dev/null +++ b/vtkm/filter/NewFilterParticleAdvection.h @@ -0,0 +1,148 @@ +//============================================================================ +// 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_NewFilterParticleAdvection_h +#define vtk_m_filter_NewFilterParticleAdvection_h + +#include +#include +#include +#include +#include +#include +#include + +namespace vtkm +{ +namespace filter +{ +/// \brief base class for advecting particles in a vector field. + +/// Takes as input a vector field and seed locations and advects the seeds +/// through the flow field. + +class NewFilterParticleAdvection : public vtkm::filter::NewFilterField +{ +public: + using SupportedTypes = vtkm::TypeListFieldVec3; + + VTKM_CONT + NewFilterParticleAdvection(vtkm::filter::particleadvection::ParticleAdvectionResultType rType); + + VTKM_CONT + void SetStepSize(vtkm::FloatDefault s) { this->StepSize = s; } + + VTKM_CONT + void SetNumberOfSteps(vtkm::Id n) { this->NumberOfSteps = n; } + + template + VTKM_CONT void SetSeeds(vtkm::cont::ArrayHandle& seeds) + { + this->Seeds = seeds; + } + + template + VTKM_CONT void SetSeeds(const std::vector& seeds, + vtkm::CopyFlag copyFlag = vtkm::CopyFlag::On) + { + this->Seeds = vtkm::cont::make_ArrayHandle(seeds, copyFlag); + } + + VTKM_CONT + void SetSolverRK4() + { + this->SolverType = vtkm::filter::particleadvection::IntegrationSolverType::RK4_TYPE; + } + VTKM_CONT + void SetSolverEuler() + { + this->SolverType = vtkm::filter::particleadvection::IntegrationSolverType::EULER_TYPE; + } + + VTKM_CONT + bool GetUseThreadedAlgorithm() { return this->UseThreadedAlgorithm; } + + VTKM_CONT + void SetUseThreadedAlgorithm(bool val) { this->UseThreadedAlgorithm = val; } + +protected: + VTKM_CONT vtkm::cont::DataSet DoExecute(const vtkm::cont::DataSet& inData) override + { + auto out = this->DoExecutePartitions(inData); + if (out.GetNumberOfPartitions() != 1) + throw vtkm::cont::ErrorFilterExecution("Wrong number of results"); + + return out.GetPartition(0); + } + + VTKM_CONT virtual void ValidateOptions() const; + + vtkm::Id NumberOfSteps; + vtkm::filter::particleadvection::ParticleAdvectionResultType ResultType = + vtkm::filter::particleadvection::ParticleAdvectionResultType::UNKNOWN_TYPE; + vtkm::cont::UnknownArrayHandle Seeds; + vtkm::filter::particleadvection::IntegrationSolverType SolverType; + vtkm::FloatDefault StepSize; + bool UseThreadedAlgorithm; + vtkm::filter::particleadvection::VectorFieldType VecFieldType; + +private: +}; + + +class NewFilterSteadyStateParticleAdvection : public NewFilterParticleAdvection +{ +public: + VTKM_CONT + NewFilterSteadyStateParticleAdvection( + vtkm::filter::particleadvection::ParticleAdvectionResultType rType) + : NewFilterParticleAdvection(rType) + { + } + +protected: + VTKM_CONT std::vector CreateDataSetIntegrators( + const vtkm::cont::PartitionedDataSet& input, + const vtkm::filter::particleadvection::BoundsMap& boundsMap) const; +}; + +class NewFilterUnsteadyStateParticleAdvection : public NewFilterParticleAdvection +{ +public: + VTKM_CONT + NewFilterUnsteadyStateParticleAdvection( + vtkm::filter::particleadvection::ParticleAdvectionResultType rType) + : NewFilterParticleAdvection(rType) + { + } + + void SetPreviousTime(vtkm::FloatDefault t1) { this->Time1 = t1; } + void SetNextTime(vtkm::FloatDefault t2) { this->Time2 = t2; } + void SetNextDataSet(const vtkm::cont::DataSet& ds) { this->Input2 = { ds }; } + void SetNextDataSet(const vtkm::cont::PartitionedDataSet& pds) { this->Input2 = pds; } + +protected: + VTKM_CONT std::vector + CreateDataSetIntegrators(const vtkm::cont::PartitionedDataSet& input, + const vtkm::filter::particleadvection::BoundsMap& boundsMap) const; + + vtkm::cont::PartitionedDataSet Input2; + vtkm::FloatDefault Time1; + vtkm::FloatDefault Time2; +}; + +} +} // namespace vtkm::filter + +#ifndef vtk_m_filter_NewFilterParticleAdvection_hxx +#include +#endif + +#endif // vtk_m_filter_NewFilterParticleAdvection_h diff --git a/vtkm/filter/NewFilterParticleAdvection.hxx b/vtkm/filter/NewFilterParticleAdvection.hxx new file mode 100644 index 000000000..e3abdd032 --- /dev/null +++ b/vtkm/filter/NewFilterParticleAdvection.hxx @@ -0,0 +1,104 @@ +//============================================================================ +// 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_NewFilterParticleAdvection_hxx +#define vtk_m_filter_FilterParticleAdvection_hxx + +#include + +namespace vtkm +{ +namespace filter +{ + +//----------------------------------------------------------------------------- +inline VTKM_CONT NewFilterParticleAdvection::NewFilterParticleAdvection( + vtkm::filter::particleadvection::ParticleAdvectionResultType rType) + : vtkm::filter::NewFilterField() + , NumberOfSteps(0) + , ResultType(rType) + , SolverType(vtkm::filter::particleadvection::IntegrationSolverType::RK4_TYPE) + , StepSize(0) + , UseThreadedAlgorithm(false) + , VecFieldType(vtkm::filter::particleadvection::VectorFieldType::VELOCITY_FIELD_TYPE) +{ +} + +void NewFilterParticleAdvection::ValidateOptions() const +{ + if (this->GetUseCoordinateSystemAsField()) + throw vtkm::cont::ErrorFilterExecution("Coordinate system as field not supported"); + if (this->Seeds.GetNumberOfValues() == 0) + throw vtkm::cont::ErrorFilterExecution("No seeds provided."); + if (this->NumberOfSteps == 0) + throw vtkm::cont::ErrorFilterExecution("Number of steps not specified."); + if (this->StepSize == 0) + throw vtkm::cont::ErrorFilterExecution("Step size not specified."); +} + +VTKM_CONT std::vector +NewFilterSteadyStateParticleAdvection::CreateDataSetIntegrators( + const vtkm::cont::PartitionedDataSet& input, + const vtkm::filter::particleadvection::BoundsMap& boundsMap) const +{ + using DSIType = vtkm::filter::particleadvection::DSISteadyState; + + std::string activeField = this->GetActiveFieldName(); + + std::vector dsi; + for (vtkm::Id i = 0; i < input.GetNumberOfPartitions(); i++) + { + vtkm::Id blockId = boundsMap.GetLocalBlockId(i); + auto ds = input.GetPartition(i); + if (!ds.HasPointField(activeField) && !ds.HasCellField(activeField)) + throw vtkm::cont::ErrorFilterExecution("Unsupported field assocation"); + + dsi.push_back(new DSIType( + ds, blockId, activeField, this->SolverType, this->VecFieldType, this->ResultType)); + } + + return dsi; +} + +VTKM_CONT std::vector +NewFilterUnsteadyStateParticleAdvection::CreateDataSetIntegrators( + const vtkm::cont::PartitionedDataSet& input, + const vtkm::filter::particleadvection::BoundsMap& boundsMap) const +{ + using DSIType = vtkm::filter::particleadvection::DSIUnsteadyState; + + std::string activeField = this->GetActiveFieldName(); + + std::vector dsi; + for (vtkm::Id i = 0; i < input.GetNumberOfPartitions(); i++) + { + vtkm::Id blockId = boundsMap.GetLocalBlockId(i); + auto ds1 = input.GetPartition(i); + auto ds2 = this->Input2.GetPartition(i); + if ((!ds1.HasPointField(activeField) && !ds1.HasCellField(activeField)) || + (!ds2.HasPointField(activeField) && !ds2.HasCellField(activeField))) + throw vtkm::cont::ErrorFilterExecution("Unsupported field assocation"); + + dsi.push_back(new DSIType(ds1, + ds2, + this->Time1, + this->Time2, + blockId, + activeField, + this->SolverType, + this->VecFieldType, + this->ResultType)); + } + + return dsi; +} + +} +} // namespace vtkm::filter +#endif diff --git a/vtkm/filter/ParticleAdvection2.h b/vtkm/filter/ParticleAdvection2.h index 6fad92ed9..8c8ab9aa9 100644 --- a/vtkm/filter/ParticleAdvection2.h +++ b/vtkm/filter/ParticleAdvection2.h @@ -12,8 +12,8 @@ #define vtk_m_filter_ParticleAdvection2_h #include -#include #include +#include #include namespace vtkm @@ -25,31 +25,18 @@ namespace filter /// Takes as input a vector field and seed locations and generates the /// end points for each seed through the vector field. -class ParticleAdvection2 - : public vtkm::filter::FilterParticleAdvection +class ParticleAdvection2 : public vtkm::filter::NewFilterSteadyStateParticleAdvection { public: - VTKM_CONT ParticleAdvection2(); - - template - vtkm::cont::PartitionedDataSet PrepareForExecution( - const vtkm::cont::PartitionedDataSet& input, - const vtkm::filter::PolicyBase& policy); - - vtkm::cont::UnknownArrayHandle SeedArray; -}; - -class ParticleAdvection3 : public vtkm::filter::NewFilterField -{ -public: - // VTKM_CONT ParticleAdvection3() {} + VTKM_CONT ParticleAdvection2() + : NewFilterSteadyStateParticleAdvection( + vtkm::filter::particleadvection::ParticleAdvectionResultType::PARTICLE_ADVECT_TYPE) + { + } protected: - VTKM_CONT vtkm::cont::DataSet DoExecute(const vtkm::cont::DataSet& inData) override; VTKM_CONT vtkm::cont::PartitionedDataSet DoExecutePartitions( const vtkm::cont::PartitionedDataSet& inData) override; - - vtkm::cont::UnknownArrayHandle SeedArray; }; } diff --git a/vtkm/filter/ParticleAdvection2.hxx b/vtkm/filter/ParticleAdvection2.hxx index eb4903e7c..16fba05ef 100644 --- a/vtkm/filter/ParticleAdvection2.hxx +++ b/vtkm/filter/ParticleAdvection2.hxx @@ -13,8 +13,6 @@ #include #include #include -#include -#include #include #include @@ -25,63 +23,19 @@ namespace vtkm namespace filter { -//----------------------------------------------------------------------------- -inline VTKM_CONT ParticleAdvection2::ParticleAdvection2() - : vtkm::filter::FilterParticleAdvection() +VTKM_CONT vtkm::cont::PartitionedDataSet ParticleAdvection2::DoExecutePartitions( + const vtkm::cont::PartitionedDataSet& input) { - this->ResultType = - vtkm::filter::particleadvection::ParticleAdvectionResultType::PARTICLE_ADVECT_TYPE; -} - -//----------------------------------------------------------------------------- -template -inline VTKM_CONT vtkm::cont::PartitionedDataSet ParticleAdvection2::PrepareForExecution( - const vtkm::cont::PartitionedDataSet& input, - const vtkm::filter::PolicyBase&) -{ - // return input; - -#if 1 using DSIType = vtkm::filter::particleadvection::DSISteadyState; - this->ValidateOptions(); - //Make sure everything matches up ok. - this->VecFieldType = vtkm::filter::particleadvection::VectorFieldType::VELOCITY_FIELD_TYPE; vtkm::filter::particleadvection::BoundsMap boundsMap(input); - std::string activeField = this->GetActiveFieldName(); + auto dsi = this->CreateDataSetIntegrators(input, boundsMap); - std::vector dsi; - for (vtkm::Id i = 0; i < input.GetNumberOfPartitions(); i++) - { - vtkm::Id blockId = boundsMap.GetLocalBlockId(i); - auto ds = input.GetPartition(i); - if (!ds.HasPointField(activeField) && !ds.HasCellField(activeField)) - throw vtkm::cont::ErrorFilterExecution("Unsupported field assocation"); - - dsi.push_back(new DSIType( - ds, blockId, activeField, this->SolverType, this->VecFieldType, this->ResultType)); - } - - this->SeedArray = this->Seeds; vtkm::filter::particleadvection::PAV pav( boundsMap, dsi, this->UseThreadedAlgorithm, this->ResultType); - return pav.Execute(this->NumberOfSteps, this->StepSize, this->SeedArray); -#endif -} -VTKM_CONT vtkm::cont::DataSet ParticleAdvection3::DoExecute(const vtkm::cont::DataSet& inData) -{ - std::cout << "Meow DS" << std::endl; - auto result = this->DoExecutePartitions(inData); - return result.GetPartition(0); -} - -VTKM_CONT vtkm::cont::PartitionedDataSet ParticleAdvection3::DoExecutePartitions( - const vtkm::cont::PartitionedDataSet& inData) -{ - std::cout << "Meow pDS" << std::endl; - return inData; + return pav.Execute(this->NumberOfSteps, this->StepSize, this->Seeds); } } diff --git a/vtkm/filter/Pathline2.h b/vtkm/filter/Pathline2.h index 35f49799d..cd70bcf3c 100644 --- a/vtkm/filter/Pathline2.h +++ b/vtkm/filter/Pathline2.h @@ -12,8 +12,8 @@ #define vtk_m_filter_Pathline2_h #include -#include #include +#include #include namespace vtkm @@ -25,39 +25,18 @@ namespace filter /// Takes as input a vector field and seed locations and generates the /// end points for each seed through the vector field. -class Pathline2 : public vtkm::filter::FilterTemporalParticleAdvection +class Pathline2 : public vtkm::filter::NewFilterUnsteadyStateParticleAdvection { public: - VTKM_CONT Pathline2(); - - template - vtkm::cont::PartitionedDataSet PrepareForExecution( - const vtkm::cont::PartitionedDataSet& input, - const vtkm::filter::PolicyBase& policy); - - vtkm::cont::UnknownArrayHandle SeedArray; - - void SetPreviousTime(vtkm::FloatDefault t1) { this->Time1 = t1; } - void SetNextTime(vtkm::FloatDefault t2) { this->Time2 = t2; } - void SetNextDataSet(const vtkm::cont::DataSet& ds) { this->DataSet2 = { ds }; } - void SetNextDataSet(const vtkm::cont::PartitionedDataSet& pds) { this->DataSet2 = pds; } - - vtkm::cont::PartitionedDataSet DataSet2; - vtkm::FloatDefault Time1; - vtkm::FloatDefault Time2; -}; - -class Pathline3 : public vtkm::filter::NewFilterField -{ -public: - // VTKM_CONT Pathline3() {} + VTKM_CONT Pathline2() + : NewFilterUnsteadyStateParticleAdvection( + vtkm::filter::particleadvection::ParticleAdvectionResultType::STREAMLINE_TYPE) + { + } protected: - VTKM_CONT vtkm::cont::DataSet DoExecute(const vtkm::cont::DataSet& inData) override; VTKM_CONT vtkm::cont::PartitionedDataSet DoExecutePartitions( const vtkm::cont::PartitionedDataSet& inData) override; - - vtkm::cont::UnknownArrayHandle SeedArray; }; } diff --git a/vtkm/filter/Pathline2.hxx b/vtkm/filter/Pathline2.hxx index 2bb7e3d7c..106642a18 100644 --- a/vtkm/filter/Pathline2.hxx +++ b/vtkm/filter/Pathline2.hxx @@ -25,70 +25,19 @@ namespace vtkm namespace filter { -//----------------------------------------------------------------------------- -inline VTKM_CONT Pathline2::Pathline2() - : vtkm::filter::FilterTemporalParticleAdvection() +VTKM_CONT vtkm::cont::PartitionedDataSet Pathline2::DoExecutePartitions( + const vtkm::cont::PartitionedDataSet& input) { - this->ResultType = vtkm::filter::particleadvection::ParticleAdvectionResultType::STREAMLINE_TYPE; -} - -//----------------------------------------------------------------------------- -template -inline VTKM_CONT vtkm::cont::PartitionedDataSet Pathline2::PrepareForExecution( - const vtkm::cont::PartitionedDataSet& input, - const vtkm::filter::PolicyBase&) -{ -// return input; -#if 1 using DSIType = vtkm::filter::particleadvection::DSIUnsteadyState; - this->ValidateOptions(); - //Make sure everything matches up ok. - this->VecFieldType = vtkm::filter::particleadvection::VectorFieldType::VELOCITY_FIELD_TYPE; vtkm::filter::particleadvection::BoundsMap boundsMap(input); - std::string activeField = this->GetActiveFieldName(); + auto dsi = this->CreateDataSetIntegrators(input, boundsMap); - std::vector dsi; - for (vtkm::Id i = 0; i < input.GetNumberOfPartitions(); i++) - { - vtkm::Id blockId = boundsMap.GetLocalBlockId(i); - auto ds1 = input.GetPartition(i); - auto ds2 = this->DataSet2.GetPartition(i); - if ((!ds1.HasPointField(activeField) && !ds1.HasCellField(activeField)) || - (!ds2.HasPointField(activeField) && !ds2.HasCellField(activeField))) - throw vtkm::cont::ErrorFilterExecution("Unsupported field assocation"); - - dsi.push_back(new DSIType(ds1, - ds2, - this->Time1, - this->Time2, - blockId, - activeField, - this->SolverType, - this->VecFieldType, - this->ResultType)); - } - - this->SeedArray = this->Seeds; vtkm::filter::particleadvection::PAV pav( boundsMap, dsi, this->UseThreadedAlgorithm, this->ResultType); - return pav.Execute(this->NumberOfSteps, this->StepSize, this->SeedArray); -#endif -} -VTKM_CONT vtkm::cont::DataSet Pathline3::DoExecute(const vtkm::cont::DataSet& inData) -{ - std::cout << "Meow DS" << std::endl; - auto result = this->DoExecutePartitions(inData); - return result.GetPartition(0); -} - -VTKM_CONT vtkm::cont::PartitionedDataSet Pathline3::DoExecutePartitions( - const vtkm::cont::PartitionedDataSet& inData) -{ - std::cout << "Meow pDS" << std::endl; - return inData; + return pav.Execute(this->NumberOfSteps, this->StepSize, this->Seeds); } } diff --git a/vtkm/filter/Streamline2.h b/vtkm/filter/Streamline2.h index f369753bc..59b9cba95 100644 --- a/vtkm/filter/Streamline2.h +++ b/vtkm/filter/Streamline2.h @@ -12,8 +12,8 @@ #define vtk_m_filter_Streamline2_h #include -#include #include +#include #include namespace vtkm @@ -25,30 +25,18 @@ namespace filter /// Takes as input a vector field and seed locations and generates the /// end points for each seed through the vector field. -class Streamline2 : public vtkm::filter::FilterParticleAdvection +class Streamline2 : public vtkm::filter::NewFilterSteadyStateParticleAdvection { public: - VTKM_CONT Streamline2(); - - template - vtkm::cont::PartitionedDataSet PrepareForExecution( - const vtkm::cont::PartitionedDataSet& input, - const vtkm::filter::PolicyBase& policy); - - vtkm::cont::UnknownArrayHandle SeedArray; -}; - -class Streamline3 : public vtkm::filter::NewFilterField -{ -public: - // VTKM_CONT Streamline3() {} + VTKM_CONT Streamline2() + : NewFilterSteadyStateParticleAdvection( + vtkm::filter::particleadvection::ParticleAdvectionResultType::STREAMLINE_TYPE) + { + } protected: - VTKM_CONT vtkm::cont::DataSet DoExecute(const vtkm::cont::DataSet& inData) override; VTKM_CONT vtkm::cont::PartitionedDataSet DoExecutePartitions( const vtkm::cont::PartitionedDataSet& inData) override; - - vtkm::cont::UnknownArrayHandle SeedArray; }; } diff --git a/vtkm/filter/Streamline2.hxx b/vtkm/filter/Streamline2.hxx index ebcbe90c0..aeff1569b 100644 --- a/vtkm/filter/Streamline2.hxx +++ b/vtkm/filter/Streamline2.hxx @@ -13,8 +13,6 @@ #include #include #include -#include -#include #include #include @@ -25,62 +23,19 @@ namespace vtkm namespace filter { -//----------------------------------------------------------------------------- -inline VTKM_CONT Streamline2::Streamline2() - : vtkm::filter::FilterParticleAdvection() +VTKM_CONT vtkm::cont::PartitionedDataSet Streamline2::DoExecutePartitions( + const vtkm::cont::PartitionedDataSet& input) { - this->ResultType = vtkm::filter::particleadvection::ParticleAdvectionResultType::STREAMLINE_TYPE; -} - -//----------------------------------------------------------------------------- -template -inline VTKM_CONT vtkm::cont::PartitionedDataSet Streamline2::PrepareForExecution( - const vtkm::cont::PartitionedDataSet& input, - const vtkm::filter::PolicyBase&) -{ - // return input; - -#if 1 using DSIType = vtkm::filter::particleadvection::DSISteadyState; - this->ValidateOptions(); - //Make sure everything matches up ok. - this->VecFieldType = vtkm::filter::particleadvection::VectorFieldType::VELOCITY_FIELD_TYPE; vtkm::filter::particleadvection::BoundsMap boundsMap(input); - std::string activeField = this->GetActiveFieldName(); + auto dsi = this->CreateDataSetIntegrators(input, boundsMap); - std::vector dsi; - for (vtkm::Id i = 0; i < input.GetNumberOfPartitions(); i++) - { - vtkm::Id blockId = boundsMap.GetLocalBlockId(i); - auto ds = input.GetPartition(i); - if (!ds.HasPointField(activeField) && !ds.HasCellField(activeField)) - throw vtkm::cont::ErrorFilterExecution("Unsupported field assocation"); - - dsi.push_back(new DSIType( - ds, blockId, activeField, this->SolverType, this->VecFieldType, this->ResultType)); - } - - this->SeedArray = this->Seeds; vtkm::filter::particleadvection::PAV pav( boundsMap, dsi, this->UseThreadedAlgorithm, this->ResultType); - return pav.Execute(this->NumberOfSteps, this->StepSize, this->SeedArray); -#endif -} -VTKM_CONT vtkm::cont::DataSet Streamline3::DoExecute(const vtkm::cont::DataSet& inData) -{ - std::cout << "Meow DS" << std::endl; - auto result = this->DoExecutePartitions(inData); - return result.GetPartition(0); -} - -VTKM_CONT vtkm::cont::PartitionedDataSet Streamline3::DoExecutePartitions( - const vtkm::cont::PartitionedDataSet& inData) -{ - std::cout << "Meow pDS" << std::endl; - return inData; + return pav.Execute(this->NumberOfSteps, this->StepSize, this->Seeds); } }