Bounded sampling

Src/Base/AMReX_MultiFabUtil.H

@@ -169,7 +169,8 @@ namespace amrex
     std::unique_ptr<MultiFab> get_slice_data(int dir, Real coord,
                                              const MultiFab& cc,
                                              const Geometry& geom, int start_comp, int ncomp,
-                                             bool interpolate=false);
+                                             bool interpolate=false,
+                                             RealBox bnd_rbx = RealBox());
 
     /**
      * \brief Get data in a cell of MultiFab/FabArray
@@ -188,7 +189,7 @@ namespace amrex
      * specified by a direction and a cell.
      */
     template <typename MF, std::enable_if_t<IsFabArray<MF>::value,int> FOO = 0>
-    MF get_line_data (MF const& mf, int dir, IntVect const& cell);
+    MF get_line_data (MF const& mf, int dir, IntVect const& cell, Box bnd_bx = Box());
 
     //! Return an iMultiFab that has the same BoxArray and DistributionMapping
     //! as the coarse MultiFab cmf. Cells covered by the coarsened fine grids
@@ -996,26 +997,45 @@ Vector<typename MF::value_type> get_cell_data (MF const& mf, IntVect const& cell
 }
 
 template <typename MF, std::enable_if_t<IsFabArray<MF>::value,int> FOO>
-MF get_line_data (MF const& mf, int dir, IntVect const& cell)
+MF get_line_data (MF const& mf, int dir, IntVect const& cell, Box bnd_bx)
 {
+    bool do_bnd = (!bnd_bx.isEmpty());
+
     BoxArray const& ba = mf.boxArray();
     DistributionMapping const& dm = mf.DistributionMap();
     const auto nboxes = static_cast<int>(ba.size());
 
     BoxList bl(ba.ixType());
     Vector<int> procmap;
     Vector<int> index_map;
-    for (int i = 0; i < nboxes; ++i) {
-        Box const& b = ba[i];
-        IntVect lo = cell;
-        lo[dir] = b.smallEnd(dir);
-        if (b.contains(lo)) {
-            IntVect hi = lo;
-            hi[dir] = b.bigEnd(dir);
-            Box b1d(lo,hi,b.ixType());
-            bl.push_back(b1d);
-            procmap.push_back(dm[i]);
-            index_map.push_back(i);
+    if (!do_bnd) {
+        for (int i = 0; i < nboxes; ++i) {
+            Box const& b = ba[i];
+            IntVect lo = cell;
+            lo[dir] = b.smallEnd(dir);
+            if (b.contains(lo)) {
+                IntVect hi = lo;
+                hi[dir] = b.bigEnd(dir);
+                Box b1d(lo,hi,b.ixType());
+                bl.push_back(b1d);
+                procmap.push_back(dm[i]);
+                index_map.push_back(i);
+            }
+        }
+    } else {
+        for (int i = 0; i < nboxes; ++i) {
+            Box const& b = ba[i];
+            IntVect lo = bnd_bx.smallEnd();
+            IntVect hi = bnd_bx.bigEnd();
+            if (!b.contains(lo)) { lo[dir] = b.smallEnd(dir); }
+            if (b.contains(lo)) {
+                if (!b.contains(hi)) { hi[dir] = b.bigEnd(dir); }
+                if (!bnd_bx.contains(hi)) { continue; }
+                Box b1d(lo,hi,b.ixType());
+                bl.push_back(b1d);
+                procmap.push_back(dm[i]);
+                index_map.push_back(i);
+            }
         }
     }
 

Src/Base/AMReX_MultiFabUtil.cpp

@@ -9,7 +9,7 @@ namespace {
     using namespace amrex;
 
     Box
-    getIndexBox(const RealBox& real_box, const Geometry& geom) {
+    getIndexBox (const RealBox& real_box, const Geometry& geom) {
         IntVect slice_lo, slice_hi;
 
         AMREX_D_TERM(slice_lo[0]=static_cast<int>(std::floor((real_box.lo(0) - geom.ProbLo(0))/geom.CellSize(0)));,
@@ -24,12 +24,11 @@ namespace {
     }
 
 
-    std::unique_ptr<MultiFab> allocateSlice(int dir, const MultiFab& cell_centered_data,
-                                            int ncomp, const Geometry& geom, Real dir_coord,
-                                            Vector<int>& slice_to_full_ba_map) {
+    std::unique_ptr<MultiFab> allocateSlice (int dir, const MultiFab& cell_centered_data,
+                                             int ncomp, const Geometry& geom, Real dir_coord,
+                                             Vector<int>& slice_to_full_ba_map, RealBox real_slice) {
 
         // Get our slice and convert to index space
-        RealBox real_slice = geom.ProbDomain();
         real_slice.setLo(dir, dir_coord);
         real_slice.setHi(dir, dir_coord);
         Box slice_box = getIndexBox(real_slice, geom);
@@ -550,7 +549,7 @@ namespace amrex
         return amrex::cast<FabArray<BaseFab<Long> > > (imf);
     }
 
-    std::unique_ptr<MultiFab> get_slice_data(int dir, Real coord, const MultiFab& cc, const Geometry& geom, int start_comp, int ncomp, bool interpolate) {
+    std::unique_ptr<MultiFab> get_slice_data(int dir, Real coord, const MultiFab& cc, const Geometry& geom, int start_comp, int ncomp, bool interpolate, RealBox bnd_rbx) {
 
         BL_PROFILE("amrex::get_slice_data");
 
@@ -559,9 +558,15 @@ namespace amrex
         }
 
         const auto geomdata = geom.data();
+        RealBox real_slice;
+        if (bnd_rbx.ok()) {
+            real_slice = bnd_rbx;
+        } else {
+            real_slice = geom.ProbDomain();
+        }
 
         Vector<int> slice_to_full_ba_map;
-        std::unique_ptr<MultiFab> slice = allocateSlice(dir, cc, ncomp, geom, coord, slice_to_full_ba_map);
+        std::unique_ptr<MultiFab> slice = allocateSlice(dir, cc, ncomp, geom, coord, slice_to_full_ba_map, real_slice);
 
         if (!slice) {
             return nullptr;