Use types for different options in PA filters.

vtkm/VecTraits.h

@@ -498,6 +498,7 @@ struct VTKM_NEVER_EXPORT VecTraits<vtkm::VecCConst<T>>
 
 namespace internal
 {
+
 /// Used for overriding VecTraits for basic scalar types.
 ///
 template <typename ScalarType>
@@ -539,6 +540,44 @@ struct VTKM_NEVER_EXPORT VecTraitsBasic
     dest[0] = src;
   }
 };
+
+namespace detail
+{
+
+template <typename T, typename = vtkm::HasVecTraits<T>>
+struct VTKM_NEVER_EXPORT SafeVecTraitsImpl;
+
+template <typename T>
+struct VTKM_NEVER_EXPORT SafeVecTraitsImpl<T, std::true_type> : vtkm::VecTraits<T>
+{
+};
+
+template <typename T>
+struct VTKM_NEVER_EXPORT SafeVecTraitsImpl<T, std::false_type> : vtkm::internal::VecTraitsBasic<T>
+{
+};
+
+} // namespace detail
+
+/// \brief A version of VecTraits that will be available for any type.
+///
+/// The `VecTraits` template is only defined for types that have a specific specialization
+/// for it. That means if you use `VecTraits` in a template, that template will likely
+/// fail to build for types that are not defined for `VecTraits`.
+///
+/// To use `VecTraits` in a class that should support all types, not just those with
+/// defined `VecTraits`, you can use this "safe" version. `SafeVecTraits` is the same as
+/// `VecTraits` if the latter is defined. If the `VecTraits` are not defined, then
+/// `SafeVecTraits` treats the type as a simple scalar value.
+///
+/// This template ensures that it will work reasonably well for all types. But be careful
+/// as if `VecTraits` is later defined, the template is likely to change.
+///
+template <typename T>
+struct VTKM_NEVER_EXPORT SafeVecTraits : detail::SafeVecTraitsImpl<T>
+{
+};
+
 } // namespace internal
 
 /// \brief VecTraits for Pair types
@@ -554,7 +593,7 @@ struct VTKM_NEVER_EXPORT VecTraits<vtkm::Pair<T, U>>
 {
 };
 
-} // anonymous namespace
+} // namespace vtkm
 
 #define VTKM_BASIC_TYPE_VECTOR(type)                                                     \
   namespace vtkm                                                                         \

vtkm/cont/ArrayExtractComponent.h

@@ -37,7 +37,7 @@ namespace internal
 // is defined rather than where it is resolved. This causes problems when extracting
 // components of, say, an ArrayHandleMultiplexer holding an ArrayHandleSOA.
 template <typename T, typename S>
-vtkm::cont::ArrayHandleStride<typename vtkm::VecTraits<T>::BaseComponentType>
+vtkm::cont::ArrayHandleStride<typename vtkm::internal::SafeVecTraits<T>::BaseComponentType>
 ArrayExtractComponentFallback(const vtkm::cont::ArrayHandle<T, S>& src,
                               vtkm::IdComponent componentIndex,
                               vtkm::CopyFlag allowCopy)
@@ -53,7 +53,7 @@ ArrayExtractComponentFallback(const vtkm::cont::ArrayHandle<T, S>& src,
                                      << vtkm::cont::TypeToString<vtkm::cont::ArrayHandle<T, S>>()
                                      << " requires an inefficient memory copy.");
 
-  using BaseComponentType = typename vtkm::VecTraits<T>::BaseComponentType;
+  using BaseComponentType = typename vtkm::internal::SafeVecTraits<T>::BaseComponentType;
   vtkm::Id numValues = src.GetNumberOfValues();
   vtkm::cont::ArrayHandleBasic<BaseComponentType> dest;
   dest.Allocate(numValues);
@@ -78,10 +78,10 @@ template <typename S>
 struct ArrayExtractComponentImpl : ArrayExtractComponentImplInefficient
 {
   template <typename T>
-  vtkm::cont::ArrayHandleStride<typename vtkm::VecTraits<T>::BaseComponentType> operator()(
+  vtkm::cont::ArrayHandleStride<typename vtkm::internal::SafeVecTraits<T>::BaseComponentType>
-    const vtkm::cont::ArrayHandle<T, S>& src,
+  operator()(const vtkm::cont::ArrayHandle<T, S>& src,
-    vtkm::IdComponent componentIndex,
+             vtkm::IdComponent componentIndex,
-    vtkm::CopyFlag allowCopy) const
+             vtkm::CopyFlag allowCopy) const
   {
     // This is the slow "default" implementation. ArrayHandle implementations should provide
     // more efficient overloads where applicable.
@@ -93,13 +93,15 @@ template <>
 struct ArrayExtractComponentImpl<vtkm::cont::StorageTagStride>
 {
   template <typename T>
-  vtkm::cont::ArrayHandleStride<typename vtkm::VecTraits<T>::BaseComponentType> operator()(
+  vtkm::cont::ArrayHandleStride<typename vtkm::internal::SafeVecTraits<T>::BaseComponentType>
-    const vtkm::cont::ArrayHandle<T, vtkm::cont::StorageTagStride>& src,
+  operator()(const vtkm::cont::ArrayHandle<T, vtkm::cont::StorageTagStride>& src,
-    vtkm::IdComponent componentIndex,
+             vtkm::IdComponent componentIndex,
-    vtkm::CopyFlag allowCopy) const
+             vtkm::CopyFlag allowCopy) const
   {
-    return this->DoExtract(
+    return this->DoExtract(src,
-      src, componentIndex, allowCopy, typename vtkm::VecTraits<T>::HasMultipleComponents{});
+                           componentIndex,
+                           allowCopy,
+                           typename vtkm::internal::SafeVecTraits<T>::HasMultipleComponents{});
   }
 
 private:
@@ -110,7 +112,7 @@ private:
                  vtkm::VecTraitsTagSingleComponent) const
   {
     VTKM_ASSERT(componentIndex == 0);
-    using VTraits = vtkm::VecTraits<T>;
+    using VTraits = vtkm::internal::SafeVecTraits<T>;
     using TBase = typename VTraits::BaseComponentType;
     VTKM_STATIC_ASSERT(VTraits::NUM_COMPONENTS == 1);
 
@@ -133,7 +135,7 @@ private:
                  vtkm::CopyFlag allowCopy,
                  vtkm::VecTraitsTagMultipleComponents) const
   {
-    using VTraits = vtkm::VecTraits<VecType>;
+    using VTraits = vtkm::internal::SafeVecTraits<VecType>;
     using T = typename VTraits::ComponentType;
     constexpr vtkm::IdComponent N = VTraits::NUM_COMPONENTS;
 
@@ -252,10 +254,10 @@ using ArrayExtractComponentIsInefficient = typename std::is_base_of<
 /// `vtkm::cont::internal::ArrayExtractComponentImpl`.
 ///
 template <typename T, typename S>
-vtkm::cont::ArrayHandleStride<typename vtkm::VecTraits<T>::BaseComponentType> ArrayExtractComponent(
+vtkm::cont::ArrayHandleStride<typename vtkm::internal::SafeVecTraits<T>::BaseComponentType>
-  const vtkm::cont::ArrayHandle<T, S>& src,
+ArrayExtractComponent(const vtkm::cont::ArrayHandle<T, S>& src,
-  vtkm::IdComponent componentIndex,
+                      vtkm::IdComponent componentIndex,
-  vtkm::CopyFlag allowCopy = vtkm::CopyFlag::On)
+                      vtkm::CopyFlag allowCopy = vtkm::CopyFlag::On)
 {
   return internal::ArrayExtractComponentImpl<S>{}(src, componentIndex, allowCopy);
 }

vtkm/cont/ArrayHandle.h

@@ -831,9 +831,9 @@ VTKM_NEVER_EXPORT VTKM_CONT inline void printSummary_ArrayHandle_Value(
   std::ostream& out,
   vtkm::VecTraitsTagMultipleComponents)
 {
-  using Traits = vtkm::VecTraits<T>;
+  using Traits = vtkm::internal::SafeVecTraits<T>;
   using ComponentType = typename Traits::ComponentType;
-  using IsVecOfVec = typename vtkm::VecTraits<ComponentType>::HasMultipleComponents;
+  using IsVecOfVec = typename vtkm::internal::SafeVecTraits<ComponentType>::HasMultipleComponents;
   vtkm::IdComponent numComponents = Traits::GetNumberOfComponents(value);
   out << "(";
   printSummary_ArrayHandle_Value(Traits::GetComponent(value, 0), out, IsVecOfVec());
@@ -853,10 +853,10 @@ VTKM_NEVER_EXPORT VTKM_CONT inline void printSummary_ArrayHandle_Value(
 {
   out << "{";
   printSummary_ArrayHandle_Value(
-    value.first, out, typename vtkm::VecTraits<T1>::HasMultipleComponents());
+    value.first, out, typename vtkm::internal::SafeVecTraits<T1>::HasMultipleComponents());
   out << ",";
   printSummary_ArrayHandle_Value(
-    value.second, out, typename vtkm::VecTraits<T2>::HasMultipleComponents());
+    value.second, out, typename vtkm::internal::SafeVecTraits<T2>::HasMultipleComponents());
   out << "}";
 }
 
@@ -872,7 +872,7 @@ VTKM_NEVER_EXPORT VTKM_CONT inline void printSummary_ArrayHandle(
 {
   using ArrayType = vtkm::cont::ArrayHandle<T, StorageT>;
   using PortalType = typename ArrayType::ReadPortalType;
-  using IsVec = typename vtkm::VecTraits<T>::HasMultipleComponents;
+  using IsVec = typename vtkm::internal::SafeVecTraits<T>::HasMultipleComponents;
 
   vtkm::Id sz = array.GetNumberOfValues();
 

vtkm/cont/UnknownArrayHandle.h

@@ -55,12 +55,12 @@ vtkm::Id UnknownAHNumberOfValues(void* mem)
   return arrayHandle->GetNumberOfValues();
 }
 
-template <typename T, typename StaticSize = typename vtkm::VecTraits<T>::IsSizeStatic>
+template <typename T, typename StaticSize = typename vtkm::internal::SafeVecTraits<T>::IsSizeStatic>
 struct UnknownAHNumberOfComponentsImpl;
 template <typename T>
 struct UnknownAHNumberOfComponentsImpl<T, vtkm::VecTraitsTagSizeStatic>
 {
-  static constexpr vtkm::IdComponent Value = vtkm::VecTraits<T>::NUM_COMPONENTS;
+  static constexpr vtkm::IdComponent Value = vtkm::internal::SafeVecTraits<T>::NUM_COMPONENTS;
 };
 template <typename T>
 struct UnknownAHNumberOfComponentsImpl<T, vtkm::VecTraitsTagSizeVariable>
@@ -74,18 +74,25 @@ vtkm::IdComponent UnknownAHNumberOfComponents()
   return UnknownAHNumberOfComponentsImpl<T>::Value;
 }
 
-template <typename T, typename StaticSize = typename vtkm::VecTraits<T>::IsSizeStatic>
+template <typename T,
+          typename = typename vtkm::internal::SafeVecTraits<T>::IsSizeStatic,
+          typename = vtkm::HasVecTraits<T>>
 struct UnknownAHNumberOfComponentsFlatImpl;
 template <typename T>
-struct UnknownAHNumberOfComponentsFlatImpl<T, vtkm::VecTraitsTagSizeStatic>
+struct UnknownAHNumberOfComponentsFlatImpl<T, vtkm::VecTraitsTagSizeStatic, std::true_type>
 {
   static constexpr vtkm::IdComponent Value = vtkm::VecFlat<T>::NUM_COMPONENTS;
 };
 template <typename T>
-struct UnknownAHNumberOfComponentsFlatImpl<T, vtkm::VecTraitsTagSizeVariable>
+struct UnknownAHNumberOfComponentsFlatImpl<T, vtkm::VecTraitsTagSizeVariable, std::true_type>
 {
   static constexpr vtkm::IdComponent Value = 0;
 };
+template <typename T>
+struct UnknownAHNumberOfComponentsFlatImpl<T, vtkm::VecTraitsTagSizeStatic, std::false_type>
+{
+  static constexpr vtkm::IdComponent Value = 1;
+};
 
 template <typename T>
 vtkm::IdComponent UnknownAHNumberOfComponentsFlat()
@@ -306,13 +313,14 @@ std::shared_ptr<UnknownAHContainer> UnknownAHNewInstanceBasic(vtkm::VecTraitsTag
 template <typename T>
 std::shared_ptr<UnknownAHContainer> UnknownAHNewInstanceBasic()
 {
-  return UnknownAHNewInstanceBasic<T>(typename vtkm::VecTraits<T>::IsSizeStatic{});
+  return UnknownAHNewInstanceBasic<T>(typename vtkm::internal::SafeVecTraits<T>::IsSizeStatic{});
 }
 
 template <typename T>
 std::shared_ptr<UnknownAHContainer> UnknownAHNewInstanceFloatBasic(vtkm::VecTraitsTagSizeStatic)
 {
-  using FloatT = typename vtkm::VecTraits<T>::template ReplaceBaseComponentType<vtkm::FloatDefault>;
+  using FloatT = typename vtkm::internal::SafeVecTraits<T>::template ReplaceBaseComponentType<
+    vtkm::FloatDefault>;
   return UnknownAHContainer::Make(vtkm::cont::ArrayHandleBasic<FloatT>{});
 }
 template <typename T>
@@ -324,7 +332,8 @@ std::shared_ptr<UnknownAHContainer> UnknownAHNewInstanceFloatBasic(vtkm::VecTrai
 template <typename T>
 std::shared_ptr<UnknownAHContainer> UnknownAHNewInstanceFloatBasic()
 {
-  return UnknownAHNewInstanceFloatBasic<T>(typename vtkm::VecTraits<T>::IsSizeStatic{});
+  return UnknownAHNewInstanceFloatBasic<T>(
+    typename vtkm::internal::SafeVecTraits<T>::IsSizeStatic{});
 }
 
 template <typename T, typename S>
@@ -333,7 +342,7 @@ inline UnknownAHContainer::UnknownAHContainer(const vtkm::cont::ArrayHandle<T, S
   , ValueType(typeid(T))
   , StorageType(typeid(S))
   , BaseComponentType(
-      UnknownAHComponentInfo::Make<typename vtkm::VecTraits<T>::BaseComponentType>())
+      UnknownAHComponentInfo::Make<typename vtkm::internal::SafeVecTraits<T>::BaseComponentType>())
   , DeleteFunction(detail::UnknownAHDelete<T, S>)
   , NewInstance(detail::UnknownAHNewInstance<T, S>)
   , NewInstanceBasic(detail::UnknownAHNewInstanceBasic<T>)

vtkm/filter/FilterParticleAdvection.h

@@ -90,7 +90,8 @@ protected:
   */
 
   vtkm::Id NumberOfSteps;
-  vtkm::filter::particleadvection::ParticleAdvectionResultType ResultType;
+  vtkm::filter::particleadvection::ParticleAdvectionResultType ResultType =
+    vtkm::filter::particleadvection::UNKNOWN_TYPE;
   vtkm::cont::ArrayHandle<ParticleType> Seeds;
   vtkm::filter::particleadvection::IntegrationSolverType SolverType;
   vtkm::FloatDefault StepSize;

vtkm/filter/ParticleAdvection2.h

@@ -13,6 +13,8 @@
 
 #include <vtkm/Particle.h>
 #include <vtkm/filter/FilterParticleAdvection.h>
+#include <vtkm/filter/NewFilterField.h>
+#include <vtkm/filter/particleadvection/ParticleAdvectionTypes.h>
 
 namespace vtkm
 {
@@ -38,6 +40,19 @@ public:
   vtkm::cont::UnknownArrayHandle SeedArray;
 };
 
+class ParticleAdvection3 : public vtkm::filter::NewFilterField
+{
+public:
+  //  VTKM_CONT ParticleAdvection3() {}
+
+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;
+};
+
 }
 } // namespace vtkm::filter
 

vtkm/filter/ParticleAdvection2.hxx

@@ -29,6 +29,7 @@ namespace filter
 inline VTKM_CONT ParticleAdvection2::ParticleAdvection2()
   : vtkm::filter::FilterParticleAdvection<ParticleAdvection2, vtkm::Particle>()
 {
+  this->ResultType = vtkm::filter::particleadvection::PARTICLE_ADVECT_TYPE;
 }
 
 //-----------------------------------------------------------------------------
@@ -56,21 +57,17 @@ inline VTKM_CONT vtkm::cont::PartitionedDataSet ParticleAdvection2::PrepareForEx
     if (!ds.HasPointField(activeField) && !ds.HasCellField(activeField))
       throw vtkm::cont::ErrorFilterExecution("Unsupported field assocation");
 
-    dsi.push_back(DSIType(ds, blockId, activeField, this->SolverType, this->VecFieldType));
+    dsi.push_back(
+      DSIType(ds, blockId, activeField, this->SolverType, this->VecFieldType, this->ResultType));
   }
 
-  std::vector<vtkm::Particle> seeds;
+  this->SeedArray = this->Seeds;
-  seeds.push_back(this->Seeds.ReadPortal().Get(0));
+  vtkm::filter::particleadvection::PAV pav(
-  seeds.push_back(this->Seeds.ReadPortal().Get(1));
+    boundsMap,
-  seeds.push_back(this->Seeds.ReadPortal().Get(2));
+    dsi,
-  std::cout << "SEEDS= " << seeds[0].Pos << " " << seeds[1].Pos << " " << seeds[2].Pos << std::endl;
+    this->UseThreadedAlgorithm,
-
+    vtkm::filter::particleadvection::ParticleAdvectionResultType::PARTICLE_ADVECT_TYPE);
-  //this->SeedArray = this->Seeds;
+  return pav.Execute(this->NumberOfSteps, this->StepSize, this->SeedArray);
-  //vtkm::cont::UncertainArrayHandle<vtkm::List<vtkm::Particle>, vtkm::cont::StorageListBasic> arr;
-  //arr = this->Seeds;
-
-  vtkm::filter::particleadvection::PAV pav(boundsMap, dsi, this->UseThreadedAlgorithm);
-  return pav.Execute(this->NumberOfSteps, this->StepSize, this->Seeds);
 
 #if 0
   //std::vector<DSIType> ddsi;
@@ -97,6 +94,20 @@ inline VTKM_CONT vtkm::cont::PartitionedDataSet ParticleAdvection2::PrepareForEx
 #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;
+}
+
 }
 } // namespace vtkm::filter
 #endif

vtkm/filter/particleadvection/ABA.h

@@ -48,8 +48,13 @@ public:
 
   vtkm::cont::PartitionedDataSet GetOutput()
   {
+    vtkm::cont::PartitionedDataSet output;
+
+    for (const auto& b : this->Blocks)
+      output.AppendPartition(b.GetOutput<ParticleType>());
+
     std::cout << "GetOutput: " << __FILE__ << " " << __LINE__ << std::endl;
-    auto output = internal::ResultHelper<ResultType, ParticleType>::GetOutput(this->Results);
+    //auto output = internal::ResultHelper<ResultType, ParticleType>::GetOutput(this->Results);
     output.PrintSummary(std::cout);
 
     return output;
@@ -100,15 +105,17 @@ public:
 
     while (this->TotalNumTerminatedParticles < this->TotalNumParticles)
     {
-      std::cout << __FILE__ << " " << __LINE__ << std::endl;
       std::vector<ParticleType> v;
       vtkm::Id numTerm = 0, blockId = -1;
-      if (GetActiveParticles(v, blockId))
+      if (this->GetActiveParticles(v, blockId))
       {
-        const auto& block = this->GetDataSet(blockId);
+        //make this a pointer to avoid the copy?
-        ResultType<ParticleType> r;
+        auto& block = this->GetDataSet(blockId);
-        block.Advect(v, this->StepSize, this->NumberOfSteps, r);
+        //block.Advect(v, this->StepSize, this->NumberOfSteps, r);
-        numTerm = this->UpdateResult(r, blockId);
+        DSIStuff<ParticleType> stuff(this->BoundsMap, this->ParticleBlockIDsMap);
+        block.Advect(v, this->StepSize, this->NumberOfSteps, stuff);
+        numTerm = this->UpdateResult(stuff);
+        //numTerm = this->UpdateResult(r, blockId);
         std::cout << " Advect: " << v.size() << " NT= " << numTerm << std::endl;
       }
 
@@ -138,9 +145,9 @@ public:
     this->TotalNumParticles = static_cast<vtkm::Id>(total);
   }
 
-  const DSI& GetDataSet(vtkm::Id id) const
+  DSI& GetDataSet(vtkm::Id id)
   {
-    for (const auto& it : this->Blocks)
+    for (auto& it : this->Blocks)
       if (it.GetID() == id)
         return it;
     throw vtkm::cont::ErrorFilterExecution("Bad block");
@@ -231,6 +238,7 @@ public:
       this->ParticleBlockIDsMap[it.first] = it.second;
   }
 
+  /*
   void UpdateTerminated(const vtkm::cont::ArrayHandle<ParticleType>& particles,
                         const std::vector<vtkm::Id>& idxTerm)
   {
@@ -238,7 +246,9 @@ public:
     for (const auto& idx : idxTerm)
       this->ParticleBlockIDsMap.erase(portal.Get(idx).ID);
   }
+  */
 
+  /*
   void ClassifyParticles(const vtkm::cont::ArrayHandle<ParticleType>& particles,
                          std::unordered_map<vtkm::Id, std::vector<vtkm::Id>>& idsMapA,
                          std::unordered_map<vtkm::Id, std::vector<vtkm::Id>>& idsMapI,
@@ -319,10 +329,28 @@ public:
     //Make sure we didn't miss anything. Every particle goes into a single bucket.
     VTKM_ASSERT(static_cast<std::size_t>(n) == (A.size() + I.size() + termIdx.size()));
   }
+  */
 
+  vtkm::Id UpdateResult(const DSIStuff<ParticleType>& stuff)
+  {
+    this->UpdateActive(stuff.A, stuff.IdMapA);
+    this->UpdateInactive(stuff.I, stuff.IdMapI);
+
+    vtkm::Id numTerm = static_cast<vtkm::Id>(stuff.TermID.size());
+    //Update terminated particles.
+    if (numTerm > 0)
+    {
+      for (const auto& id : stuff.TermID)
+        this->ParticleBlockIDsMap.erase(id);
+    }
+
+    return numTerm;
+  }
+
+  /*
   vtkm::Id UpdateResult(ResultType<ParticleType>& res, vtkm::Id blockId)
   {
-    std::cout << "UpdateResult: blockId= " << blockId << std::endl;
+    std::cout<<"UpdateResult: blockId= "<<blockId<<std::endl;
     std::unordered_map<vtkm::Id, std::vector<vtkm::Id>> idsMapI, idsMapA;
     std::vector<ParticleType> A, I;
     std::vector<vtkm::Id> termIdx;
@@ -340,6 +368,7 @@ public:
 
     return numTerm;
   }
+  */
 
   virtual bool GetBlockAndWait(const vtkm::Id& numLocalTerm)
   {

vtkm/filter/particleadvection/DSI.h

@@ -19,6 +19,8 @@
 #include <vtkm/worklet/particleadvection/RK4Integrator.h>
 #include <vtkm/worklet/particleadvection/Stepper.h>
 
+#include <vtkm/exec/internal/Variant.h>
+
 namespace vtkm
 {
 namespace filter
@@ -26,6 +28,23 @@ namespace filter
 namespace particleadvection
 {
 
+template <typename ParticleType>
+struct DSIStuff
+{
+  DSIStuff(const vtkm::filter::particleadvection::BoundsMap& boundsMap,
+           const std::unordered_map<vtkm::Id, std::vector<vtkm::Id>>& particleBlockIDsMap)
+    : BoundsMap(boundsMap)
+    , ParticleBlockIDsMap(particleBlockIDsMap)
+  {
+  }
+
+  std::vector<ParticleType> A, I;
+  std::unordered_map<vtkm::Id, std::vector<vtkm::Id>> IdMapA, IdMapI;
+  std::vector<vtkm::Id> TermID;
+  const vtkm::filter::particleadvection::BoundsMap BoundsMap;
+  const std::unordered_map<vtkm::Id, std::vector<vtkm::Id>> ParticleBlockIDsMap;
+};
+
 class DSI
 {
 public:
@@ -33,10 +52,13 @@ public:
       vtkm::Id id,
       const std::string& fieldNm,
       vtkm::filter::particleadvection::IntegrationSolverType solverType,
-      vtkm::filter::particleadvection::VectorFieldType vecFieldType)
+      vtkm::filter::particleadvection::VectorFieldType vecFieldType,
+      vtkm::filter::particleadvection::ParticleAdvectionResultType resultType)
+    : Rank(this->Comm.rank())
   {
     this->SolverType = solverType;
     this->VecFieldType = vecFieldType;
+    this->ResType = resultType;
     this->DataSet = ds;
     this->Id = id;
     this->FieldName = fieldNm;
@@ -47,12 +69,19 @@ public:
   VTKM_CONT vtkm::Id GetID() const { return this->Id; }
   VTKM_CONT void SetCopySeedFlag(bool val) { this->CopySeedArray = val; }
 
-  //template <typename ParticleType, typename ResultType>
+
-  template <typename ParticleType, template <typename> class ResultType>
+  template <typename ParticleType>
   VTKM_CONT void Advect(std::vector<ParticleType>& v,
                         vtkm::FloatDefault stepSize, //move these to member data(?)
                         vtkm::Id maxSteps,
-                        ResultType<ParticleType>& result) const;
+                        DSIStuff<ParticleType>& stuff);
+
+  template <typename ParticleType>
+  VTKM_CONT vtkm::cont::DataSet GetOutput() const;
+
+protected:
+  template <typename ParticleType, template <typename> class ResultType>
+  VTKM_CONT void UpdateResult(ResultType<ParticleType>& result, DSIStuff<ParticleType>& stuff);
 
   template <typename ArrayType>
   VTKM_CONT void GetVelocityField(
@@ -69,6 +98,7 @@ public:
   VTKM_CONT void GetElectroMagneticField(
     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(),
@@ -79,14 +109,31 @@ public:
       vtkm::worklet::particleadvection::ElectroMagneticField<ArrayType>(arr1, arr2, assoc);
   }
 
-  //private:
+  template <typename ParticleType>
+  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;
   vtkm::filter::particleadvection::VectorFieldType VecFieldType;
+  vtkm::filter::particleadvection::ParticleAdvectionResultType ResType =
+    vtkm::filter::particleadvection::UNKNOWN_TYPE;
   bool CopySeedArray = false;
-  //Add results array?
+
+  using ResultVariantType =
+    vtkm::exec::internal::Variant<vtkm::worklet::ParticleAdvectionResult<vtkm::Particle>,
+                                  vtkm::worklet::StreamlineResult<vtkm::Particle>,
+                                  vtkm::worklet::ParticleAdvectionResult<vtkm::ChargedParticle>,
+                                  vtkm::worklet::StreamlineResult<vtkm::ChargedParticle>>;
+
+  std::vector<ResultVariantType> Results;
 };
 
 }

vtkm/filter/particleadvection/DSI.hxx

@@ -17,12 +17,264 @@ namespace filter
 {
 namespace particleadvection
 {
+#if 0
+namespace internal
+{
+//Helper class to store the different result types.
+//template <typename ResultType>
+template <template <typename> class ResultType, typename ParticleType>
+class ResultHelper2;
+
+template <typename ParticleType>
+class ResultHelper2<vtkm::worklet::ParticleAdvectionResult, ParticleType>
+{
+  using PAType = vtkm::worklet::ParticleAdvectionResult<ParticleType>;
+
+public:
+  static void Store(const PAType& result)
+  {
+  }
+
+  static vtkm::cont::DataSet
+  GetOutput()
+  {
+    vtkm::cont::DataSet ds;
+
+    return ds;
+  }
+};
+
+template <typename ParticleType>
+class ResultHelper2<vtkm::worklet::StreamlineResult, ParticleType>
+{
+  using PAType = vtkm::worklet::StreamlineResult<ParticleType>;
+
+public:
+  static void Store(const PAType& result)
+  {
+  }
+
+  static vtkm::cont::DataSet
+  GetOutput()
+  {
+    vtkm::cont::DataSet ds;
+
+    return ds;
+  }
+};
+
+}
+#endif
+
+template <typename ParticleType>
+VTKM_CONT vtkm::cont::DataSet DSI::GetOutput() const
+{
+  vtkm::cont::DataSet ds;
+
+  std::size_t nResults = this->Results.size();
+  if (nResults == 0)
+    return ds;
+
+  std::cout << "Check the variant type!" << std::endl;
+
+  if (this->ResType == PARTICLE_ADVECT_TYPE)
+  {
+    using RType = vtkm::worklet::ParticleAdvectionResult<ParticleType>;
+
+    std::vector<vtkm::cont::ArrayHandle<ParticleType>> allParticles;
+    allParticles.reserve(nResults);
+
+    //Get all the particles and put them into a single ArrayHandle: pts
+    for (const auto& vres : this->Results)
+      allParticles.push_back(vres.Get<RType>().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->ResType == STREAMLINE_TYPE)
+  {
+    using ResultType = vtkm::worklet::StreamlineResult<ParticleType>;
+
+    //Easy case with one result.
+    if (nResults == 1)
+    {
+      const auto& res = this->Results[0].Get<ResultType>();
+      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].Get<ResultType>();
+        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].Get<ResultType>();
+        // 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].Get<ResultType>();
+        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 ds;
+}
 
 template <typename ParticleType, template <typename> class ResultType>
+VTKM_CONT void DSI::UpdateResult(ResultType<ParticleType>& result, DSIStuff<ParticleType>& stuff)
+{
+  this->ClassifyParticles(result.Particles, stuff);
+  this->Results.push_back(result);
+}
+
+template <typename ParticleType>
+VTKM_CONT void DSI::ClassifyParticles(const vtkm::cont::ArrayHandle<ParticleType>& particles,
+                                      DSIStuff<ParticleType>& stuff) const
+{
+  stuff.A.clear();
+  stuff.I.clear();
+  stuff.TermID.clear();
+  stuff.IdMapI.clear();
+  stuff.IdMapA.clear();
+
+  auto portal = particles.WritePortal();
+  vtkm::Id n = portal.GetNumberOfValues();
+
+  for (vtkm::Id i = 0; i < n; i++)
+  {
+    auto p = portal.Get(i);
+
+    if (p.Status.CheckTerminate())
+      stuff.TermID.push_back(p.ID);
+    else
+    {
+      const auto& it = stuff.ParticleBlockIDsMap.find(p.ID);
+      VTKM_ASSERT(it != stuff.ParticleBlockIDsMap.end());
+      auto currBIDs = it->second;
+      VTKM_ASSERT(!currBIDs.empty());
+
+      std::vector<vtkm::Id> newIDs;
+      if (p.Status.CheckSpatialBounds() && !p.Status.CheckTookAnySteps())
+        newIDs.assign(std::next(currBIDs.begin(), 1), currBIDs.end());
+      else
+        newIDs = stuff.BoundsMap.FindBlocks(p.Pos, currBIDs);
+
+      //reset the particle status.
+      p.Status = vtkm::ParticleStatus();
+
+      if (newIDs.empty()) //No blocks, we're done.
+      {
+        p.Status.SetTerminate();
+        stuff.TermID.push_back(p.ID);
+      }
+      else
+      {
+        //If we have more than blockId, we want to minimize communication
+        //and put any blocks owned by this rank first.
+        if (newIDs.size() > 1)
+        {
+          for (auto idit = newIDs.begin(); idit != newIDs.end(); idit++)
+          {
+            vtkm::Id bid = *idit;
+            if (stuff.BoundsMap.FindRank(bid) == this->Rank)
+            {
+              newIDs.erase(idit);
+              newIDs.insert(newIDs.begin(), bid);
+              break;
+            }
+          }
+        }
+
+        int dstRank = stuff.BoundsMap.FindRank(newIDs[0]);
+        if (dstRank == this->Rank)
+        {
+          stuff.A.push_back(p);
+          stuff.IdMapA[p.ID] = newIDs;
+        }
+        else
+        {
+          stuff.I.push_back(p);
+          stuff.IdMapI[p.ID] = newIDs;
+        }
+      }
+      portal.Set(i, p);
+    }
+  }
+
+  //Make sure we didn't miss anything. Every particle goes into a single bucket.
+  VTKM_ASSERT(static_cast<std::size_t>(n) ==
+              (stuff.A.size() + stuff.I.size() + stuff.TermID.size()));
+}
+
+template <typename ParticleType>
 VTKM_CONT void DSI::Advect(std::vector<ParticleType>& v,
                            vtkm::FloatDefault stepSize,
                            vtkm::Id maxSteps,
-                           ResultType<ParticleType>& result) const
+                           DSIStuff<ParticleType>& stuff)
 {
   //using Association = vtkm::cont::Field::Association;
   using ArrayType = vtkm::cont::ArrayHandle<vtkm::Vec3f>;
@@ -46,9 +298,22 @@ VTKM_CONT void DSI::Advect(std::vector<ParticleType>& v,
 
       GridEvType eval(this->DataSet, velField);
       StepperType stepper(eval, stepSize);
-      vtkm::worklet::ParticleAdvection Worklet;
+
-      result = Worklet.Run(stepper, seedArray, maxSteps);
+      //make this a template
-      //Put results in unknown array??
+      if (this->ResType == PARTICLE_ADVECT_TYPE)
+      {
+        vtkm::worklet::ParticleAdvection Worklet;
+        auto r = Worklet.Run(stepper, seedArray, maxSteps);
+        this->UpdateResult(r, stuff);
+        // /*result =*/ Worklet.Run(stepper, seedArray, maxSteps);
+      }
+      else
+      {
+        vtkm::worklet::Streamline Worklet;
+        auto r = Worklet.Run(stepper, seedArray, maxSteps);
+        this->UpdateResult(r, stuff);
+        //Put results in unknown array??
+      }
     }
     else if (this->VecFieldType == ELECTRO_MAGNETIC_FIELD_TYPE) //vtkm::ChargedParticle
     {
@@ -62,8 +327,21 @@ VTKM_CONT void DSI::Advect(std::vector<ParticleType>& v,
 
       GridEvType eval(this->DataSet, velField);
       StepperType stepper(eval, stepSize);
-      vtkm::worklet::ParticleAdvection Worklet;
+
-      result = Worklet.Run(stepper, seedArray, maxSteps);
+      if (this->ResType == PARTICLE_ADVECT_TYPE)
+      {
+        vtkm::worklet::ParticleAdvection Worklet;
+        auto r = Worklet.Run(stepper, seedArray, maxSteps);
+        this->UpdateResult(r, stuff);
+        ///*result =*/ auto r = Worklet.Run(stepper, seedArray, maxSteps);
+      }
+      else
+      {
+        vtkm::worklet::Streamline Worklet;
+        auto r = Worklet.Run(stepper, seedArray, maxSteps);
+        this->UpdateResult(r, stuff);
+        ///*result =*/ Worklet.Run(stepper, seedArray, maxSteps);
+      }
     }
   }
 
@@ -81,8 +359,8 @@ VTKM_CONT void DSI::Advect(std::vector<ParticleType>& v,
 
       GridEvType eval(this->DataSet, velField);
       StepperType stepper(eval, stepSize);
-      vtkm::worklet::ParticleAdvection Worklet;
+      //      vtkm::worklet::ParticleAdvection Worklet;
-      result = Worklet.Run(stepper, seedArray, maxSteps);
+      //      result = Worklet.Run(stepper, seedArray, maxSteps);
       //Put results in unknown array??
     }
   }

vtkm/filter/particleadvection/PAV.h

@@ -27,38 +27,78 @@ class PAV
 public:
   PAV(const vtkm::filter::particleadvection::BoundsMap& bm,
       const std::vector<DSI>& blocks,
-      const bool& useThreaded)
+      const bool& useThreaded,
+      const vtkm::filter::particleadvection::ParticleAdvectionResultType& parType)
     : Blocks(blocks)
     , BoundsMap(bm)
+    , ResultType(parType)
     , UseThreadedAlgorithm(useThreaded)
   {
   }
 
   vtkm::cont::PartitionedDataSet Execute(vtkm::Id numSteps,
                                          vtkm::FloatDefault stepSize,
-                                         const vtkm::cont::ArrayHandle<vtkm::Particle>& seeds)
+                                         const vtkm::cont::UnknownArrayHandle& seeds)
+  {
+    using ParticleArray = vtkm::cont::ArrayHandle<vtkm::Particle>;
+    using ChargedParticleArray = vtkm::cont::ArrayHandle<vtkm::ChargedParticle>;
+
+    if (seeds.IsBaseComponentType<vtkm::Particle>())
+      return this->Execute(numSteps, stepSize, seeds.AsArrayHandle<ParticleArray>());
+    else if (seeds.IsBaseComponentType<vtkm::ChargedParticle>())
+      return this->Execute(numSteps, stepSize, seeds.AsArrayHandle<ChargedParticleArray>());
+
+    throw vtkm::cont::ErrorFilterExecution("Unsupported options in ABA");
+  }
+
+private:
+  template <typename ParticleType>
+  vtkm::cont::PartitionedDataSet Execute(vtkm::Id numSteps,
+                                         vtkm::FloatDefault stepSize,
+                                         const vtkm::cont::ArrayHandle<ParticleType>& seeds)
   {
     if (!this->UseThreadedAlgorithm)
     {
+      //make a templated algorithm execution()
+      if (this->ResultType ==
+          vtkm::filter::particleadvection::ParticleAdvectionResultType::PARTICLE_ADVECT_TYPE)
+      {
+        using AlgorithmType =
+          vtkm::filter::particleadvection::ABA<vtkm::worklet::ParticleAdvectionResult,
+                                               ParticleType>;
+
+        AlgorithmType algo(this->BoundsMap, this->Blocks);
+        algo.Execute(numSteps, stepSize, seeds);
+        return algo.GetOutput();
+      }
+      else
+      {
+        using AlgorithmType =
+          vtkm::filter::particleadvection::ABA<vtkm::worklet::StreamlineResult, ParticleType>;
+
+        AlgorithmType algo(this->BoundsMap, this->Blocks);
+        algo.Execute(numSteps, stepSize, seeds);
+        return algo.GetOutput();
+      }
+    }
+    else
+    {
+      std::cout << "Change me to threaded ABA" << std::endl;
       using AlgorithmType =
-        vtkm::filter::particleadvection::ABA<vtkm::worklet::ParticleAdvectionResult,
+        vtkm::filter::particleadvection::ABA<vtkm::worklet::ParticleAdvectionResult, ParticleType>;
-                                             vtkm::Particle>;
 
       AlgorithmType algo(this->BoundsMap, this->Blocks);
       algo.Execute(numSteps, stepSize, seeds);
       return algo.GetOutput();
     }
-    else
-    {
-      //using ThreadedAlgorithmType = vtkm::filter::particleadvection::ParticleAdvectionThreadedAlgorithm;
-      //AlgorithmType algo;
-    }
-
-    throw vtkm::cont::ErrorFilterExecution("Unsupported options in ABA");
   }
 
+
+
+
   std::vector<DSI> Blocks;
   vtkm::filter::particleadvection::BoundsMap BoundsMap;
+  ParticleAdvectionResultType ResultType;
   bool UseThreadedAlgorithm;
 };
 

vtkm/filter/particleadvection/ParticleAdvectionTypes.h

@@ -30,7 +30,8 @@ enum VectorFieldType
 
 enum ParticleAdvectionResultType
 {
-  PARTICLE_ADVECT_TYPE = 0,
+  UNKNOWN_TYPE = -1,
+  PARTICLE_ADVECT_TYPE,
   STREAMLINE_TYPE,
 };
 

vtkm/filter/testing/UnitTestStreamlineFilter.cxx

@@ -65,7 +65,6 @@ void TestStreamline()
                                      vtkm::Particle(vtkm::Vec3f(.2f, 3.0f, .2f), 2) });
 
     vtkm::filter::ParticleAdvection2 streamline;
-
     streamline.SetStepSize(0.1f);
     streamline.SetNumberOfSteps(20);
     streamline.SetSeeds(seedArray);
@@ -584,7 +583,7 @@ void TestStreamlineFilters()
   vtkm::filter::ParticleAdvection2 filter;
   TestStreamline();
 
-  /*
+
   std::vector<bool> flags = { true, false };
   std::vector<FilterType> fTypes = { FilterType::PARTICLE_ADVECTION,
                                      FilterType::STREAMLINE,
@@ -601,7 +600,7 @@ void TestStreamlineFilters()
 
   TestStreamline();
   TestPathline();
-
+  /*
   for (auto useSL : flags)
     TestAMRStreamline(useSL);