migrate particle density filters

vtkm/filter/CMakeLists.txt

@@ -26,6 +26,8 @@ set(deprecated_headers
   Histogram.h
   NDEntropy.h
   NDHistogram.h
+  ParticleDensityCloudInCell.h
+  ParticleDensityNearestGridPoint.h
   )
 
 vtkm_declare_headers(${deprecated_headers})
@@ -102,9 +104,6 @@ set(extra_headers
   MaskPoints.h
   MeshQuality.h
   MIRFilter.h
-  ParticleDensityBase.h
-  ParticleDensityCloudInCell.h
-  ParticleDensityNearestGridPoint.h
   ParticleAdvection.h
   Pathline.h
   PathParticle.h
@@ -152,8 +151,6 @@ set(extra_header_template_sources
   Mask.hxx
   MeshQuality.hxx
   MIRFilter.hxx
-  ParticleDensityCloudInCell.hxx
-  ParticleDensityNearestGridPoint.hxx
   ParticleAdvection.hxx
   Pathline.hxx
   PathParticle.hxx

vtkm/filter/ParticleDensityBase.h

@@ -1,124 +0,0 @@
-//============================================================================
-//  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/FilterDataSetWithField.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::FilterDataSetWithField<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->FilterDataSetWithField<Derived>::PrepareForExecution(input, policy);
-    }
-  }
-
-  template <typename T, typename StorageType, typename Policy>
-  VTKM_CONT bool DoMapField(vtkm::cont::DataSet&,
-                            const vtkm::cont::ArrayHandle<T, StorageType>&,
-                            const vtkm::filter::FieldMetadata&,
-                            vtkm::filter::PolicyBase<Policy>)
-  {
-    return false;
-  }
-
-  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;
-
-public:
-  // conceptually protected but CUDA needs this to be public
-  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

@@ -7,51 +7,34 @@
 //  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
 //  PURPOSE.  See the above copyright notice for more information.
 //============================================================================
+#ifndef vtk_m_filter_ParticleDensityCloudInCell_h
+#define vtk_m_filter_ParticleDensityCloudInCell_h
 
-#ifndef vtk_m_filter_particle_density_cic_h
+#include <vtkm/Deprecated.h>
-#define vtk_m_filter_particle_density_cic_h
+#include <vtkm/filter/density_estimate/ParticleDensityCloudInCell.h>
-
-#include <vtkm/filter/ParticleDensityBase.h>
 
 namespace vtkm
 {
 namespace filter
 {
-/// \brief Estimate the density of particles using the Cloud-in-Cell method
+
-/// This filter treats the CoordinateSystem of a DataSet as positions of particles.
+VTKM_DEPRECATED(1.8,
-/// The particles are infinitesimal in size with finite mass (or other scalar attributes
+                "Use vtkm/filter/density_estimate/ParticleDensityCloudInCell.h instead of "
-/// such as charge). The filter estimates density by imposing a regular grid as
+                "vtkm/filter/ParticleDensityCloudInCell.h.")
-/// specified in the constructor. It spreads the mass of each particle to its 8 nearest
+inline void ParticleDensityCloudInCell_deprecated() {}
-/// neighboring grid points and summing the contribution of particles for each point
+
-/// in the grid.
+inline void ParticleDensityCloudInCell_deprecated_warning()
-/// 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 DivideByVolume 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:
+  ParticleDensityCloudInCell_deprecated();
-  using Superclass = ParticleDensityBase<ParticleDensityCloudInCell>;
+}
 
-  ParticleDensityCloudInCell(const vtkm::Id3& dimension,
+class VTKM_DEPRECATED(1.8, "Use vtkm::filter::density_estimate::ParticleDensityCloudInCell.")
-                             const vtkm::Vec3f& origin,
+  ParticleDensityCloudInCell : public vtkm::filter::density_estimate::ParticleDensityCloudInCell
-                             const vtkm::Vec3f& spacing);
+{
-
+  using density_estimate::ParticleDensityCloudInCell::ParticleDensityCloudInCell;
-  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);
 };
-} // filter
-} // vtkm
 
-#include <vtkm/filter/ParticleDensityCloudInCell.hxx>
+}
-#endif // vtk_m_filter_particle_density_cic_h
+} // namespace vtkm::filter
+
+#endif //vtk_m_filter_ParticleDensityCloudInCell_h

vtkm/filter/ParticleDensityNearestGridPoint.h

@@ -7,51 +7,35 @@
 //  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
 //  PURPOSE.  See the above copyright notice for more information.
 //============================================================================
+#ifndef vtk_m_filter_ParticleDensityNearestGridPoint_h
+#define vtk_m_filter_ParticleDensityNearestGridPoint_h
 
-#ifndef vtk_m_filter_particle_density_ngp_h
+#include <vtkm/Deprecated.h>
-#define vtk_m_filter_particle_density_ngp_h
+#include <vtkm/filter/density_estimate/ParticleDensityNearestGridPoint.h>
-
-#include <vtkm/filter/ParticleDensityBase.h>
 
 namespace vtkm
 {
 namespace filter
 {
-/// \brief Estimate the density of particles using the Nearest Grid Point method
+
-/// This filter treats the CoordinateSystem of a DataSet as positions of particles.
+VTKM_DEPRECATED(1.8,
-/// The particles are infinitesimal in size with finite mass (or other scalar attributes
+                "Use vtkm/filter/density_estimate/ParticleDensityNearestGridPoint.h instead of "
-/// such as charge). The filter estimates density by imposing a regular grid as
+                "vtkm/filter/ParticleDensityNearestGridPoint.h.")
-/// specified in the constructor and summing the mass of particles within each cell
+inline void ParticleDensityNearestGridPoint_deprecated() {}
-/// in the grid.
+
-/// The mass of particles is established by setting the active field (using SetActiveField).
+inline void ParticleDensityNearestGridPoint_deprecated_warning()
-/// 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 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.
-/// In addition, you can also simply count the number of particles in each cell by calling
-/// SetComputeNumberDensity(true).
-class ParticleDensityNearestGridPoint : public ParticleDensityBase<ParticleDensityNearestGridPoint>
 {
-public:
+  ParticleDensityNearestGridPoint_deprecated();
-  using Superclass = ParticleDensityBase<ParticleDensityNearestGridPoint>;
+}
 
-  ParticleDensityNearestGridPoint(const vtkm::Id3& dimension,
+class VTKM_DEPRECATED(1.8, "Use vtkm::filter::density_estimate::ParticleDensityNearestGridPoint.")
-                                  const vtkm::Vec3f& origin,
+  ParticleDensityNearestGridPoint
-                                  const vtkm::Vec3f& spacing);
+  : public vtkm::filter::density_estimate::ParticleDensityNearestGridPoint
-
+{
-  ParticleDensityNearestGridPoint(const vtkm::Id3& dimension, const vtkm::Bounds& bounds);
+  using density_estimate::ParticleDensityNearestGridPoint::ParticleDensityNearestGridPoint;
-
-  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);
 };
-}
-}
 
-#include <vtkm/filter/ParticleDensityNearestGridPoint.hxx>
+}
+} // namespace vtkm::filter
 
-#endif //vtk_m_filter_particle_density_ngp_h
+#endif //vtk_m_filter_ParticleDensityNearestGridPoint_h

vtkm/filter/density_estimate/CMakeLists.txt

@@ -12,13 +12,20 @@ set(density_estimate_headers
   Histogram.h
   NDEntropy.h
   NDHistogram.h
+  ParticleDensityBase.h
+  ParticleDensityCloudInCell.h
+  ParticleDensityNearestGridPoint.h
   )
 
 set(density_estimate_sources_device
   Entropy.cxx
   Histogram.cxx
   NDEntropy.cxx
-  NDHistogram.cxx)
+  NDHistogram.cxx
+  ParticleDensityBase.cxx
+  ParticleDensityCloudInCell.cxx
+  ParticleDensityNearestGridPoint.cxx
+  )
 
 vtkm_library(
   NAME vtkm_filter_density_estimate

vtkm/filter/density_estimate/ParticleDensityBase.cxx

@@ -0,0 +1,54 @@
+//============================================================================
+//  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.
+//============================================================================
+
+#include <vtkm/filter/density_estimate/ParticleDensityBase.h>
+#include <vtkm/worklet/WorkletMapField.h>
+
+namespace
+{
+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
+}
+
+namespace vtkm
+{
+namespace filter
+{
+namespace density_estimate
+{
+
+VTKM_CONT void ParticleDensityBase::DoDivideByVolume(
+  const vtkm::cont::UnknownArrayHandle& density) const
+{
+  auto volume = this->Spacing[0] * this->Spacing[1] * this->Spacing[2];
+  this->Invoke(DivideByVolumeWorklet{ volume }, density);
+}
+} // namespace density_estimate
+} // namespace filter
+} // namespace vtkm

vtkm/filter/density_estimate/ParticleDensityBase.h

@@ -0,0 +1,75 @@
+//============================================================================
+//  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_density_estimate_ParticleDensityBase_h
+#define vtk_m_filter_density_estimate_ParticleDensityBase_h
+
+#include <vtkm/filter/NewFilterField.h>
+#include <vtkm/filter/density_estimate/vtkm_filter_density_estimate_export.h>
+
+namespace vtkm
+{
+namespace filter
+{
+namespace density_estimate
+{
+class VTKM_FILTER_DENSITY_ESTIMATE_EXPORT ParticleDensityBase : public vtkm::filter::NewFilterField
+{
+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:
+  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:
+  // Note: we are using the paradoxical "const ArrayHandle&" parameter whose content can actually
+  // be change by the function.
+  VTKM_CONT void DoDivideByVolume(const vtkm::cont::UnknownArrayHandle& array) const;
+
+  vtkm::Id3 Dimension; // Cell dimension
+  vtkm::Vec3f Origin;
+  vtkm::Vec3f Spacing;
+  bool ComputeNumberDensity;
+  bool DivideByVolume;
+};
+} // namespace density_estimate
+} // namespace filter
+} // namespace vtkm
+
+#endif //vtk_m_filter_density_estimate_ParticleDensityBase_h

vtkm/filter/density_estimate/ParticleDensityCloudInCell.cxx

@@ -8,13 +8,10 @@
 //  PURPOSE.  See the above copyright notice for more information.
 //============================================================================
 
-#ifndef vtk_m_filter_particle_density_cic_hxx
-#define vtk_m_filter_particle_density_cic_hxx
-
 #include <vtkm/cont/ArrayCopy.h>
 #include <vtkm/cont/CellLocatorUniformGrid.h>
 #include <vtkm/cont/DataSetBuilderUniform.h>
-#include <vtkm/filter/PolicyBase.h>
+#include <vtkm/filter/density_estimate/ParticleDensityCloudInCell.h>
 #include <vtkm/worklet/WorkletMapField.h>
 
 namespace vtkm
@@ -74,25 +71,22 @@ namespace vtkm
 {
 namespace filter
 {
-inline VTKM_CONT ParticleDensityCloudInCell::ParticleDensityCloudInCell(const vtkm::Id3& dimension,
+namespace density_estimate
+{
+VTKM_CONT ParticleDensityCloudInCell::ParticleDensityCloudInCell(const vtkm::Id3& dimension,
                                                                  const vtkm::Vec3f& origin,
                                                                  const vtkm::Vec3f& spacing)
   : Superclass(dimension, origin, spacing)
 {
 }
 
-inline VTKM_CONT ParticleDensityCloudInCell::ParticleDensityCloudInCell(const Id3& dimension,
+VTKM_CONT ParticleDensityCloudInCell::ParticleDensityCloudInCell(const Id3& dimension,
                                                                  const vtkm::Bounds& bounds)
   : Superclass(dimension, bounds)
 {
 }
 
-template <typename T, typename StorageType, typename Policy>
+VTKM_CONT vtkm::cont::DataSet ParticleDensityCloudInCell::DoExecute(const cont::DataSet& input)
-inline VTKM_CONT vtkm::cont::DataSet ParticleDensityCloudInCell::DoExecute(
-  const cont::DataSet& dataSet,
-  const cont::ArrayHandle<T, StorageType>& field,
-  const vtkm::filter::FieldMetadata&,
-  PolicyBase<Policy>)
 {
   // Unlike ParticleDensityNGP, particle deposit mass on the grid points, thus it is natural to
   // return the density as PointField;
@@ -104,30 +98,56 @@ inline VTKM_CONT vtkm::cont::DataSet ParticleDensityCloudInCell::DoExecute(
   locator.SetCoordinates(uniform.GetCoordinateSystem());
   locator.Update();
 
-  auto coords = dataSet.GetCoordinateSystem().GetDataAsMultiplexer();
+  auto coords = input.GetCoordinateSystem().GetDataAsMultiplexer();
 
+  auto resolveType = [&, this](const auto& concrete) {
+    // use std::decay to remove const ref from the decltype of concrete.
+    using T = typename std::decay_t<decltype(concrete)>::ValueType;
+
+    // We create an ArrayHandle and pass it to the Worklet as AtomicArrayInOut.
+    // However, the ArrayHandle needs to be allocated and initialized first. The
+    // easiest way to do it is to copy from an ArrayHandleConstant
     vtkm::cont::ArrayHandle<T> density;
     vtkm::cont::ArrayCopy(vtkm::cont::ArrayHandleConstant<T>(0, uniform.GetNumberOfPoints()),
                           density);
 
     this->Invoke(vtkm::worklet::CICWorklet{},
                  coords,
-               field,
+                 concrete,
                  locator,
                  uniform.GetCellSet().template AsCellSet<vtkm::cont::CellSetStructured<3>>(),
                  density);
 
     if (DivideByVolume)
     {
-    auto volume = this->Spacing[0] * this->Spacing[1] * this->Spacing[2];
+      this->DoDivideByVolume(density);
-    this->Invoke(DivideByVolumeWorklet{ volume }, density);
     }
 
     uniform.AddField(vtkm::cont::make_FieldPoint("density", density));
+  };
+
+  // Note: This is the so called Immediately-Invoked Function Expression (IIFE). Here we define
+  // a lambda expression and immediately call it at the end. This allows us to not declare an
+  // UnknownArrayHandle first and then assign it in the if-else statement. If I really want to
+  // show-off, I can even inline the `fieldArray` variable and turn it into a long expression.
+  auto fieldArray = [&, this]() -> vtkm::cont::UnknownArrayHandle {
+    if (this->ComputeNumberDensity)
+    {
+      return vtkm::cont::make_ArrayHandleConstant(vtkm::FloatDefault{ 1 },
+                                                  input.GetNumberOfPoints());
+    }
+    else
+    {
+      return this->GetFieldFromDataSet(input).GetData();
+    }
+  }();
+  fieldArray.CastAndCallForTypes<
+    vtkm::TypeListFieldScalar,
+    vtkm::ListAppend<VTKM_DEFAULT_STORAGE_LIST, vtkm::List<vtkm::cont::StorageTagConstant>>>(
+    resolveType);
 
   return uniform;
 }
-
+} // namespace density_estimate
-}
+} // namespace filter
-}
+} // namespace vtkm
-#endif // vtk_m_filter_particle_density_cic_hxx

vtkm/filter/density_estimate/ParticleDensityCloudInCell.h

@@ -0,0 +1,56 @@
+//============================================================================
+//  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_density_estimate_ParticleDensityCIC_h
+#define vtk_m_filter_density_estimate_ParticleDensityCIC_h
+
+#include <vtkm/filter/density_estimate/ParticleDensityBase.h>
+
+namespace vtkm
+{
+namespace filter
+{
+namespace density_estimate
+{
+/// \brief Estimate the density of particles using the Cloud-in-Cell method
+/// 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 spreads 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 DivideByVolume feature of this filter.
+/// In addition, you can also simply count the number of particles in each cell by calling
+/// SetComputeNumberDensity(true).
+class VTKM_FILTER_DENSITY_ESTIMATE_EXPORT ParticleDensityCloudInCell : public ParticleDensityBase
+{
+public:
+  using Superclass = ParticleDensityBase;
+
+  ParticleDensityCloudInCell(const vtkm::Id3& dimension,
+                             const vtkm::Vec3f& origin,
+                             const vtkm::Vec3f& spacing);
+
+  ParticleDensityCloudInCell(const Id3& dimension, const vtkm::Bounds& bounds);
+
+private:
+  VTKM_CONT vtkm::cont::DataSet DoExecute(const vtkm::cont::DataSet& input) override;
+};
+} // namespace density_estimate
+} // namespace filter
+} // namespace vtkm
+
+#endif // vtk_m_filter_density_estimate_ParticleDensityCIC_h

vtkm/filter/density_estimate/ParticleDensityNearestGridPoint.cxx

@@ -8,14 +8,11 @@
 //  PURPOSE.  See the above copyright notice for more information.
 //============================================================================
 
-#ifndef vtk_m_filter_particle_density_ngp_hxx
-#define vtk_m_filter_particle_density_ngp_hxx
-
 #include <vtkm/cont/ArrayCopy.h>
 #include <vtkm/cont/ArrayHandleConstant.h>
 #include <vtkm/cont/CellLocatorUniformGrid.h>
 #include <vtkm/cont/DataSetBuilderUniform.h>
-#include <vtkm/filter/PolicyBase.h>
+#include <vtkm/filter/density_estimate/ParticleDensityNearestGridPoint.h>
 #include <vtkm/worklet/WorkletMapField.h>
 
 namespace vtkm
@@ -57,7 +54,9 @@ namespace vtkm
 {
 namespace filter
 {
-inline VTKM_CONT ParticleDensityNearestGridPoint::ParticleDensityNearestGridPoint(
+namespace density_estimate
+{
+VTKM_CONT ParticleDensityNearestGridPoint::ParticleDensityNearestGridPoint(
   const vtkm::Id3& dimension,
   const vtkm::Vec3f& origin,
   const vtkm::Vec3f& spacing)
@@ -65,20 +64,15 @@ inline VTKM_CONT ParticleDensityNearestGridPoint::ParticleDensityNearestGridPoin
 {
 }
 
-inline VTKM_CONT ParticleDensityNearestGridPoint::ParticleDensityNearestGridPoint(
+VTKM_CONT ParticleDensityNearestGridPoint::ParticleDensityNearestGridPoint(
   const Id3& dimension,
   const vtkm::Bounds& bounds)
   : Superclass(dimension, bounds)
 {
 }
 
-template <typename T, typename StorageType, typename Policy>
+VTKM_CONT vtkm::cont::DataSet ParticleDensityNearestGridPoint::DoExecute(
-inline VTKM_CONT vtkm::cont::DataSet ParticleDensityNearestGridPoint::DoExecute(
+  const vtkm::cont::DataSet& input)
-  const vtkm::cont::DataSet& dataSet,
-  const vtkm::cont::ArrayHandle<T, StorageType>&
-    field, // particles' scala field to be deposited to the mesh, e.g. mass or charge
-  const vtkm::filter::FieldMetadata&,
-  vtkm::filter::PolicyBase<Policy>)
 {
   // TODO: it really doesn't need to be a UniformGrid, any CellSet with CellLocator will work.
   //  Make it another input rather an output generated.
@@ -96,27 +90,51 @@ inline VTKM_CONT vtkm::cont::DataSet ParticleDensityNearestGridPoint::DoExecute(
   locator.SetCoordinates(uniform.GetCoordinateSystem());
   locator.Update();
 
-  auto coords = dataSet.GetCoordinateSystem().GetDataAsMultiplexer();
+  auto coords = input.GetCoordinateSystem().GetDataAsMultiplexer();
+
+  auto resolveType = [&, this](const auto& concrete) {
+    // use std::decay to remove const ref from the decltype of concrete.
+    using T = typename std::decay_t<decltype(concrete)>::ValueType;
 
     // We create an ArrayHandle and pass it to the Worklet as AtomicArrayInOut.
-  // However the ArrayHandle needs to be allocated and initialized first. The
+    // However, the ArrayHandle needs to be allocated and initialized first. The
     // easiest way to do it is to copy from an ArrayHandleConstant
     vtkm::cont::ArrayHandle<T> density;
-  vtkm::cont::ArrayCopy(vtkm::cont::ArrayHandleConstant<T>(0, uniform.GetNumberOfCells()), density);
+    vtkm::cont::ArrayCopy(vtkm::cont::ArrayHandleConstant<T>(0, uniform.GetNumberOfCells()),
+                          density);
 
-  this->Invoke(vtkm::worklet::NGPWorklet{}, coords, field, locator, density);
+    this->Invoke(vtkm::worklet::NGPWorklet{}, coords, concrete, locator, density);
 
     if (DivideByVolume)
     {
-    auto volume = this->Spacing[0] * this->Spacing[1] * this->Spacing[2];
+      this->DoDivideByVolume(density);
-    this->Invoke(DivideByVolumeWorklet{ volume }, density);
     }
 
     uniform.AddField(vtkm::cont::make_FieldCell("density", density));
+  };
+
+  // Note: This is the so called Immediately-Invoked Function Expression (IIFE). Here we define
+  // a lambda expression and immediately call it at the end. This allows us to not declare an
+  // UnknownArrayHandle first and then assign it in the if-else statement. If I really want to
+  // show-off, I can even inline the `fieldArray` variable and turn it into a long expression.
+  auto fieldArray = [&, this]() -> vtkm::cont::UnknownArrayHandle {
+    if (this->ComputeNumberDensity)
+    {
+      return vtkm::cont::make_ArrayHandleConstant(vtkm::FloatDefault{ 1 },
+                                                  input.GetNumberOfPoints());
+    }
+    else
+    {
+      return this->GetFieldFromDataSet(input).GetData();
+    }
+  }();
+  fieldArray.CastAndCallForTypes<
+    vtkm::TypeListFieldScalar,
+    vtkm::ListAppend<VTKM_DEFAULT_STORAGE_LIST, vtkm::List<vtkm::cont::StorageTagConstant>>>(
+    resolveType);
 
   return uniform;
 }
-
 }
 }
-#endif //vtk_m_filter_particle_density_ngp_hxx
+}

vtkm/filter/density_estimate/ParticleDensityNearestGridPoint.h

@@ -0,0 +1,55 @@
+//============================================================================
+//  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_density_estimate_ParticleDensityNGP_h
+#define vtk_m_filter_density_estimate_ParticleDensityNGP_h
+
+#include <vtkm/filter/density_estimate/ParticleDensityBase.h>
+
+namespace vtkm
+{
+namespace filter
+{
+namespace density_estimate
+{
+/// \brief Estimate the density of particles using the Nearest Grid Point method
+/// 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 and summing the mass of particles within each cell
+/// 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 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.
+/// In addition, you can also simply count the number of particles in each cell by calling
+/// SetComputeNumberDensity(true).
+class VTKM_FILTER_DENSITY_ESTIMATE_EXPORT ParticleDensityNearestGridPoint
+  : public ParticleDensityBase
+{
+public:
+  using Superclass = ParticleDensityBase;
+
+  ParticleDensityNearestGridPoint(const vtkm::Id3& dimension,
+                                  const vtkm::Vec3f& origin,
+                                  const vtkm::Vec3f& spacing);
+
+  ParticleDensityNearestGridPoint(const vtkm::Id3& dimension, const vtkm::Bounds& bounds);
+
+private:
+  VTKM_CONT vtkm::cont::DataSet DoExecute(const vtkm::cont::DataSet& input) override;
+};
+} // namespace density_estimate
+} // namespace filter
+} // namespace vtkm
+#endif //vtk_m_filter_density_estimate_ParticleDensityNGP_h

vtkm/filter/density_estimate/testing/CMakeLists.txt

@@ -13,7 +13,8 @@ set(unit_tests
   UnitTestHistogramFilter.cxx
   UnitTestNDEntropyFilter.cxx
   UnitTestNDHistogramFilter.cxx
-  UnitTestPartitionedDataSetHistogramFilter.cxx)
+  UnitTestPartitionedDataSetHistogramFilter.cxx
+  UnitTestParticleDensity.cxx)
 
 set(libraries
   vtkm_filter_density_estimate

vtkm/filter/density_estimate/testing/UnitTestParticleDensity.cxx

@@ -11,8 +11,8 @@
 #include <vtkm/cont/ArrayHandleRandomUniformReal.h>
 #include <vtkm/cont/DataSetBuilderExplicit.h>
 #include <vtkm/cont/testing/Testing.h>
-#include <vtkm/filter/ParticleDensityCloudInCell.h>
+#include <vtkm/filter/density_estimate/ParticleDensityCloudInCell.h>
-#include <vtkm/filter/ParticleDensityNearestGridPoint.h>
+#include <vtkm/filter/density_estimate/ParticleDensityNearestGridPoint.h>
 #include <vtkm/worklet/DescriptiveStatistics.h>
 
 void TestNGP()
@@ -41,7 +41,7 @@ void TestNGP()
   dataSet.AddCellField("mass", mass);
 
   auto cellDims = vtkm::Id3{ 3, 3, 3 };
-  vtkm::filter::ParticleDensityNearestGridPoint filter{
+  vtkm::filter::density_estimate::ParticleDensityNearestGridPoint filter{
     cellDims, { 0.f, 0.f, 0.f }, vtkm::Vec3f{ 1.f / 3.f, 1.f / 3.f, 1.f / 3.f }
   };
   filter.SetActiveField("mass");
@@ -91,9 +91,9 @@ void TestCIC()
   dataSet.AddCellField("mass", mass);
 
   auto cellDims = vtkm::Id3{ 3, 3, 3 };
-  vtkm::filter::ParticleDensityCloudInCell filter{ cellDims,
+  vtkm::filter::density_estimate::ParticleDensityCloudInCell filter{
-                                                   { 0.f, 0.f, 0.f },
+    cellDims, { 0.f, 0.f, 0.f }, vtkm::Vec3f{ 1.f / 3.f, 1.f / 3.f, 1.f / 3.f }
-                                                   vtkm::Vec3f{ 1.f / 3.f, 1.f / 3.f, 1.f / 3.f } };
+  };
   filter.SetActiveField("mass");
   auto density = filter.Execute(dataSet);
 

vtkm/filter/testing/CMakeLists.txt

@@ -45,7 +45,6 @@ set(unit_tests
   UnitTestMeshQualityFilter.cxx
   UnitTestMIRFilter.cxx
   UnitTestMultiBlockFilter.cxx
-  UnitTestParticleDensity.cxx
   UnitTestPartitionedDataSetFilters.cxx
   UnitTestPointAverageFilter.cxx
   UnitTestPointAverageCellSetExtrude.cxx