CarpetIOHDF5: re-add ghost and symmetry points when removing buffers

NOTE: this assumes (like other parts of CarpetIOHDF5) that the number of
symmetry points is the number of ghost points.
NOTE: it likely outputs too many points when RotatingSymmetry is used and only
buffer points touch the symmetry boundary.

1 files changed, 26 insertions, 1 deletions

diff --git a/Carpet/CarpetIOHDF5/src/OutputSlice.cc b/Carpet/CarpetIOHDF5/src/OutputSlice.cc
index 1dfa369f5..4f1620c0e 100644
--- a/Carpet/CarpetIOHDF5/src/OutputSlice.cc
+++ b/Carpet/CarpetIOHDF5/src/OutputSlice.cc
@@ -1068,7 +1068,32 @@ namespace CarpetIOHDF5 {
     ibset exts(ibbox(lo,hi,str));
     // do grid arrays have buffer zones?
     if (not output_buffer_points) {
-      exts &= allactive;
+      ivect loghosts(0);
+      ivect highosts(0);
+      // possibly re-add ghost points and symmetry points 
+      // TODO: come up with a nicer way to do this
+      for (int d=0; d<groupdim; ++d) {
+        bool const output_lower_ghosts =
+          obnds[0][d]
+          ? (is_symbnd[2*d]
+             ? output_symmetry_points
+             : (output_boundary_points and out3D_outer_ghosts))
+          : (output_ghost_points and out3D_ghosts);
+        bool const output_upper_ghosts =
+          obnds[1][d]
+          ? (is_symbnd[2*d+1]
+             ? output_symmetry_points
+             : (output_boundary_points and out3D_outer_ghosts))
+          : (output_ghost_points and out3D_ghosts);
+
+        if (output_lower_ghosts) {
+          loghosts[d] = ghost_width[0][d];
+        }
+        if (output_upper_ghosts) {
+          highosts[d] = ghost_width[1][d];
+        }
+      }
+      exts = exts & allactive.expand(loghosts, highosts);
     }
 
     return exts;