consolidate two particle density filters

vtkm/filter/ParticleDensityBase.h

@@ -0,0 +1,113 @@
+//============================================================================
+//  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_particle_density_base_h
+#define vtk_m_filter_particle_density_base_h
+
+#include <vtkm/filter/FilterField.h>
+#include <vtkm/worklet/WorkletMapField.h>
+
+namespace vtkm
+{
+namespace filter
+{
+// We only need the CoordinateSystem and scalar fields of the input dataset thus a FilterField
+template <typename Derived>
+class ParticleDensityBase : public vtkm::filter::FilterField<Derived>
+{
+public:
+  // deposit scalar field associated with particles, e.g. mass/charge to mesh cells
+  using SupportedTypes = vtkm::TypeListFieldScalar;
+
+protected:
+  ParticleDensityBase(const vtkm::Id3& dimension,
+                      const vtkm::Vec3f& origin,
+                      const vtkm::Vec3f& spacing)
+    : Dimension(dimension)
+    , Origin(origin)
+    , Spacing(spacing)
+    , ComputeNumberDensity(false)
+    , DivideByVolume(true)
+  {
+  }
+
+  ParticleDensityBase(const vtkm::Id3& dimension, const vtkm::Bounds& bounds)
+    : Dimension(dimension)
+    , Origin({ static_cast<vtkm::FloatDefault>(bounds.X.Min),
+               static_cast<vtkm::FloatDefault>(bounds.Y.Min),
+               static_cast<vtkm::FloatDefault>(bounds.Z.Min) })
+    , Spacing(vtkm::Vec3f{ static_cast<vtkm::FloatDefault>(bounds.X.Length()),
+                           static_cast<vtkm::FloatDefault>(bounds.Y.Length()),
+                           static_cast<vtkm::FloatDefault>(bounds.Z.Length()) } /
+              Dimension)
+    , ComputeNumberDensity(false)
+    , DivideByVolume(true)
+  {
+  }
+
+public:
+  template <typename DerivedPolicy>
+  VTKM_CONT vtkm::cont::DataSet PrepareForExecution(const vtkm::cont::DataSet& input,
+                                                    vtkm::filter::PolicyBase<DerivedPolicy> policy)
+  {
+    if (this->ComputeNumberDensity)
+    {
+      return static_cast<Derived*>(this)->DoExecute(
+        input,
+        vtkm::cont::make_ArrayHandleConstant(vtkm::FloatDefault{ 1 }, input.GetNumberOfPoints()),
+        vtkm::filter::FieldMetadata{}, // Ignored
+        policy);
+    }
+    else
+    {
+      return this->FilterField<Derived>::PrepareForExecution(input, policy);
+    }
+  }
+
+  VTKM_CONT void SetComputeNumberDensity(bool yes) { this->ComputeNumberDensity = yes; }
+
+  VTKM_CONT bool GetComputeNumberDensity() const { return this->ComputeNumberDensity; }
+
+  VTKM_CONT void SetDivideByVolume(bool yes) { this->DivideByVolume = yes; }
+
+  VTKM_CONT bool GetDivideByVolume() const { return this->DivideByVolume; }
+
+protected:
+  vtkm::Id3 Dimension; // Cell dimension
+  vtkm::Vec3f Origin;
+  vtkm::Vec3f Spacing;
+  bool ComputeNumberDensity;
+  bool DivideByVolume;
+
+  class DivideByVolumeWorklet : public vtkm::worklet::WorkletMapField
+  {
+  public:
+    using ControlSignature = void(FieldInOut field);
+    using ExecutionSignature = void(_1);
+
+    VTKM_EXEC_CONT
+    explicit DivideByVolumeWorklet(vtkm::Float64 volume)
+      : Volume(volume)
+    {
+    }
+
+    template <typename T>
+    VTKM_EXEC void operator()(T& value) const
+    {
+      value = static_cast<T>(value / Volume);
+    }
+
+  private:
+    vtkm::Float64 Volume;
+  }; // class DivideByVolumeWorklet
+};
+}
+}
+#endif //vtk_m_filter_particle_density_base_h

vtkm/filter/ParticleDensityCloudInCell.h

@@ -11,37 +11,44 @@
 #ifndef vtk_m_filter_particle_density_cic_h
 #define vtk_m_filter_particle_density_cic_h
 
-#include <vtkm/filter/FilterField.h>
+#include <vtkm/filter/ParticleDensityBase.h>
 
 namespace vtkm
 {
 namespace filter
 {
 /// \brief Estimate the density of particles using the Cloud-in-Cell method
-
-// We only need the CoordinateSystem of the input dataset thus a FilterField
-class ParticleDensityCloudInCell : public vtkm::filter::FilterField<ParticleDensityCloudInCell>
+/// This filter treats the CoordinateSystem of a DataSet as positions of particles.
+/// The particles are infinitesimal in size with finite mass (or other scalar attributes
+/// such as charge). The filter estimates density by imposing a regular grid as
+/// specified in the constructor. It spread the mass of each particle to its 8 nearest
+/// neighboring grid points and summing the contribution of particles for each point
+/// in the grid.
+/// The mass of particles is established by setting the active field (using SetActiveField).
+/// Note that the "mass" can actually be another quantity. For example, you could use
+/// electrical charge in place of mass to compute the charge density.
+/// Once the sum of the mass is computed for each grid point, the mass is divided by the
+/// volume of the cell. Thus, the density will be computed as the units of the mass field
+/// per the cubic units of the coordinate system. If you just want a sum of the mass in each
+/// cell, turn off the DivideByVolumeWorklet feature of this filter.
+/// In addition, you can also simply count the number of particles in each cell by calling
+/// SetComputeNumberDensity(true).
+class ParticleDensityCloudInCell : public ParticleDensityBase<ParticleDensityCloudInCell>
 {
 public:
-  // ParticleDensity only support turning 2D/3D particle positions into density
-  // FIXME: 2D?
-  //using SupportedTypes = vtkm::TypeListFieldVec3;
-  using SupportedTypes = vtkm::TypeListFieldScalar;
+  using Superclass = ParticleDensityBase<ParticleDensityCloudInCell>;
 
   ParticleDensityCloudInCell(const vtkm::Id3& dimension,
                              const vtkm::Vec3f& origin,
                              const vtkm::Vec3f& spacing);
 
+  ParticleDensityCloudInCell(const Id3& dimension, const vtkm::Bounds& bounds);
+
   template <typename T, typename StorageType, typename Policy>
   VTKM_CONT vtkm::cont::DataSet DoExecute(const vtkm::cont::DataSet& input,
                                           const vtkm::cont::ArrayHandle<T, StorageType>& field,
                                           const vtkm::filter::FieldMetadata& fieldMeta,
                                           vtkm::filter::PolicyBase<Policy> policy);
-
-private:
-  vtkm::Id3 Dimension; // Point dimension
-  vtkm::Vec3f Origin;
-  vtkm::Vec3f Spacing;
 };
 } // filter
 } // vtkm

vtkm/filter/ParticleDensityCloudInCell.hxx

@@ -80,12 +80,15 @@ namespace filter
 inline VTKM_CONT ParticleDensityCloudInCell::ParticleDensityCloudInCell(const vtkm::Id3& dimension,
                                                                         const vtkm::Vec3f& origin,
                                                                         const vtkm::Vec3f& spacing)
-  : Dimension(dimension)
-  , Origin(origin)
-  , Spacing(spacing)
+  : Superclass(dimension, origin, spacing)
 {
 }
 
+inline VTKM_CONT ParticleDensityCloudInCell::ParticleDensityCloudInCell(const Id3& dimension,
+                                                                        const vtkm::Bounds& bounds)
+  : Superclass(dimension, bounds)
+{
+}
 
 template <typename T, typename StorageType, typename Policy>
 inline VTKM_CONT vtkm::cont::DataSet ParticleDensityCloudInCell::DoExecute(
@@ -113,6 +116,12 @@ inline VTKM_CONT vtkm::cont::DataSet ParticleDensityCloudInCell::DoExecute(
 
   this->Invoke(vtkm::worklet::CICWorklet{}, coords, field, locator, uniform.GetCellSet(), density);
 
+  if (DivideByVolume)
+  {
+    auto volume = this->Spacing[0] * this->Spacing[1] * this->Spacing[2];
+    this->Invoke(DivideByVolumeWorklet{ volume }, density);
+  }
+
   uniform.AddField(vtkm::cont::make_FieldPoint("density", density));
 
   return uniform;

vtkm/filter/ParticleDensityNearestGridPoint.h

@@ -11,7 +11,7 @@
 #ifndef vtk_m_filter_particle_density_ngp_h
 #define vtk_m_filter_particle_density_ngp_h
 
-#include <vtkm/filter/FilterField.h>
+#include <vtkm/filter/ParticleDensityBase.h>
 
 namespace vtkm
 {
@@ -29,17 +29,13 @@ namespace filter
 /// Once the sum of the mass is computed for each grid cell, the mass is divided by the
 /// volume of the cell. Thus, the density will be computed as the units of the mass field
 /// per the cubic units of the coordinate system. If you just want a sum of the mass in each
-/// cell, turn off the DivideByVolume feature of this filter.
+/// cell, turn off the DivideByVolumeWorklet feature of this filter.
 /// In addition, you can also simply count the number of particles in each cell by calling
 /// SetComputeNumberDensity(true).
-
-// We only need the CoordinateSystem and scalar fields of the input dataset thus a FilterField
-class ParticleDensityNearestGridPoint
-  : public vtkm::filter::FilterField<ParticleDensityNearestGridPoint>
+class ParticleDensityNearestGridPoint : public ParticleDensityBase<ParticleDensityNearestGridPoint>
 {
 public:
-  // deposit scalar field associated with particles, e.g. mass/charge to mesh cells
-  using SupportedTypes = vtkm::TypeListFieldScalar;
+  using Superclass = ParticleDensityBase<ParticleDensityNearestGridPoint>;
 
   ParticleDensityNearestGridPoint(const vtkm::Id3& dimension,
                                   const vtkm::Vec3f& origin,
@@ -47,30 +43,11 @@ public:
 
   ParticleDensityNearestGridPoint(const vtkm::Id3& dimension, const vtkm::Bounds& bounds);
 
-  template <typename DerivedPolicy>
-  VTKM_CONT vtkm::cont::DataSet PrepareForExecution(const vtkm::cont::DataSet& input,
-                                                    vtkm::filter::PolicyBase<DerivedPolicy> policy);
-
   template <typename T, typename StorageType, typename Policy>
   VTKM_CONT vtkm::cont::DataSet DoExecute(const vtkm::cont::DataSet& input,
                                           const vtkm::cont::ArrayHandle<T, StorageType>& field,
                                           const vtkm::filter::FieldMetadata& fieldMeta,
                                           vtkm::filter::PolicyBase<Policy> policy);
-
-  VTKM_CONT void SetComputeNumberDensity(bool yes) { this->ComputeNumberDensity = yes; }
-
-  VTKM_CONT bool GetComputeNumberDensity() const { return this->ComputeNumberDensity; }
-
-  VTKM_CONT void SetDivideByVolume(bool yes) { this->DivideByVolume = yes; }
-
-  VTKM_CONT bool GetDivideByVolume() const { return this->DivideByVolume; }
-
-private:
-  vtkm::Id3 Dimension; // Cell dimension
-  vtkm::Vec3f Origin;
-  vtkm::Vec3f Spacing;
-  bool ComputeNumberDensity;
-  bool DivideByVolume;
 };
 }
 }

vtkm/filter/ParticleDensityNearestGridPoint.hxx

@@ -50,34 +50,9 @@ public:
     // We simply ignore that particular particle when it is not in the mesh.
   }
 }; //NGPWorklet
-namespace detail
-{
-class DividByVolume : public vtkm::worklet::WorkletMapField
-{
-public:
-  using ControlSignature = void(FieldInOut field);
-  using ExecutionSignature = void(_1);
-
-  VTKM_EXEC_CONT
-  explicit DividByVolume(vtkm::Float64 volume)
-    : Volume(volume)
-  {
-  }
-
-  template <typename T>
-  VTKM_EXEC void operator()(T& value) const
-  {
-    value = static_cast<T>(value / Volume);
-  }
-
-private:
-  vtkm::Float64 Volume;
-};
-} //detail
 } //worklet
 } //vtkm
 
-
 namespace vtkm
 {
 namespace filter
@@ -86,49 +61,17 @@ inline VTKM_CONT ParticleDensityNearestGridPoint::ParticleDensityNearestGridPoin
   const vtkm::Id3& dimension,
   const vtkm::Vec3f& origin,
   const vtkm::Vec3f& spacing)
-  : Dimension(dimension)
-  , Origin(origin)
-  , Spacing(spacing)
-  , ComputeNumberDensity(false)
-  , DivideByVolume(true)
+  : Superclass(dimension, origin, spacing)
 {
 }
 
 inline VTKM_CONT ParticleDensityNearestGridPoint::ParticleDensityNearestGridPoint(
   const Id3& dimension,
   const vtkm::Bounds& bounds)
-  : Dimension(dimension)
-  , Origin({ static_cast<vtkm::FloatDefault>(bounds.X.Min),
-             static_cast<vtkm::FloatDefault>(bounds.Y.Min),
-             static_cast<vtkm::FloatDefault>(bounds.Z.Min) })
-  , Spacing(vtkm::Vec3f{ static_cast<vtkm::FloatDefault>(bounds.X.Length()),
-                         static_cast<vtkm::FloatDefault>(bounds.Y.Length()),
-                         static_cast<vtkm::FloatDefault>(bounds.Z.Length()) } /
-            Dimension)
-  , ComputeNumberDensity(false)
-  , DivideByVolume(true)
+  : Superclass(dimension, bounds)
 {
 }
 
-template <typename DerivedPolicy>
-VTKM_CONT vtkm::cont::DataSet ParticleDensityNearestGridPoint::PrepareForExecution(
-  const vtkm::cont::DataSet& input,
-  vtkm::filter::PolicyBase<DerivedPolicy> policy)
-{
-  if (this->ComputeNumberDensity)
-  {
-    return this->DoExecute(
-      input,
-      vtkm::cont::make_ArrayHandleConstant(vtkm::FloatDefault{ 1 }, input.GetNumberOfPoints()),
-      vtkm::filter::FieldMetadata{}, // Ignored
-      policy);
-  }
-  else
-  {
-    return this->FilterField::PrepareForExecution(input, policy);
-  }
-}
-
 template <typename T, typename StorageType, typename Policy>
 inline VTKM_CONT vtkm::cont::DataSet ParticleDensityNearestGridPoint::DoExecute(
   const vtkm::cont::DataSet& dataSet,
@@ -166,7 +109,7 @@ inline VTKM_CONT vtkm::cont::DataSet ParticleDensityNearestGridPoint::DoExecute(
   if (DivideByVolume)
   {
     auto volume = this->Spacing[0] * this->Spacing[1] * this->Spacing[2];
-    this->Invoke(vtkm::worklet::detail::DividByVolume{ volume }, density);
+    this->Invoke(DivideByVolumeWorklet{ volume }, density);
   }
 
   uniform.AddField(vtkm::cont::make_FieldCell("density", density));

vtkm/filter/testing/UnitTestParticleDensity.cxx

@@ -103,7 +103,17 @@ void TestCIC()
   auto mass_result = vtkm::worklet::DescriptiveStatistics::Run(mass);
   auto density_result = vtkm::worklet::DescriptiveStatistics::Run(field);
   // Unfortunately, floating point atomics suffer from precision error more than everything else.
-  VTKM_TEST_ASSERT(test_equal(density_result.Sum(), mass_result.Sum(), 10));
+  VTKM_TEST_ASSERT(test_equal(density_result.Sum(), mass_result.Sum() * 27.0, 10));
+
+  filter.SetComputeNumberDensity(true);
+  filter.SetDivideByVolume(false);
+  auto counts = filter.Execute(dataSet);
+
+  vtkm::cont::ArrayHandle<vtkm::FloatDefault> field1;
+  counts.GetPointField("density").GetData().AsArrayHandle<vtkm::FloatDefault>(field1);
+
+  auto counts_result = vtkm::worklet::DescriptiveStatistics::Run(field1);
+  VTKM_TEST_ASSERT(test_equal(counts_result.Sum(), mass_result.N(), 0.1));
 }
 
 void TestParticleDensity()