Merge branch 'master' into particle_density_cic

.gitlab-ci.yml

@@ -131,6 +131,7 @@ stages:
     - export PATH=$PWD/.gitlab:$PATH
     - SCCACHE_IDLE_TIMEOUT=0 sccache --start-server
     - sccache --show-stats
+    - .gitlab/ci/config/google_benchmarks.sh
     - "cmake --version"
     - "cmake -V -P .gitlab/ci/config/gitlab_ci_setup.cmake"
     - "ctest -VV -S .gitlab/ci/ctest_configure.cmake"

.gitlab/ci/config/google_benchmarks.sh

@@ -0,0 +1,27 @@
+#!/bin/bash
+
+set -xe
+
+readonly version="v1.5.2"
+readonly tarball="$version.tar.gz"
+readonly url="https://github.com/google/benchmark/archive/$tarball"
+readonly sha256sum="dccbdab796baa1043f04982147e67bb6e118fe610da2c65f88912d73987e700c"
+readonly install_dir="$HOME/gbench"
+
+if ! [[ "$VTKM_SETTINGS" =~ "benchmarks" ]]; then
+  exit 0
+fi
+
+cd "$HOME"
+
+echo "$sha256sum  $tarball" > gbenchs.sha256sum
+curl --insecure -OL "$url"
+sha256sum --check gbenchs.sha256sum
+tar xf "$tarball"
+
+mkdir build
+mkdir "$install_dir"
+
+cmake -GNinja -S benchmark* -B build -DBENCHMARK_DOWNLOAD_DEPENDENCIES=ON
+cmake --build build
+cmake --install build --prefix "$install_dir"

.gitlab/ci/config/initial_config.cmake

@@ -64,6 +64,7 @@ foreach(option IN LISTS options)
 
   elseif(benchmarks STREQUAL option)
     set(VTKm_ENABLE_BENCHMARKS "ON" CACHE STRING "")
+    set(ENV{CMAKE_PREFIX_PATH} "$ENV{HOME}/gbench")
 
   elseif(mpi STREQUAL option)
     set(VTKm_ENABLE_MPI "ON" CACHE STRING "")

.gitlab/ci/ubuntu1804.yml

@@ -16,7 +16,7 @@ build:ubuntu1804_gcc9:
     CC: "gcc-9"
     CXX: "g++-9"
     CMAKE_BUILD_TYPE: Debug
-    VTKM_SETTINGS: "tbb+openmp+mpi+shared+hdf5"
+    VTKM_SETTINGS: "benchmarks+tbb+openmp+mpi+shared+hdf5"
 
 test:ubuntu1804_gcc9:
   tags:
@@ -56,7 +56,7 @@ build:ubuntu1804_gcc7:
     CC: "gcc-7"
     CXX: "g++-7"
     CUDAHOSTCXX: "g++-7"
-    VTKM_SETTINGS: "cuda+turing+mpi+64bit_floats+no_virtual"
+    VTKM_SETTINGS: "benchmarks+cuda+turing+mpi+64bit_floats+no_virtual"
 
 test:ubuntu1804_gcc7:
   tags:
@@ -202,7 +202,7 @@ build:ubuntu1804_kokkos:
   variables:
     CMAKE_GENERATOR: "Ninja"
     CMAKE_BUILD_TYPE: Release
-    VTKM_SETTINGS: "kokkos+turing+static+64bit_floats"
+    VTKM_SETTINGS: "benchmarks+kokkos+turing+static+64bit_floats"
 
 test:ubuntu1804_kokkos:
   tags:

vtkm/cont/CellSetExtrude.cxx

@@ -182,6 +182,31 @@ void CellSetExtrude::GetCellPointIds(vtkm::Id id, vtkm::Id* ptids) const
   }
 }
 
+template <vtkm::IdComponent NumIndices>
+VTKM_CONT void CellSetExtrude::GetIndices(vtkm::Id index,
+                                          vtkm::Vec<vtkm::Id, NumIndices>& ids) const
+{
+  static_assert(NumIndices == 6, "There are always 6 points in a wedge.");
+  this->GetCellPointIds(index, ids.data());
+}
+
+VTKM_CONT void CellSetExtrude::GetIndices(vtkm::Id index,
+                                          vtkm::cont::ArrayHandle<vtkm::Id>& ids) const
+{
+  ids.Allocate(6);
+  auto outIdPortal = ids.WritePortal();
+  vtkm::cont::Token token;
+  auto conn = this->PrepareForInput(vtkm::cont::DeviceAdapterTagSerial{},
+                                    vtkm::TopologyElementTagCell{},
+                                    vtkm::TopologyElementTagPoint{},
+                                    token);
+  auto indices = conn.GetIndices(index);
+  for (vtkm::IdComponent i = 0; i < 6; i++)
+  {
+    outIdPortal.Set(i, indices[i]);
+  }
+}
+
 std::shared_ptr<CellSet> CellSetExtrude::NewInstance() const
 {
   return std::make_shared<CellSetExtrude>();

vtkm/cont/CellSetExtrude.h

@@ -100,6 +100,11 @@ public:
 
   bool GetIsPeriodic() const { return this->IsPeriodic; }
 
+  template <vtkm::IdComponent NumIndices>
+  VTKM_CONT void GetIndices(vtkm::Id index, vtkm::Vec<vtkm::Id, NumIndices>& ids) const;
+
+  VTKM_CONT void GetIndices(vtkm::Id index, vtkm::cont::ArrayHandle<vtkm::Id>& ids) const;
+
   template <typename VisitTopology, typename IncidentTopology>
   using ExecConnectivityType =
     typename detail::CellSetExtrudeConnectivityChooser<VisitTopology,

vtkm/filter/ParticleDensityNearestGridPoint.h

@@ -18,6 +18,20 @@ 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.
+/// 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).
 
 // We only need the CoordinateSystem and scalar fields of the input dataset thus a FilterField
 class ParticleDensityNearestGridPoint
@@ -33,16 +47,30 @@ 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

@@ -12,6 +12,7 @@
 #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>
@@ -49,6 +50,30 @@ 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
 
@@ -64,6 +89,8 @@ inline VTKM_CONT ParticleDensityNearestGridPoint::ParticleDensityNearestGridPoin
   : Dimension(dimension)
   , Origin(origin)
   , Spacing(spacing)
+  , ComputeNumberDensity(false)
+  , DivideByVolume(true)
 {
 }
 
@@ -78,9 +105,30 @@ inline VTKM_CONT ParticleDensityNearestGridPoint::ParticleDensityNearestGridPoin
                          static_cast<vtkm::FloatDefault>(bounds.Y.Length()),
                          static_cast<vtkm::FloatDefault>(bounds.Z.Length()) } /
             Dimension)
+  , ComputeNumberDensity(false)
+  , DivideByVolume(true)
 {
 }
 
+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,
@@ -105,17 +153,22 @@ inline VTKM_CONT vtkm::cont::DataSet ParticleDensityNearestGridPoint::DoExecute(
   locator.SetCoordinates(uniform.GetCoordinateSystem());
   locator.Update();
 
-  vtkm::cont::ArrayHandle<vtkm::Vec3f> coords;
-  dataSet.GetCoordinateSystem().GetData().AsArrayHandle<vtkm::Vec3f>(coords);
+  auto coords = dataSet.GetCoordinateSystem().GetDataAsMultiplexer();
 
   // We create an ArrayHandle and pass it to the Worklet as AtomicArrayInOut.
   // However the ArrayHandle needs to be allocated and initialized first. The
-  // easily way to do it is to copy from an ArrayHandleConstant
+  // 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);
 
   this->Invoke(vtkm::worklet::NGPWorklet{}, coords, field, locator, density);
 
+  if (DivideByVolume)
+  {
+    auto volume = this->Spacing[0] * this->Spacing[1] * this->Spacing[2];
+    this->Invoke(vtkm::worklet::detail::DividByVolume{ volume }, density);
+  }
+
   uniform.AddField(vtkm::cont::make_FieldCell("density", density));
 
   return uniform;

vtkm/filter/testing/UnitTestParticleDensity.cxx

@@ -35,8 +35,9 @@ void TestNGP()
   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);
+  vtkm::cont::ArrayHandle<vtkm::FloatDefault> mass;
+  vtkm::cont::ArrayCopy(vtkm::cont::ArrayHandleRandomUniformReal<vtkm::FloatDefault>(N, 0xd1ce),
+                        mass);
   dataSet.AddCellField("mass", mass);
 
   auto cellDims = vtkm::Id3{ 3, 3, 3 };
@@ -46,13 +47,23 @@ void TestNGP()
   filter.SetActiveField("mass");
   auto density = filter.Execute(dataSet);
 
-  vtkm::cont::ArrayHandle<vtkm::Float32> field;
-  density.GetCellField("density").GetData().AsArrayHandle<vtkm::Float32>(field);
+  vtkm::cont::ArrayHandle<vtkm::FloatDefault> field;
+  density.GetCellField("density").GetData().AsArrayHandle<vtkm::FloatDefault>(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(), 10));
+  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.GetCellField("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 TestCIC()

vtkm/io/VTKDataSetWriter.cxx

@@ -396,6 +396,10 @@ void Write(std::ostream& out, const vtkm::cont::DataSet& dataSet)
     // these function just like explicit cell sets
     WriteDataSetAsUnstructured(out, dataSet, cellSet.Cast<vtkm::cont::CellSetSingleType<>>());
   }
+  else if (cellSet.IsType<vtkm::cont::CellSetExtrude>())
+  {
+    WriteDataSetAsUnstructured(out, dataSet, cellSet.Cast<vtkm::cont::CellSetExtrude>());
+  }
   else
   {
     throw vtkm::cont::ErrorBadType("Could not determine type to write out.");

vtkm/rendering/Camera.cxx

@@ -286,7 +286,7 @@ void Camera::ResetToBounds(const vtkm::Bounds& dataBounds,
   vtkm::Float32 diagonalLength = vtkm::Magnitude(totalExtent);
   this->SetPosition(center + directionOfProjection * diagonalLength * 1.0f);
   this->SetFieldOfView(60.0f);
-  this->SetClippingRange(0.1f * diagonalLength, diagonalLength * 5.0f);
+  this->SetClippingRange(0.1f * diagonalLength, diagonalLength * 10.0f);
 
   // Reset for 2D camera
   this->SetViewRange2D(db);