Shuenhoy/vtk-m
1
0
Fork 0
mirror of https://gitlab.kitware.com/vtk/vtk-m synced 2024-10-05 09:59:12 +00:00
Code Issues 2 Actions Packages Projects Releases Wiki Activity

Expand support for DSI

Browse Source
This commit is contained in:
Dave Pugmire 2022-06-30 09:54:43 -04:00
parent c40356bc8a
commit db4d8a181e
11 changed files with 384 additions and 185 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

35
vtkm/filter/ParticleAdvection2.hxx
View File

@ -38,8 +38,9 @@ inline VTKM_CONT vtkm::cont::PartitionedDataSet ParticleAdvection2::PrepareForEx
const vtkm::cont::PartitionedDataSet& input,
const vtkm::filter::PolicyBase<DerivedPolicy>&)
{
// using AlgorithmType = vtkm::filter::particleadvection::ParticleAdvectionAlgorithm;
// using ThreadedAlgorithmType = vtkm::filter::particleadvection::ParticleAdvectionThreadedAlgorithm;
// return input;
#if 1
using DSIType = vtkm::filter::particleadvection::DSI;
this->ValidateOptions();
@ -62,35 +63,9 @@ inline VTKM_CONT vtkm::cont::PartitionedDataSet ParticleAdvection2::PrepareForEx
}
this->SeedArray = this->Seeds;
vtkm::filter::particleadvection::PAV pav(
boundsMap,
dsi,
this->UseThreadedAlgorithm,
vtkm::filter::particleadvection::ParticleAdvectionResultType::PARTICLE_ADVECT_TYPE);
vtkm::filter::particleadvection::PAV<DSIType> pav(
boundsMap, dsi, this->UseThreadedAlgorithm, this->ResultType);
return pav.Execute(this->NumberOfSteps, this->StepSize, this->SeedArray);
#if 0
//std::vector<DSIType> ddsi;
/*
vtkm::filter::particleadvection::RunAlgo<DSIType, AlgorithmType>(
boundsMap, ddsi, this->NumberOfSteps, this->StepSize, this->Seeds);
*/
vtkm::cont::PartitionedDataSet output;
return output;
/*
//using DSIType = vtkm::filter::particleadvection::DataSetIntegrator;
//std::vector<DSIType> dsi;
auto dsi = this->CreateDataSetIntegrators(input, boundsMap);
if (this->GetUseThreadedAlgorithm())
return vtkm::filter::particleadvection::RunAlgo<DSIType, ThreadedAlgorithmType>(
boundsMap, dsi, this->NumberOfSteps, this->StepSize, this->Seeds);
else
return vtkm::filter::particleadvection::RunAlgo<DSIType, AlgorithmType>(
boundsMap, dsi, this->NumberOfSteps, this->StepSize, this->Seeds);
*/
#endif
}

9
vtkm/filter/PathParticle2.h
View File

@ -37,6 +37,15 @@ public:
const vtkm::filter::PolicyBase<DerivedPolicy>& 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 PathParticle3 : public vtkm::filter::NewFilterField

45
vtkm/filter/PathParticle2.hxx
View File

@ -38,8 +38,6 @@ inline VTKM_CONT vtkm::cont::PartitionedDataSet PathParticle2::PrepareForExecuti
const vtkm::cont::PartitionedDataSet& input,
const vtkm::filter::PolicyBase<DerivedPolicy>&)
{
// using AlgorithmType = vtkm::filter::particleadvection::ParticleAdvectionAlgorithm;
// using ThreadedAlgorithmType = vtkm::filter::particleadvection::ParticleAdvectionThreadedAlgorithm;
using DSIType = vtkm::filter::particleadvection::DSI;
this->ValidateOptions();
@ -53,42 +51,27 @@ inline VTKM_CONT vtkm::cont::PartitionedDataSet PathParticle2::PrepareForExecuti
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))
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(
ds, blockId, activeField, this->SolverType, this->VecFieldType, this->ResultType));
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(
vtkm::filter::particleadvection::PAV<DSIType> pav(
boundsMap, dsi, this->UseThreadedAlgorithm, this->ResultType);
return pav.Execute(this->NumberOfSteps, this->StepSize, this->SeedArray);
#if 0
//std::vector<DSIType> ddsi;
/*
vtkm::filter::particleadvection::RunAlgo<DSIType, AlgorithmType>(
boundsMap, ddsi, this->NumberOfSteps, this->StepSize, this->Seeds);
*/
vtkm::cont::PartitionedDataSet output;
return output;
/*
//using DSIType = vtkm::filter::particleadvection::DataSetIntegrator;
//std::vector<DSIType> dsi;
auto dsi = this->CreateDataSetIntegrators(input, boundsMap);
if (this->GetUseThreadedAlgorithm())
return vtkm::filter::particleadvection::RunAlgo<DSIType, ThreadedAlgorithmType>(
boundsMap, dsi, this->NumberOfSteps, this->StepSize, this->Seeds);
else
return vtkm::filter::particleadvection::RunAlgo<DSIType, AlgorithmType>(
boundsMap, dsi, this->NumberOfSteps, this->StepSize, this->Seeds);
*/
#endif
}
VTKM_CONT vtkm::cont::DataSet PathParticle3::DoExecute(const vtkm::cont::DataSet& inData)

9
vtkm/filter/Pathline2.h
View File

@ -36,6 +36,15 @@ public:
const vtkm::filter::PolicyBase<DerivedPolicy>& 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

45
vtkm/filter/Pathline2.hxx
View File

@ -38,8 +38,6 @@ inline VTKM_CONT vtkm::cont::PartitionedDataSet Pathline2::PrepareForExecution(
const vtkm::cont::PartitionedDataSet& input,
const vtkm::filter::PolicyBase<DerivedPolicy>&)
{
// using AlgorithmType = vtkm::filter::particleadvection::ParticleAdvectionAlgorithm;
// using ThreadedAlgorithmType = vtkm::filter::particleadvection::ParticleAdvectionThreadedAlgorithm;
using DSIType = vtkm::filter::particleadvection::DSI;
this->ValidateOptions();
@ -53,42 +51,27 @@ inline VTKM_CONT vtkm::cont::PartitionedDataSet Pathline2::PrepareForExecution(
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))
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(
ds, blockId, activeField, this->SolverType, this->VecFieldType, this->ResultType));
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(
vtkm::filter::particleadvection::PAV<DSIType> pav(
boundsMap, dsi, this->UseThreadedAlgorithm, this->ResultType);
return pav.Execute(this->NumberOfSteps, this->StepSize, this->SeedArray);
#if 0
//std::vector<DSIType> ddsi;
/*
vtkm::filter::particleadvection::RunAlgo<DSIType, AlgorithmType>(
boundsMap, ddsi, this->NumberOfSteps, this->StepSize, this->Seeds);
*/
vtkm::cont::PartitionedDataSet output;
return output;
/*
//using DSIType = vtkm::filter::particleadvection::DataSetIntegrator;
//std::vector<DSIType> dsi;
auto dsi = this->CreateDataSetIntegrators(input, boundsMap);
if (this->GetUseThreadedAlgorithm())
return vtkm::filter::particleadvection::RunAlgo<DSIType, ThreadedAlgorithmType>(
boundsMap, dsi, this->NumberOfSteps, this->StepSize, this->Seeds);
else
return vtkm::filter::particleadvection::RunAlgo<DSIType, AlgorithmType>(
boundsMap, dsi, this->NumberOfSteps, this->StepSize, this->Seeds);
*/
#endif
}
VTKM_CONT vtkm::cont::DataSet Pathline3::DoExecute(const vtkm::cont::DataSet& inData)

30
vtkm/filter/Streamline2.hxx
View File

@ -38,8 +38,9 @@ inline VTKM_CONT vtkm::cont::PartitionedDataSet Streamline2::PrepareForExecution
const vtkm::cont::PartitionedDataSet& input,
const vtkm::filter::PolicyBase<DerivedPolicy>&)
{
// using AlgorithmType = vtkm::filter::particleadvection::StreamlineAlgorithm;
// using ThreadedAlgorithmType = vtkm::filter::particleadvection::StreamlineThreadedAlgorithm;
// return input;
#if 1
using DSIType = vtkm::filter::particleadvection::DSI;
this->ValidateOptions();
@ -62,32 +63,9 @@ inline VTKM_CONT vtkm::cont::PartitionedDataSet Streamline2::PrepareForExecution
}
this->SeedArray = this->Seeds;
vtkm::filter::particleadvection::PAV pav(
vtkm::filter::particleadvection::PAV<DSIType> pav(
boundsMap, dsi, this->UseThreadedAlgorithm, this->ResultType);
return pav.Execute(this->NumberOfSteps, this->StepSize, this->SeedArray);
#if 0
//std::vector<DSIType> ddsi;
/*
vtkm::filter::particleadvection::RunAlgo<DSIType, AlgorithmType>(
boundsMap, ddsi, this->NumberOfSteps, this->StepSize, this->Seeds);
*/
vtkm::cont::PartitionedDataSet output;
return output;
/*
//using DSIType = vtkm::filter::particleadvection::DataSetIntegrator;
//std::vector<DSIType> dsi;
auto dsi = this->CreateDataSetIntegrators(input, boundsMap);
if (this->GetUseThreadedAlgorithm())
return vtkm::filter::particleadvection::RunAlgo<DSIType, ThreadedAlgorithmType>(
boundsMap, dsi, this->NumberOfSteps, this->StepSize, this->Seeds);
else
return vtkm::filter::particleadvection::RunAlgo<DSIType, AlgorithmType>(
boundsMap, dsi, this->NumberOfSteps, this->StepSize, this->Seeds);
*/
#endif
}

6
vtkm/filter/particleadvection/ABA.h
View File

@ -24,11 +24,11 @@ namespace filter
namespace particleadvection
{
template <template <typename> class ResultType, typename ParticleType>
template <typename DSIType, template <typename> class ResultType, typename ParticleType>
class ABA
{
public:
ABA(const vtkm::filter::particleadvection::BoundsMap& bm, std::vector<DSI*>& blocks)
ABA(const vtkm::filter::particleadvection::BoundsMap& bm, std::vector<DSIType*>& blocks)
: Blocks(blocks)
, BoundsMap(bm)
, NumRanks(this->Comm.size())
@ -393,7 +393,7 @@ public:
//Member data
std::vector<ParticleType> Active;
std::vector<DSI*> Blocks;
std::vector<DSIType*> Blocks;
vtkm::filter::particleadvection::BoundsMap BoundsMap;
vtkmdiy::mpi::communicator Comm = vtkm::cont::EnvironmentTracker::GetCommunicator();
std::vector<ParticleType> Inactive;

121
vtkm/filter/particleadvection/DSI.h
View File

@ -54,18 +54,47 @@ public:
vtkm::filter::particleadvection::IntegrationSolverType solverType,
vtkm::filter::particleadvection::VectorFieldType vecFieldType,
vtkm::filter::particleadvection::ParticleAdvectionResultType resultType)
: Rank(this->Comm.rank())
: Data(ds)
, FieldName(fieldNm)
, Id(id)
, SolverType(solverType)
, Rank(this->Comm.rank())
, ResType(resultType)
, VecFieldType(vecFieldType)
{
this->SolverType = solverType;
this->VecFieldType = vecFieldType;
this->ResType = resultType;
this->DataSet = ds;
this->Id = id;
this->FieldName = fieldNm;
//check that things are valid.
}
DSI(const vtkm::cont::DataSet& ds1,
const vtkm::cont::DataSet& ds2,
vtkm::FloatDefault t1,
vtkm::FloatDefault t2,
vtkm::Id id,
const std::string& fieldNm,
vtkm::filter::particleadvection::IntegrationSolverType solverType,
vtkm::filter::particleadvection::VectorFieldType vecFieldType,
vtkm::filter::particleadvection::ParticleAdvectionResultType resultType)
: FieldName(fieldNm)
, Id(id)
, SolverType(solverType)
, Rank(this->Comm.rank())
, ResType(resultType)
, VecFieldType(vecFieldType)
{
std::cout << "****** Create Unsteady DSI." << std::endl;
this->Data = UnsteadyStateDataType(ds1, ds2, t1, t2);
}
VTKM_CONT bool IsSteadyState() const
{
return this->Data.GetIndex() == this->Data.GetIndexOf<SteadyStateDataType>();
}
VTKM_CONT bool IsUnsteadyState() const
{
return this->Data.GetIndex() == this->Data.GetIndexOf<UnsteadyStateDataType>();
}
VTKM_CONT vtkm::Id GetID() const { return this->Id; }
VTKM_CONT void SetCopySeedFlag(bool val) { this->CopySeedArray = val; }
@ -86,12 +115,33 @@ public:
DSIStuff<ParticleType>& stuff);
template <typename ArrayType>
VTKM_CONT void GetVelocityField(
VTKM_CONT void GetSteadyStateVelocityField(
vtkm::worklet::particleadvection::VelocityField<ArrayType>& velocityField) const
{
auto assoc = this->DataSet.GetField(this->FieldName).GetAssociation();
VTKM_ASSERT(this->Data.GetIndex() == this->Data.GetIndexOf<SteadyStateDataType>());
const auto& data = this->Data.Get<SteadyStateDataType>();
this->GetVelocityField(data, velocityField);
}
template <typename ArrayType>
VTKM_CONT void GetUnsteadyStateVelocityField(
vtkm::worklet::particleadvection::VelocityField<ArrayType>& velocityField1,
vtkm::worklet::particleadvection::VelocityField<ArrayType>& velocityField2) const
{
VTKM_ASSERT(this->Data.GetIndex() == this->Data.GetIndexOf<UnsteadyStateDataType>());
const auto& data = this->Data.Get<UnsteadyStateDataType>();
this->GetVelocityField(data.DataSet1, velocityField1);
this->GetVelocityField(data.DataSet2, velocityField2);
}
template <typename ArrayType>
VTKM_CONT void GetVelocityField(
const vtkm::cont::DataSet& ds,
vtkm::worklet::particleadvection::VelocityField<ArrayType>& velocityField) const
{
auto assoc = ds.GetField(this->FieldName).GetAssociation();
ArrayType arr;
vtkm::cont::ArrayCopyShallowIfPossible(this->DataSet.GetField(this->FieldName).GetData(), arr);
vtkm::cont::ArrayCopyShallowIfPossible(ds.GetField(this->FieldName).GetData(), arr);
velocityField = vtkm::worklet::particleadvection::VelocityField<ArrayType>(arr, assoc);
}
@ -101,14 +151,14 @@ public:
vtkm::worklet::particleadvection::ElectroMagneticField<ArrayType>& elecMagField) const
{
std::cout << "FIX ME: need fieldname2" << std::endl;
/*
ArrayType arr1, arr2;
vtkm::cont::ArrayCopyShallowIfPossible(this->DataSet.GetField(this->FieldName).GetData(), arr1);
vtkm::cont::ArrayCopyShallowIfPossible(this->DataSet.GetField(this->FieldName).GetData(),
arr2); //fieldname 2
vtkm::cont::ArrayCopyShallowIfPossible(this->DataSet.GetField(this->FieldName).GetData(), arr2); //fieldname 2
auto assoc = this->DataSet.GetField(this->FieldName).GetAssociation();
elecMagField =
vtkm::worklet::particleadvection::ElectroMagneticField<ArrayType>(arr1, arr2, assoc);
elecMagField = vtkm::worklet::particleadvection::ElectroMagneticField<ArrayType>(arr1, arr2, assoc);
*/
}
//template <typename ParticleType>
@ -118,12 +168,8 @@ public:
VTKM_CONT void ClassifyParticles(const vtkm::cont::ArrayHandle<ParticleType>& particles,
DSIStuff<ParticleType>& stuff) const;
template <typename ParticleType, template <typename> class ResultType>
VTKM_CONT vtkm::cont::DataSet ResultToDataSet() const;
vtkmdiy::mpi::communicator Comm = vtkm::cont::EnvironmentTracker::GetCommunicator();
vtkm::Id Rank;
vtkm::cont::DataSet DataSet;
vtkm::Id Id;
std::string FieldName;
vtkm::filter::particleadvection::IntegrationSolverType SolverType;
@ -132,15 +178,48 @@ public:
vtkm::filter::particleadvection::UNKNOWN_TYPE;
bool CopySeedArray = false;
/*
struct SteadyStateDataType
{
SteadyStateDataType(const vtkm::cont::DataSet& ds) : DataSet(ds) {}
vtkm::cont::DataSet DataSet;
};
*/
using ResultType =
using SteadyStateDataType = vtkm::cont::DataSet;
struct UnsteadyStateDataType
{
UnsteadyStateDataType(const vtkm::cont::DataSet& ds1,
const vtkm::cont::DataSet& ds2,
vtkm::FloatDefault t1,
vtkm::FloatDefault t2)
: DataSet1(ds1)
, DataSet2(ds2)
, Time1(t1)
, Time2(t2)
{
}
vtkm::cont::DataSet DataSet1;
vtkm::cont::DataSet DataSet2;
vtkm::FloatDefault Time1;
vtkm::FloatDefault Time2;
};
using DSType = vtkm::cont::internal::Variant<SteadyStateDataType, UnsteadyStateDataType>;
vtkm::cont::internal::Variant<SteadyStateDataType, UnsteadyStateDataType> Data;
using RType =
vtkm::cont::internal::Variant<vtkm::worklet::ParticleAdvectionResult<vtkm::Particle>,
vtkm::worklet::ParticleAdvectionResult<vtkm::ChargedParticle>,
vtkm::worklet::StreamlineResult<vtkm::Particle>,
vtkm::worklet::StreamlineResult<vtkm::ChargedParticle>>;
std::vector<ResultType> Results;
std::vector<RType> Results;
};
}
}
}

208
vtkm/filter/particleadvection/DSI.hxx
View File

@ -17,6 +17,88 @@ namespace filter
{
namespace particleadvection
{
namespace internal
{
template <typename GridEvalType,
typename WorkletType,
template <typename>
class ResultType,
typename ParticleType,
template <typename>
class IntType>
class AdvectHelper;
using ArrayType = vtkm::cont::ArrayHandle<vtkm::Vec3f>;
using FieldType = vtkm::worklet::particleadvection::VelocityField<ArrayType>;
using GridEvType = vtkm::worklet::particleadvection::GridEvaluator<FieldType>;
using TempGridEvType = vtkm::worklet::particleadvection::TemporalGridEvaluator<FieldType>;
//Steady state
template <typename WorkletType,
template <typename>
class ResultType,
typename ParticleType,
template <typename>
class IntType>
class AdvectHelper<GridEvType, WorkletType, ResultType, ParticleType, IntType>
{
public:
static void Advect(const FieldType& velField,
const vtkm::cont::DataSet& ds,
//const vtkm::filter::particleadvection::DSI::SteadyStateDataType& data,
vtkm::cont::ArrayHandle<ParticleType>& seedArray,
vtkm::FloatDefault stepSize,
vtkm::Id maxSteps,
ResultType<ParticleType>& result)
{
std::cout << "ADVECT HELPER!!!!!!!" << std::endl;
using StepperType = vtkm::worklet::particleadvection::Stepper<IntType<GridEvType>, GridEvType>;
GridEvType eval(ds, velField);
StepperType stepper(eval, stepSize);
WorkletType worklet;
result = worklet.Run(stepper, seedArray, maxSteps);
}
};
//unSteady state
template <typename WorkletType,
template <typename>
class ResultType,
typename ParticleType,
template <typename>
class IntType>
class AdvectHelper<TempGridEvType, WorkletType, ResultType, ParticleType, IntType>
{
public:
static void Advect(const FieldType& velField1,
const FieldType& velField2,
const vtkm::filter::particleadvection::DSI::UnsteadyStateDataType& data,
vtkm::cont::ArrayHandle<ParticleType>& seedArray,
vtkm::FloatDefault stepSize,
vtkm::Id maxSteps,
ResultType<ParticleType>& result)
{
std::cout << "TEMP ADVECT HELPER!!!!!!!" << std::endl;
using StepperType =
vtkm::worklet::particleadvection::Stepper<IntType<TempGridEvType>, TempGridEvType>;
TempGridEvType eval(data.DataSet1, data.Time1, velField1, data.DataSet2, data.Time2, velField2);
StepperType stepper(eval, stepSize);
WorkletType worklet;
result = worklet.Run(stepper, seedArray, maxSteps);
}
};
}
#if 0
namespace internal
{
@ -79,13 +161,13 @@ DSI::Meow(const char* func, const int& lineNum) const
std::cout << " ******************************";
std::cout << func << " " << lineNum << " ";
using RType = vtkm::worklet::ParticleAdvectionResult<vtkm::Particle>;
using ResType = vtkm::worklet::ParticleAdvectionResult<vtkm::Particle>;
//const auto& R0 = this->Results[0]->Get<RType>();
//const auto& P0 = R0.Particles.ReadPortal().Get(0);
//RType* r = static_cast<RType*>(this->Results[0]);
//const auto& P0 = r->Particles.ReadPortal().Get(0);
const auto& P0 = this->Results[0].Get<RType>().Particles.ReadPortal().Get(0);
const auto& P0 = this->Results[0].Get<ResType>().Particles.ReadPortal().Get(0);
std::cout << " PT0= " << P0.Pos << std::endl;
//std::cout<<" PT0= "<<P0.Pos<<std::endl;
@ -111,12 +193,12 @@ VTKM_CONT bool DSI::GetOutput(vtkm::cont::DataSet& ds) const
if (this->ResType == PARTICLE_ADVECT_TYPE)
{
using RType = vtkm::worklet::ParticleAdvectionResult<ParticleType>;
using ResType = vtkm::worklet::ParticleAdvectionResult<ParticleType>;
std::vector<vtkm::cont::ArrayHandle<ParticleType>> allParticles;
allParticles.reserve(nResults);
for (const auto& vres : this->Results)
allParticles.push_back(vres.Get<RType>().Particles);
allParticles.push_back(vres.Get<ResType>().Particles);
vtkm::cont::ArrayHandle<vtkm::Vec3f> pts;
vtkm::cont::ParticleArrayCopy(allParticles, pts);
@ -140,12 +222,12 @@ VTKM_CONT bool DSI::GetOutput(vtkm::cont::DataSet& ds) const
}
else if (this->ResType == STREAMLINE_TYPE)
{
using RType = vtkm::worklet::StreamlineResult<ParticleType>;
using ResType = vtkm::worklet::StreamlineResult<ParticleType>;
//Easy case with one result.
if (nResults == 1)
{
const auto& res = this->Results[0].Get<RType>();
const auto& res = this->Results[0].Get<ResType>();
ds.AddCoordinateSystem(vtkm::cont::CoordinateSystem("coordinates", res.Positions));
ds.SetCellSet(res.PolyLines);
}
@ -155,7 +237,7 @@ VTKM_CONT bool DSI::GetOutput(vtkm::cont::DataSet& ds) const
vtkm::Id totalNumCells = 0, totalNumPts = 0;
for (std::size_t i = 0; i < nResults; i++)
{
const auto& res = this->Results[i].Get<RType>();
const auto& res = this->Results[i].Get<ResType>();
if (i == 0)
posOffsets[i] = 0;
else
@ -170,7 +252,7 @@ VTKM_CONT bool DSI::GetOutput(vtkm::cont::DataSet& ds) const
appendPts.Allocate(totalNumPts);
for (std::size_t i = 0; i < nResults; i++)
{
const auto& res = this->Results[i].Get<RType>();
const auto& res = this->Results[i].Get<ResType>();
// copy all values into appendPts starting at offset.
vtkm::cont::Algorithm::CopySubRange(
res.Positions, 0, res.Positions.GetNumberOfValues(), appendPts, posOffsets[i]);
@ -182,7 +264,7 @@ VTKM_CONT bool DSI::GetOutput(vtkm::cont::DataSet& ds) const
std::size_t off = 0;
for (std::size_t i = 0; i < nResults; i++)
{
const auto& res = this->Results[i].Get<RType>();
const auto& res = this->Results[i].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));
@ -382,18 +464,117 @@ VTKM_CONT void DSI::Advect(std::vector<ParticleType>& v,
std::cout << "DSI::Advect() " << v.size() << std::endl;
//Assume all RK4.
if (this->VecFieldType == VELOCITY_FIELD_TYPE)
{
using FieldType = vtkm::worklet::particleadvection::VelocityField<ArrayType>;
if (this->IsSteadyState())
{
const auto& data = this->Data.Get<SteadyStateDataType>();
using GridEvType = vtkm::worklet::particleadvection::GridEvaluator<FieldType>;
FieldType velField;
this->GetSteadyStateVelocityField(velField);
if (this->ResType == PARTICLE_ADVECT_TYPE)
{
using AHType = internal::AdvectHelper<GridEvType,
vtkm::worklet::ParticleAdvection,
vtkm::worklet::ParticleAdvectionResult,
ParticleType,
vtkm::worklet::particleadvection::RK4Integrator>;
vtkm::worklet::ParticleAdvectionResult<ParticleType> result;
AHType::Advect(velField, data, seedArray, stepSize, maxSteps, result);
this->UpdateResult(result, stuff);
}
else
{
using AHType = internal::AdvectHelper<GridEvType,
vtkm::worklet::Streamline,
vtkm::worklet::StreamlineResult,
ParticleType,
vtkm::worklet::particleadvection::RK4Integrator>;
vtkm::worklet::StreamlineResult<ParticleType> result;
AHType::Advect(velField, data, seedArray, stepSize, maxSteps, result);
this->UpdateResult(result, stuff);
}
}
else if (this->IsUnsteadyState())
{
const auto& data = this->Data.Get<UnsteadyStateDataType>();
using GridEvType = vtkm::worklet::particleadvection::TemporalGridEvaluator<FieldType>;
FieldType velField1, velField2;
this->GetUnsteadyStateVelocityField(velField1, velField2);
if (this->ResType == PARTICLE_ADVECT_TYPE)
{
using AHType = internal::AdvectHelper<GridEvType,
vtkm::worklet::ParticleAdvection,
vtkm::worklet::ParticleAdvectionResult,
ParticleType,
vtkm::worklet::particleadvection::RK4Integrator>;
vtkm::worklet::ParticleAdvectionResult<ParticleType> result;
AHType::Advect(velField1, velField2, data, seedArray, stepSize, maxSteps, result);
this->UpdateResult(result, stuff);
}
else
{
using AHType = internal::AdvectHelper<GridEvType,
vtkm::worklet::Streamline,
vtkm::worklet::StreamlineResult,
ParticleType,
vtkm::worklet::particleadvection::RK4Integrator>;
vtkm::worklet::StreamlineResult<ParticleType> result;
AHType::Advect(velField1, velField2, data, seedArray, stepSize, maxSteps, result);
this->UpdateResult(result, stuff);
}
}
else
throw vtkm::cont::ErrorFilterExecution("Unsupported Data Type in DSI");
}
#if 0
if (this->SolverType == IntegrationSolverType::RK4_TYPE)
{
if (this->VecFieldType == VELOCITY_FIELD_TYPE) //vtkm::Particle, VelocityField
{
using FieldType = vtkm::worklet::particleadvection::VelocityField<ArrayType>;
using GridEvType = vtkm::worklet::particleadvection::GridEvaluator<FieldType>;
using RK4_Type = vtkm::worklet::particleadvection::RK4Integrator<GridEvType>;
using StepperType = vtkm::worklet::particleadvection::Stepper<RK4_Type, GridEvType>;
//using GridEvType = vtkm::worklet::particleadvection::GridEvaluator<FieldType>;
using GridEvType = vtkm::worklet::particleadvection::TemporalGridEvaluator<FieldType>;
FieldType velField;
this->GetVelocityField(velField);
if (this->ResType == PARTICLE_ADVECT_TYPE)
{
using AHType = internal::AdvectHelper<GridEvType,
vtkm::worklet::ParticleAdvection,
vtkm::worklet::ParticleAdvectionResult,
ParticleType,
vtkm::worklet::particleadvection::RK4Integrator>;
vtkm::worklet::ParticleAdvectionResult<ParticleType> result;
//AHType::Advect(velField, this->DataSet, seedArray, stepSize, maxSteps, result);
AHType::Advect(velField, velField, this->DataSet, this->DataSet, 0.0, 1.0, seedArray, stepSize, maxSteps, result);
this->UpdateResult(result, stuff);
}
else
{
using AHType = internal::AdvectHelper<GridEvType,
vtkm::worklet::Streamline,
vtkm::worklet::StreamlineResult,
ParticleType,
vtkm::worklet::particleadvection::RK4Integrator>;
vtkm::worklet::StreamlineResult<ParticleType> result;
//AHType::Advect(velField, this->DataSet, seedArray, stepSize, maxSteps, result);
AHType::Advect(velField, velField, this->DataSet, this->DataSet, 0.0, 1.0, seedArray, stepSize, maxSteps, result);
this->UpdateResult(result, stuff);
}
/*
GridEvType eval(this->DataSet, velField);
StepperType stepper(eval, stepSize);
@ -411,6 +592,7 @@ VTKM_CONT void DSI::Advect(std::vector<ParticleType>& v,
this->UpdateResult(r, stuff);
//Put results in unknown array??
}
*/
}
else if (this->VecFieldType == ELECTRO_MAGNETIC_FIELD_TYPE) //vtkm::ChargedParticle
{
@ -461,8 +643,10 @@ VTKM_CONT void DSI::Advect(std::vector<ParticleType>& v,
//Put results in unknown array??
}
}
#endif
}
}
}
}

18
vtkm/filter/particleadvection/PAV.h
View File

@ -22,11 +22,12 @@ namespace filter
namespace particleadvection
{
template <typename DSIType>
class PAV
{
public:
PAV(const vtkm::filter::particleadvection::BoundsMap& bm,
const std::vector<DSI*> blocks,
const std::vector<DSIType*> blocks,
const bool& useThreaded,
const vtkm::filter::particleadvection::ParticleAdvectionResultType& parType)
: Blocks(blocks)
@ -63,9 +64,8 @@ private:
if (this->ResultType ==
vtkm::filter::particleadvection::ParticleAdvectionResultType::PARTICLE_ADVECT_TYPE)
{
using AlgorithmType =
vtkm::filter::particleadvection::ABA<vtkm::worklet::ParticleAdvectionResult,
ParticleType>;
using AlgorithmType = vtkm::filter::particleadvection::
ABA<DSIType, vtkm::worklet::ParticleAdvectionResult, ParticleType>;
AlgorithmType algo(this->BoundsMap, this->Blocks);
algo.Execute(numSteps, stepSize, seeds);
@ -80,8 +80,8 @@ private:
}
else
{
using AlgorithmType =
vtkm::filter::particleadvection::ABA<vtkm::worklet::StreamlineResult, ParticleType>;
using AlgorithmType = vtkm::filter::particleadvection::
ABA<DSIType, vtkm::worklet::StreamlineResult, ParticleType>;
AlgorithmType algo(this->BoundsMap, this->Blocks);
algo.Execute(numSteps, stepSize, seeds);
@ -91,8 +91,8 @@ private:
else
{
std::cout << "Change me to threaded ABA" << std::endl;
using AlgorithmType =
vtkm::filter::particleadvection::ABA<vtkm::worklet::ParticleAdvectionResult, ParticleType>;
using AlgorithmType = vtkm::filter::particleadvection::
ABA<DSIType, vtkm::worklet::ParticleAdvectionResult, ParticleType>;
AlgorithmType algo(this->BoundsMap, this->Blocks);
algo.Execute(numSteps, stepSize, seeds);
@ -101,7 +101,7 @@ private:
}
std::vector<DSI*> Blocks;
std::vector<DSIType*> Blocks;
vtkm::filter::particleadvection::BoundsMap BoundsMap;
ParticleAdvectionResultType ResultType;
bool UseThreadedAlgorithm;

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

@ -68,7 +68,7 @@ void TestStreamline()
vtkm::Particle(vtkm::Vec3f(.2f, 2.0f, .2f), 1),
vtkm::Particle(vtkm::Vec3f(.2f, 3.0f, .2f), 2) });
vtkm::filter::Streamline streamline;
vtkm::filter::Streamline2 streamline;
streamline.SetStepSize(0.1f);
streamline.SetNumberOfSteps(20);
@ -89,6 +89,12 @@ void TestStreamline()
}
}
void DumpDS(const vtkm::cont::DataSet& ds, const std::string& fname)
{
vtkm::io::VTKDataSetWriter wrt(fname);
wrt.WriteDataSet(ds);
}
void TestPathline()
{
const vtkm::Id3 dims(5, 5, 5);
@ -127,7 +133,7 @@ void TestPathline()
vtkm::Id numExpectedPoints;
if (fType == 0)
{
vtkm::filter::Pathline filt;
vtkm::filter::Pathline2 filt;
filt.SetActiveField(var);
filt.SetStepSize(stepSize);
filt.SetNumberOfSteps(numSteps);
@ -137,10 +143,15 @@ void TestPathline()
filt.SetNextDataSet(ds2);
output = filt.Execute(ds1);
numExpectedPoints = 33;
/*
DumpDS(output, "pathline.vtk");
DumpDS(ds1, "ds1.vtk");
DumpDS(ds2, "ds2.vtk");
*/
}
else
{
vtkm::filter::PathParticle filt;
vtkm::filter::PathParticle2 filt;
filt.SetActiveField(var);
filt.SetStepSize(stepSize);
filt.SetNumberOfSteps(numSteps);
@ -226,7 +237,7 @@ void TestAMRStreamline(bool useSL)
if (useSL)
{
vtkm::filter::Streamline filter;
vtkm::filter::Streamline2 filter;
filter.SetStepSize(0.1f);
filter.SetNumberOfSteps(100000);
filter.SetSeeds(seedArray);
@ -375,18 +386,6 @@ void TestPartitionedDataSet(vtkm::Id num, bool useGhost, FilterType fType)
vtkm::Particle(vtkm::Vec3f(.2f, 2.0f, .2f), 1) });
vtkm::Id numSeeds = seedArray.GetNumberOfValues();
if (idx == 0)
{
for (int j = 0; j < pds.GetNumberOfPartitions(); j++)
{
auto ds_j = pds.GetPartition(j);
std::string nm = "ds_" + std::to_string(j) + ".vtk";
vtkm::io::VTKDataSetWriter wrt(nm);
wrt.WriteDataSet(ds_j);
}
}
if (fType == FilterType::STREAMLINE || fType == FilterType::PATHLINE)
{
vtkm::cont::PartitionedDataSet out;
@ -568,7 +567,7 @@ void TestStreamlineFile(const std::string& fname,
vtkm::cont::DataSet output;
if (useSL)
{
vtkm::filter::Streamline streamline;
vtkm::filter::Streamline2 streamline;
streamline.SetStepSize(stepSize);
streamline.SetNumberOfSteps(maxSteps);
streamline.SetSeeds(seedArray);
@ -612,7 +611,7 @@ void TestStreamlineFilters()
FilterType::PATH_PARTICLE };
//fTypes = {FilterType::PARTICLE_ADVECTION,FilterType::STREAMLINE};
//fTypes = {FilterType::PATH_PARTICLE};
//fTypes = {FilterType::PATHLINE};
for (int n = 1; n < 3; n++)
{
for (auto useGhost : flags)
@ -620,14 +619,14 @@ void TestStreamlineFilters()
TestPartitionedDataSet(n, useGhost, ft);
}
return;
TestStreamline();
TestPathline();
for (auto useSL : flags)
TestAMRStreamline(useSL);
return;
//Fusion test.
std::vector<vtkm::Vec3f> fusionPts, fusionEndPts;
fusionPts.push_back(vtkm::Vec3f(0.8f, 0.6f, 0.6f));