4 files changed, 85 insertions, 0 deletions
diff --git a/Source/Utils/CMakeLists.txt b/Source/Utils/CMakeLists.txt
index ce785cc33..c9da2ae47 100644
--- a/Source/Utils/CMakeLists.txt
+++ b/Source/Utils/CMakeLists.txt
@@ -4,6 +4,7 @@ target_sources(WarpX
     CoarsenMR.cpp
     Interpolate.cpp
     IntervalsParser.cpp
+    ParticleUtils.cpp
     RelativeCellPosition.cpp
     WarpXAlgorithmSelection.cpp
     WarpXMovingWindow.cpp
diff --git a/Source/Utils/Make.package b/Source/Utils/Make.package
index 08b9aa3a9..61ced5d7f 100644
--- a/Source/Utils/Make.package
+++ b/Source/Utils/Make.package
@@ -7,5 +7,6 @@ CEXE_sources += CoarsenMR.cpp
 CEXE_sources += Interpolate.cpp
 CEXE_sources += IntervalsParser.cpp
 CEXE_sources += RelativeCellPosition.cpp
+CEXE_sources += ParticleUtils.cpp
 
 VPATH_LOCATIONS   += $(WARPX_HOME)/Source/Utils
diff --git a/Source/Utils/ParticleUtils.H b/Source/Utils/ParticleUtils.H
new file mode 100644
index 000000000..cb25607a3
--- /dev/null
+++ b/Source/Utils/ParticleUtils.H
@@ -0,0 +1,31 @@
+/* Copyright 2019-2020 Neil Zaim, Yinjian Zhao
+ *
+ * This file is part of WarpX.
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+#ifndef WARPX_PARTICLE_UTILS_H_
+#define WARPX_PARTICLE_UTILS_H_
+
+#include "Particles/WarpXParticleContainer.H"
+#include <AMReX_DenseBins.H>
+
+namespace ParticleUtils {
+
+    /**
+    * \brief Find the particles and count the particles that are in each cell. More specifically
+    * this function returns an amrex::DenseBins object containing an offset array and a permutation
+    * array which can be used to loop over all the cells in a tile and apply an algorithm to
+    * particles of a given species present in each cell.
+    * Note that this does *not* rearrange particle arrays.
+    *
+    * @param[in] lev the index of the refinement level.
+    * @param[in] mfi the MultiFAB iterator.
+    * @param[in] ptile the particle tile.
+    */
+    amrex::DenseBins<WarpXParticleContainer::ParticleType>
+    findParticlesInEachCell( int const lev, amrex::MFIter const& mfi,
+                             WarpXParticleContainer::ParticleTileType const& ptile);
+}
+
+#endif // WARPX_PARTICLE_UTILS_H_
diff --git a/Source/Utils/ParticleUtils.cpp b/Source/Utils/ParticleUtils.cpp
new file mode 100644
index 000000000..92c03bc30
--- /dev/null
+++ b/Source/Utils/ParticleUtils.cpp
@@ -0,0 +1,52 @@
+/* Copyright 2019-2020 Neil Zaim, Yinjian Zhao
+ *
+ * This file is part of WarpX.
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+#include "ParticleUtils.H"
+#include "WarpX.H"
+
+namespace ParticleUtils {
+
+    using namespace amrex;
+    // Define shortcuts for frequently-used type names
+    using ParticleType = WarpXParticleContainer::ParticleType;
+    using ParticleTileType = WarpXParticleContainer::ParticleTileType;
+    using ParticleBins = DenseBins<ParticleType>;
+    using index_type = ParticleBins::index_type;
+
+    /* Find the particles and count the particles that are in each cell.
+       Note that this does *not* rearrange particle arrays */
+    ParticleBins
+    findParticlesInEachCell( int const lev, MFIter const& mfi,
+                             ParticleTileType const& ptile) {
+
+        // Extract particle structures for this tile
+        int const np = ptile.numParticles();
+        ParticleType const* particle_ptr = ptile.GetArrayOfStructs()().data();
+
+        // Extract box properties
+        Geometry const& geom = WarpX::GetInstance().Geom(lev);
+        Box const& cbx = mfi.tilebox(IntVect::TheZeroVector()); //Cell-centered box
+        const auto lo = lbound(cbx);
+        const auto dxi = geom.InvCellSizeArray();
+        const auto plo = geom.ProbLoArray();
+
+        // Find particles that are in each cell;
+        // results are stored in the object `bins`.
+        ParticleBins bins;
+        bins.build(np, particle_ptr, cbx,
+            // Pass lambda function that returns the cell index
+            [=] AMREX_GPU_HOST_DEVICE (const ParticleType& p) noexcept -> IntVect
+            {
+                return IntVect(AMREX_D_DECL(
+                                   static_cast<int>((p.pos(0)-plo[0])*dxi[0] - lo.x),
+                                   static_cast<int>((p.pos(1)-plo[1])*dxi[1] - lo.y),
+                                   static_cast<int>((p.pos(2)-plo[2])*dxi[2] - lo.z)));
+            });
+
+        return bins;
+    }
+
+}