Filters for charged particle advection.

vtkm/filter/FilterParticleAdvection.h

@@ -25,7 +25,7 @@ namespace filter
 /// Takes as input a vector field and seed locations and advects the seeds
 /// through the flow field.
 
-template <class Derived>
+template <class Derived, typename ParticleType>
 class FilterParticleAdvection : public vtkm::filter::FilterDataSetWithField<Derived>
 {
 public:
@@ -34,6 +34,16 @@ public:
   VTKM_CONT
   FilterParticleAdvection();
 
+  VTKM_CONT
+  void SetActiveFields(
+    const std::string& name1,
+    const std::string& name2,
+    vtkm::cont::Field::Association association = vtkm::cont::Field::Association::ANY)
+  {
+    this->FilterDataSetWithField<Derived>::SetActiveField(name1, association);
+    this->ActiveFieldName2 = name2;
+  }
+
   VTKM_CONT
   void SetStepSize(vtkm::FloatDefault s) { this->StepSize = s; }
 
@@ -41,14 +51,14 @@ public:
   void SetNumberOfSteps(vtkm::Id n) { this->NumberOfSteps = n; }
 
   VTKM_CONT
-  void SetSeeds(const std::vector<vtkm::Particle>& seeds,
+  void SetSeeds(const std::vector<ParticleType>& seeds,
                 vtkm::CopyFlag copyFlag = vtkm::CopyFlag::On)
   {
     this->Seeds = vtkm::cont::make_ArrayHandle(seeds, copyFlag);
   }
 
   VTKM_CONT
-  void SetSeeds(vtkm::cont::ArrayHandle<vtkm::Particle>& seeds) { this->Seeds = seeds; }
+  void SetSeeds(vtkm::cont::ArrayHandle<ParticleType>& seeds) { this->Seeds = seeds; }
 
   VTKM_CONT
   bool GetUseThreadedAlgorithm() { return this->UseThreadedAlgorithm; }
@@ -72,13 +82,16 @@ protected:
   VTKM_CONT virtual void ValidateOptions() const;
 
   using DSIType = vtkm::filter::particleadvection::DataSetIntegrator;
+  using DSIType2 = vtkm::filter::particleadvection::DataSetIntegrator2;
+
   VTKM_CONT std::vector<DSIType> CreateDataSetIntegrators(
     const vtkm::cont::PartitionedDataSet& input,
     const vtkm::filter::particleadvection::BoundsMap& boundsMap) const;
 
+  std::string ActiveFieldName2;
   vtkm::Id NumberOfSteps;
   vtkm::FloatDefault StepSize;
-  vtkm::cont::ArrayHandle<vtkm::Particle> Seeds;
+  vtkm::cont::ArrayHandle<ParticleType> Seeds;
   bool UseThreadedAlgorithm;
 
 private:

vtkm/filter/FilterParticleAdvection.hxx

@@ -18,8 +18,8 @@ namespace filter
 {
 
 //-----------------------------------------------------------------------------
-template <typename Derived>
-inline VTKM_CONT FilterParticleAdvection<Derived>::FilterParticleAdvection()
+template <typename Derived, typename ParticleType>
+inline VTKM_CONT FilterParticleAdvection<Derived, ParticleType>::FilterParticleAdvection()
   : vtkm::filter::FilterDataSetWithField<Derived>()
   , NumberOfSteps(0)
   , StepSize(0)
@@ -27,8 +27,8 @@ inline VTKM_CONT FilterParticleAdvection<Derived>::FilterParticleAdvection()
 {
 }
 
-template <typename Derived>
-void FilterParticleAdvection<Derived>::ValidateOptions() const
+template <typename Derived, typename ParticleType>
+void FilterParticleAdvection<Derived, ParticleType>::ValidateOptions() const
 {
   if (this->GetUseCoordinateSystemAsField())
     throw vtkm::cont::ErrorFilterExecution("Coordinate system as field not supported");
@@ -40,9 +40,9 @@ void FilterParticleAdvection<Derived>::ValidateOptions() const
     throw vtkm::cont::ErrorFilterExecution("Step size not specified.");
 }
 
-template <typename Derived>
+template <typename Derived, typename ParticleType>
 std::vector<vtkm::filter::particleadvection::DataSetIntegrator>
-FilterParticleAdvection<Derived>::CreateDataSetIntegrators(
+FilterParticleAdvection<Derived, ParticleType>::CreateDataSetIntegrators(
   const vtkm::cont::PartitionedDataSet& input,
   const vtkm::filter::particleadvection::BoundsMap& boundsMap) const
 {
@@ -66,9 +66,10 @@ FilterParticleAdvection<Derived>::CreateDataSetIntegrators(
 }
 
 //-----------------------------------------------------------------------------
-template <typename Derived>
+template <typename Derived, typename ParticleType>
 template <typename DerivedPolicy>
-inline VTKM_CONT vtkm::cont::DataSet FilterParticleAdvection<Derived>::PrepareForExecution(
+inline VTKM_CONT vtkm::cont::DataSet
+FilterParticleAdvection<Derived, ParticleType>::PrepareForExecution(
   const vtkm::cont::DataSet& input,
   vtkm::filter::PolicyBase<DerivedPolicy> policy)
 {

vtkm/filter/FilterTemporalParticleAdvection.h

@@ -25,8 +25,9 @@ namespace filter
 /// Takes as input a vector field and seed locations and advects the seeds
 /// through the flow field.
 
-template <class Derived>
-class FilterTemporalParticleAdvection : public vtkm::filter::FilterParticleAdvection<Derived>
+template <class Derived, typename ParticleType>
+class FilterTemporalParticleAdvection
+  : public vtkm::filter::FilterParticleAdvection<Derived, ParticleType>
 {
 public:
   VTKM_CONT
@@ -52,7 +53,7 @@ public:
 
 protected:
   VTKM_CONT void ValidateOptions(const vtkm::cont::PartitionedDataSet& input) const;
-  using vtkm::filter::FilterParticleAdvection<Derived>::ValidateOptions;
+  using vtkm::filter::FilterParticleAdvection<Derived, ParticleType>::ValidateOptions;
 
   using DSIType = vtkm::filter::particleadvection::TemporalDataSetIntegrator;
   VTKM_CONT std::vector<DSIType> CreateDataSetIntegrators(

vtkm/filter/FilterTemporalParticleAdvection.hxx

@@ -18,19 +18,20 @@ namespace filter
 {
 
 //-----------------------------------------------------------------------------
-template <typename Derived>
-inline VTKM_CONT FilterTemporalParticleAdvection<Derived>::FilterTemporalParticleAdvection()
-  : vtkm::filter::FilterParticleAdvection<Derived>()
+template <class Derived, typename ParticleType>
+inline VTKM_CONT
+FilterTemporalParticleAdvection<Derived, ParticleType>::FilterTemporalParticleAdvection()
+  : vtkm::filter::FilterParticleAdvection<Derived, ParticleType>()
   , PreviousTime(0)
   , NextTime(0)
 {
 }
 
-template <typename Derived>
-void FilterTemporalParticleAdvection<Derived>::ValidateOptions(
+template <typename Derived, typename ParticleType>
+void FilterTemporalParticleAdvection<Derived, ParticleType>::ValidateOptions(
   const vtkm::cont::PartitionedDataSet& input) const
 {
-  this->vtkm::filter::FilterParticleAdvection<Derived>::ValidateOptions();
+  this->vtkm::filter::FilterParticleAdvection<Derived, ParticleType>::ValidateOptions();
 
   if (this->NextDataSet.GetNumberOfPartitions() != input.GetNumberOfPartitions())
     throw vtkm::cont::ErrorFilterExecution("Number of partitions do not match");
@@ -38,9 +39,9 @@ void FilterTemporalParticleAdvection<Derived>::ValidateOptions(
     throw vtkm::cont::ErrorFilterExecution("Previous time must be less than Next time.");
 }
 
-template <typename Derived>
+template <typename Derived, typename ParticleType>
 std::vector<vtkm::filter::particleadvection::TemporalDataSetIntegrator>
-FilterTemporalParticleAdvection<Derived>::CreateDataSetIntegrators(
+FilterTemporalParticleAdvection<Derived, ParticleType>::CreateDataSetIntegrators(
   const vtkm::cont::PartitionedDataSet& input,
   const vtkm::filter::particleadvection::BoundsMap& boundsMap) const
 {
@@ -65,9 +66,10 @@ FilterTemporalParticleAdvection<Derived>::CreateDataSetIntegrators(
 }
 
 //-----------------------------------------------------------------------------
-template <typename Derived>
+template <typename Derived, typename ParticleType>
 template <typename DerivedPolicy>
-inline VTKM_CONT vtkm::cont::DataSet FilterTemporalParticleAdvection<Derived>::PrepareForExecution(
+inline VTKM_CONT vtkm::cont::DataSet
+FilterTemporalParticleAdvection<Derived, ParticleType>::PrepareForExecution(
   const vtkm::cont::DataSet& input,
   vtkm::filter::PolicyBase<DerivedPolicy> policy)
 {

vtkm/filter/ParticleAdvection.h

@@ -23,7 +23,8 @@ namespace filter
 /// Takes as input a vector field and seed locations and generates the
 /// end points for each seed through the vector field.
 
-class ParticleAdvection : public vtkm::filter::FilterParticleAdvection<ParticleAdvection>
+class ParticleAdvection
+  : public vtkm::filter::FilterParticleAdvection<ParticleAdvection, vtkm::Particle>
 {
 public:
   VTKM_CONT ParticleAdvection();

vtkm/filter/ParticleAdvection.hxx

@@ -23,7 +23,7 @@ namespace filter
 
 //-----------------------------------------------------------------------------
 inline VTKM_CONT ParticleAdvection::ParticleAdvection()
-  : vtkm::filter::FilterParticleAdvection<ParticleAdvection>()
+  : vtkm::filter::FilterParticleAdvection<ParticleAdvection, vtkm::Particle>()
 {
 }
 
@@ -34,18 +34,42 @@ inline VTKM_CONT vtkm::cont::PartitionedDataSet ParticleAdvection::PrepareForExe
   const vtkm::filter::PolicyBase<DerivedPolicy>&)
 {
   using AlgorithmType = vtkm::filter::particleadvection::ParticleAdvectionAlgorithm;
+  using AlgorithmType2 = vtkm::filter::particleadvection::ParticleAdvectionAlgorithm2;
   using ThreadedAlgorithmType = vtkm::filter::particleadvection::ParticleAdvectionThreadedAlgorithm;
 
   this->ValidateOptions();
   vtkm::filter::particleadvection::BoundsMap boundsMap(input);
-  auto dsi = this->CreateDataSetIntegrators(input, boundsMap);
 
-  if (this->GetUseThreadedAlgorithm())
-    return vtkm::filter::particleadvection::RunAlgo<DSIType, ThreadedAlgorithmType>(
-      boundsMap, dsi, this->NumberOfSteps, this->StepSize, this->Seeds);
+  if (this->ActiveFieldName2.empty())
+  {
+    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);
+  }
   else
-    return vtkm::filter::particleadvection::RunAlgo<DSIType, AlgorithmType>(
+  {
+    std::vector<vtkm::filter::particleadvection::DataSetIntegrator2> dsi;
+    std::string field1 = this->GetActiveFieldName();
+    std::string field2 = this->ActiveFieldName2;
+
+    for (vtkm::Id i = 0; i < input.GetNumberOfPartitions(); i++)
+    {
+      vtkm::Id blockId = boundsMap.GetLocalBlockId(i);
+      auto ds = input.GetPartition(i);
+      if (!ds.HasPointField(field1) || !ds.HasPointField(field2))
+        throw vtkm::cont::ErrorFilterExecution("Unsupported field assocation");
+      dsi.push_back(DSIType2(ds, blockId, field1, field2));
+    }
+
+    std::cout << "Electromagnietc RunAlgo" << std::endl;
+    return vtkm::filter::particleadvection::RunAlgo<DSIType2, AlgorithmType2>(
       boundsMap, dsi, this->NumberOfSteps, this->StepSize, this->Seeds);
+  }
 }
 
 }

vtkm/filter/PathParticle.h

@@ -21,7 +21,8 @@ namespace filter
 
 /// Takes as input a vector field and seed locations and generates the
 /// paths taken by the seeds through the vector field.
-class PathParticle : public vtkm::filter::FilterTemporalParticleAdvection<PathParticle>
+class PathParticle
+  : public vtkm::filter::FilterTemporalParticleAdvection<PathParticle, vtkm::Particle>
 {
 public:
   VTKM_CONT

vtkm/filter/PathParticle.hxx

@@ -24,7 +24,7 @@ namespace filter
 
 //-----------------------------------------------------------------------------
 inline VTKM_CONT PathParticle::PathParticle()
-  : vtkm::filter::FilterTemporalParticleAdvection<PathParticle>()
+  : vtkm::filter::FilterTemporalParticleAdvection<PathParticle, vtkm::Particle>()
 {
 }
 

vtkm/filter/Pathline.h

@@ -21,7 +21,7 @@ namespace filter
 
 /// Takes as input a vector field and seed locations and generates the
 /// paths taken by the seeds through the vector field.
-class Pathline : public vtkm::filter::FilterTemporalParticleAdvection<Pathline>
+class Pathline : public vtkm::filter::FilterTemporalParticleAdvection<Pathline, vtkm::Particle>
 {
 public:
   VTKM_CONT

vtkm/filter/Pathline.hxx

@@ -24,7 +24,7 @@ namespace filter
 
 //-----------------------------------------------------------------------------
 inline VTKM_CONT Pathline::Pathline()
-  : vtkm::filter::FilterTemporalParticleAdvection<Pathline>()
+  : vtkm::filter::FilterTemporalParticleAdvection<Pathline, vtkm::Particle>()
 {
 }
 

vtkm/filter/Streamline.h

@@ -21,7 +21,7 @@ namespace filter
 
 /// Takes as input a vector field and seed locations and generates the
 /// paths taken by the seeds through the vector field.
-class Streamline : public vtkm::filter::FilterParticleAdvection<Streamline>
+class Streamline : public vtkm::filter::FilterParticleAdvection<Streamline, vtkm::Particle>
 {
 public:
   VTKM_CONT

vtkm/filter/Streamline.hxx

@@ -23,7 +23,7 @@ namespace filter
 
 //-----------------------------------------------------------------------------
 inline VTKM_CONT Streamline::Streamline()
-  : vtkm::filter::FilterParticleAdvection<Streamline>()
+  : vtkm::filter::FilterParticleAdvection<Streamline, vtkm::Particle>()
 {
 }
 

vtkm/filter/particleadvection/DataSetIntegrator.h

@@ -92,6 +92,7 @@ public:
   DataSetIntegrator(const vtkm::cont::DataSet& ds, vtkm::Id id, const std::string& fieldNm)
     : DataSetIntegratorBase<vtkm::worklet::particleadvection::GridEvaluator<
         vtkm::worklet::particleadvection::VelocityField<FieldHandleType>>>(false, id)
+  //    : DataSetIntegratorBase<vtkm::worklet::particleadvection::GridEvaluator<FieldType>>(false, id)
   {
     using Association = vtkm::cont::Field::Association;
     using FieldType = vtkm::worklet::particleadvection::VelocityField<FieldHandleType>;
@@ -103,6 +104,30 @@ public:
   }
 };
 
+class VTKM_ALWAYS_EXPORT DataSetIntegrator2
+  : public DataSetIntegratorBase<vtkm::worklet::particleadvection::GridEvaluator<
+      vtkm::worklet::particleadvection::ElectroMagneticField<vtkm::cont::ArrayHandle<vtkm::Vec3f>>>>
+{
+public:
+  DataSetIntegrator2(const vtkm::cont::DataSet& ds,
+                     vtkm::Id id,
+                     const std::string& fieldNm1,
+                     const std::string& fieldNm2)
+    : DataSetIntegratorBase<vtkm::worklet::particleadvection::GridEvaluator<
+        vtkm::worklet::particleadvection::ElectroMagneticField<FieldHandleType>>>(false, id)
+  {
+    using Association = vtkm::cont::Field::Association;
+    using FieldType = vtkm::worklet::particleadvection::ElectroMagneticField<FieldHandleType>;
+    using EvalType = vtkm::worklet::particleadvection::GridEvaluator<FieldType>;
+    Association association = ds.GetField(fieldNm1).GetAssociation();
+    auto fieldArray1 = this->GetFieldHandle(ds, fieldNm1);
+    auto fieldArray2 = this->GetFieldHandle(ds, fieldNm2);
+
+    FieldType field(fieldArray1, fieldArray2, association);
+    this->Eval = std::shared_ptr<EvalType>(new EvalType(ds, field));
+  }
+};
+
 class VTKM_ALWAYS_EXPORT TemporalDataSetIntegrator
   : public DataSetIntegratorBase<vtkm::worklet::particleadvection::TemporalGridEvaluator<
       vtkm::worklet::particleadvection::VelocityField<vtkm::cont::ArrayHandle<vtkm::Vec3f>>>>

vtkm/filter/particleadvection/DataSetIntegrator.hxx

@@ -20,6 +20,8 @@ namespace particleadvection
 
 using GridEvalType = vtkm::worklet::particleadvection::GridEvaluator<
   vtkm::worklet::particleadvection::VelocityField<vtkm::cont::ArrayHandle<vtkm::Vec3f>>>;
+using GridEvalType2 = vtkm::worklet::particleadvection::GridEvaluator<
+  vtkm::worklet::particleadvection::ElectroMagneticField<vtkm::cont::ArrayHandle<vtkm::Vec3f>>>;
 using TemporalGridEvalType = vtkm::worklet::particleadvection::TemporalGridEvaluator<
   vtkm::worklet::particleadvection::VelocityField<vtkm::cont::ArrayHandle<vtkm::Vec3f>>>;
 using RK4Type = vtkm::worklet::particleadvection::RK4Integrator<GridEvalType>;
@@ -42,6 +44,20 @@ inline void DataSetIntegratorBase<GridEvalType>::DoAdvect(
   result = Worklet.Run(stepper, seeds, maxSteps);
 }
 
+//-----
+// Specialization for ParticleAdvection worklet w/ electromagnetic field
+template <>
+template <>
+inline void DataSetIntegratorBase<GridEvalType2>::DoAdvect(
+  vtkm::cont::ArrayHandle<vtkm::Particle>& seeds,
+  const Stepper& stepper,
+  vtkm::Id maxSteps,
+  vtkm::worklet::ParticleAdvectionResult<vtkm::Particle>& result) const
+{
+  vtkm::worklet::ParticleAdvection Worklet;
+  result = Worklet.Run(stepper, seeds, maxSteps);
+}
+
 //-----
 // Specialization for Streamline worklet
 template <>

vtkm/filter/particleadvection/ParticleAdvectionAlgorithm.h

@@ -23,6 +23,7 @@ namespace particleadvection
 {
 
 using DSIType = vtkm::filter::particleadvection::DataSetIntegrator;
+using DSIType2 = vtkm::filter::particleadvection::DataSetIntegrator2;
 using TDSIType = vtkm::filter::particleadvection::TemporalDataSetIntegrator;
 
 //-------------------------------------------------------------------------------------------
@@ -41,6 +42,20 @@ public:
   }
 };
 
+//Particle advection
+class VTKM_ALWAYS_EXPORT ParticleAdvectionAlgorithm2
+  : public AdvectorBaseAlgorithm<DSIType2, vtkm::worklet::ParticleAdvectionResult<vtkm::Particle>>
+{
+public:
+  ParticleAdvectionAlgorithm2(const vtkm::filter::particleadvection::BoundsMap& bm,
+                              const std::vector<DSIType2>& blocks)
+    : AdvectorBaseAlgorithm<DSIType2, vtkm::worklet::ParticleAdvectionResult<vtkm::Particle>>(
+        bm,
+        blocks)
+  {
+  }
+};
+
 //Threaded particle advection
 class VTKM_ALWAYS_EXPORT ParticleAdvectionThreadedAlgorithm
   : public AdvectorBaseThreadedAlgorithm<DSIType,

vtkm/filter/testing/UnitTestStreamlineFilter.cxx

@@ -82,6 +82,44 @@ void TestStreamline()
     vtkm::cont::DynamicCellSet dcells = output.GetCellSet();
     VTKM_TEST_ASSERT(dcells.GetNumberOfCells() == 3, "Wrong number of cells");
   }
+
+  //Test ElectomagneticField advection
+  for (auto& ds : dataSets)
+  {
+    auto vecField = CreateConstantVectorField(ds.GetNumberOfPoints(), vecX);
+    ds.AddPointField(fieldName, vecField);
+    std::vector<vtkm::ChargedParticle> pts;
+    vtkm::FloatDefault mass = 1.0f, charge = 1.0f, weight = 1.0;
+    vtkm::Vec3f momentum(.1, .1, .1);
+    pts.push_back(
+      vtkm::ChargedParticle(vtkm::Vec3f(.2f, 1.0f, .2f), 0, mass, charge, weight, momentum));
+    pts.push_back(
+      vtkm::ChargedParticle(vtkm::Vec3f(.2f, 2.0f, .2f), 0, mass, charge, weight, momentum));
+    pts.push_back(
+      vtkm::ChargedParticle(vtkm::Vec3f(.2f, 3.0f, .2f), 0, mass, charge, weight, momentum));
+    vtkm::cont::ArrayHandle<vtkm::ChargedParticle> seedArray =
+      vtkm::cont::make_ArrayHandle(pts, vtkm::CopyFlag::On);
+
+    /*
+    vtkm::filter::ParticleAdvection filter;
+    filter.SetStepSize(0.1f);
+    filter.SetNumberOfSteps(20);
+    filter.SetSeeds(seedArray);
+
+    filter.SetActiveFields(fieldName, fieldName);
+    auto output = filter.Execute(ds);
+
+    //Validate the result is correct.
+    VTKM_TEST_ASSERT(output.GetNumberOfCoordinateSystems() == 1,
+                     "Wrong number of coordinate systems in the output dataset");
+
+    vtkm::cont::CoordinateSystem coords = output.GetCoordinateSystem();
+    VTKM_TEST_ASSERT(coords.GetNumberOfPoints() == 63, "Wrong number of coordinates");
+
+    vtkm::cont::DynamicCellSet dcells = output.GetCellSet();
+    VTKM_TEST_ASSERT(dcells.GetNumberOfCells() == 3, "Wrong number of cells");
+    */
+  }
 }
 
 void TestPathline()