Merge topic 'particle_density_cic'

508b992cb Merge branch 'particle_density_cic' of gitlab.kitware.com:ollielo/vtk-m into particle_density_cic
056ee7269 fixed another copy paste bug
f71a01b39 Fixed misunderstanding of tolerance.
afe0a3e5f correct uncareful search and replace error
b3c6ea396 install ParticleDensityBase.h
7b8cc401c public CUDA
6796b1a45 consolidate two particle density filters
4a22f54bd Merge branch 'master' into particle_density_cic
...

Acked-by: Kitware Robot <kwrobot@kitware.com>
Acked-by: Kenneth Moreland <kmorel@acm.org>
Merge-request: !2434

vtkm/filter/CMakeLists.txt

@@ -83,6 +83,8 @@ set(extra_headers
   MeshQuality.h
   NDEntropy.h
   NDHistogram.h
+  ParticleDensityBase.h
+  ParticleDensityCloudInCell.h
   ParticleDensityNearestGridPoint.h
   ParticleAdvection.h
   Pathline.h
@@ -135,6 +137,7 @@ set(extra_header_template_sources
   MeshQuality.hxx
   NDEntropy.hxx
   NDHistogram.hxx
+  ParticleDensityCloudInCell.hxx
   ParticleDensityNearestGridPoint.hxx
   ParticleAdvection.hxx
   Pathline.hxx

vtkm/filter/ParticleDensityBase.h

@@ -0,0 +1,115 @@
+//============================================================================
+//  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;
+
+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

@@ -0,0 +1,57 @@
+//============================================================================
+//  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_cic_h
+#define vtk_m_filter_particle_density_cic_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.
+/// 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 ParticleDensityCloudInCell : public ParticleDensityBase<ParticleDensityCloudInCell>
+{
+public:
+  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);
+};
+} // filter
+} // vtkm
+
+#include <vtkm/filter/ParticleDensityCloudInCell.hxx>
+#endif // vtk_m_filter_particle_density_cic_h

vtkm/filter/ParticleDensityCloudInCell.hxx

@@ -0,0 +1,129 @@
+//============================================================================
+//  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_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/worklet/WorkletMapField.h>
+
+namespace vtkm
+{
+namespace worklet
+{
+class CICWorklet : public vtkm::worklet::WorkletMapField
+{
+public:
+  using ControlSignature = void(FieldIn coords,
+                                FieldIn field,
+                                ExecObject locator,
+                                WholeCellSetIn<Cell, Point> cellSet,
+                                AtomicArrayInOut density);
+  using ExecutionSignature = void(_1, _2, _3, _4, _5);
+
+  template <typename Point,
+            typename T,
+            typename CellLocatorExecObj,
+            typename CellSet,
+            typename AtomicArray>
+  VTKM_EXEC void operator()(const Point& point,
+                            const T value,
+                            const CellLocatorExecObj& locator,
+                            const CellSet& cellSet,
+                            AtomicArray& density) const
+  {
+    vtkm::Id cellId{};
+    vtkm::Vec3f parametric;
+
+    if (locator.FindCell(point, cellId, parametric) == vtkm::ErrorCode::Success)
+    {
+      // iterate through all the points of the cell and deposit with correct weight.
+      auto indices = cellSet.GetIndices(cellId);
+      auto rparametric = vtkm::Vec3f{ 1, 1, 1 } - parametric;
+
+      // deposit the scalar field value in proportion to the volume of the sub-hexahedron
+      // the vertex is in.
+      density.Add(indices[0], value * parametric[0] * parametric[1] * parametric[2]);
+      density.Add(indices[1], value * rparametric[0] * parametric[1] * parametric[2]);
+      density.Add(indices[2], value * rparametric[0] * rparametric[1] * parametric[2]);
+      density.Add(indices[3], value * parametric[0] * rparametric[1] * parametric[2]);
+
+      density.Add(indices[4], value * parametric[0] * parametric[1] * rparametric[2]);
+      density.Add(indices[5], value * rparametric[0] * parametric[1] * rparametric[2]);
+      density.Add(indices[6], value * rparametric[0] * rparametric[1] * rparametric[2]);
+      density.Add(indices[7], value * parametric[0] * rparametric[1] * rparametric[2]);
+    }
+
+    // We simply ignore that particular particle when it is not in the mesh.
+  }
+};
+} // worklet
+} // vtkm
+
+namespace vtkm
+{
+namespace filter
+{
+inline 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,
+                                                                        const vtkm::Bounds& bounds)
+  : Superclass(dimension, bounds)
+{
+}
+
+template <typename T, typename StorageType, typename Policy>
+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;
+  auto uniform = vtkm::cont::DataSetBuilderUniform::Create(
+    this->Dimension + vtkm::Id3{ 1, 1, 1 }, this->Origin, this->Spacing);
+
+  vtkm::cont::CellLocatorUniformGrid locator;
+  locator.SetCellSet(uniform.GetCellSet());
+  locator.SetCoordinates(uniform.GetCoordinateSystem());
+  locator.Update();
+
+  vtkm::cont::ArrayHandle<vtkm::Vec3f> coords;
+  dataSet.GetCoordinateSystem().GetData().AsArrayHandle<vtkm::Vec3f>(coords);
+
+  vtkm::cont::ArrayHandle<T> density;
+  vtkm::cont::ArrayCopy(vtkm::cont::ArrayHandleConstant<T>(0, uniform.GetNumberOfPoints()),
+                        density);
+
+  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;
+}
+
+}
+}
+#endif // vtk_m_filter_particle_density_cic_hxx

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
 {
@@ -32,14 +32,10 @@ namespace filter
 /// 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).
-
-// 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,48 +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)
 {
 }
 
-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)
+inline VTKM_CONT ParticleDensityNearestGridPoint::ParticleDensityNearestGridPoint(
+  const Id3& dimension,
+  const vtkm::Bounds& bounds)
+  : 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,
@@ -165,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

@@ -11,6 +11,7 @@
 #include <vtkm/cont/ArrayHandleRandomUniformReal.h>
 #include <vtkm/cont/DataSetBuilderExplicit.h>
 #include <vtkm/cont/testing/Testing.h>
+#include <vtkm/filter/ParticleDensityCloudInCell.h>
 #include <vtkm/filter/ParticleDensityNearestGridPoint.h>
 #include <vtkm/worklet/DescriptiveStatistics.h>
 
@@ -65,9 +66,60 @@ void TestNGP()
   VTKM_TEST_ASSERT(test_equal(counts_result.Sum(), mass_result.N(), 0.1));
 }
 
+void TestCIC()
+{
+  const vtkm::Id N = 1000000;
+  // This is a better way to create this array, but I am temporarily breaking this
+  // functionality (November 2020) so that I can split up merge requests that move
+  // ArrayHandles to the new buffer types. This should be restored once
+  // ArrayHandleTransform is converted to the new type.
+  auto composite = vtkm::cont::make_ArrayHandleCompositeVector(
+    vtkm::cont::ArrayHandleRandomUniformReal<vtkm::Float32>(N, 0xceed),
+    vtkm::cont::ArrayHandleRandomUniformReal<vtkm::Float32>(N, 0xdeed),
+    vtkm::cont::ArrayHandleRandomUniformReal<vtkm::Float32>(N, 0xabba));
+  vtkm::cont::ArrayHandle<vtkm::Vec3f> positions;
+  vtkm::cont::ArrayCopy(composite, positions);
+
+  vtkm::cont::ArrayHandle<vtkm::Id> connectivity;
+  vtkm::cont::ArrayCopy(vtkm::cont::make_ArrayHandleIndex(N), connectivity);
+
+  auto dataSet = vtkm::cont::DataSetBuilderExplicit::Create(
+    positions, vtkm::CellShapeTagVertex{}, 1, connectivity);
+
+  vtkm::cont::ArrayHandle<vtkm::Float32> mass;
+  vtkm::cont::ArrayCopy(vtkm::cont::ArrayHandleRandomUniformReal<vtkm::Float32>(N, 0xd1ce), mass);
+  dataSet.AddCellField("mass", mass);
+
+  auto cellDims = vtkm::Id3{ 3, 3, 3 };
+  vtkm::filter::ParticleDensityCloudInCell 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");
+  auto density = filter.Execute(dataSet);
+
+  vtkm::cont::ArrayHandle<vtkm::Float32> field;
+  density.GetPointField("density").GetData().AsArrayHandle<vtkm::Float32>(field);
+
+  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() * 27.0, 0.1));
+
+  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()
 {
   TestNGP();
+  TestCIC();
 }
 
 int UnitTestParticleDensity(int argc, char* argv[])