Re-organise some refluxing internals

Re-organise some of the internal details of refluxing. Refluxing
requires restricting fluxes from fine to coarse grids. Previously,
Carpet would internally store bounding boxes that were offset by 1/2
grid point, and adjust (correct) these boxes in various places in a
rather ad-hoc manner. This is now cleaned up.

Remove the (unused) "accumulate" operator, and the (unused)
"accumulate" prolongation type.

2 files changed, 20 insertions, 6 deletions

diff --git a/Carpet/CarpetRegrid2/src/property.cc b/Carpet/CarpetRegrid2/src/property.cc
index 404a6e205..c1c8df42c 100644
--- a/Carpet/CarpetRegrid2/src/property.cc
+++ b/Carpet/CarpetRegrid2/src/property.cc
@@ -269,10 +269,10 @@ namespace CarpetRegrid2 {
       // again.
       
       // In detail, we first expand by reffact-1-N points, then expand
-      // to the next coarser grid, then expand back to the current
-      // grid, and expand by N points again. This sequence is correct
-      // for both vertex and cell centred grids, and N is determined
-      // by the number of buffer zones.
+      // (contract???) to the next coarser grid, then expand back to
+      // the current grid, and expand by N points again. This sequence
+      // is correct for both vertex and cell centred grids, and N is
+      // determined by the number of buffer zones.
       
       // N is the number of buffer zones modulo the refinement factor.
       // We cannot shrink the domain (since we cannot shrink
@@ -299,7 +299,17 @@ namespace CarpetRegrid2 {
     if (not snap_to_coarse) return true;
     
     ibset const snapped = snapped_regions (hh, dd, bnd, regions, rl);
-    return regions.AT(rl) == snapped;
+    
+    // We cannot test for equality, since the difference may be
+    // outside of the domain (and hence irrelevant)
+    // return regions.AT(rl) == snapped;
+    
+    // Test whether any part of the difference (i.e. that part of the
+    // level that would be added by snapping) is inside the domain. If
+    // the difference is outside, we can safely ignore it.
+    ibbox const& baseextent = hh.baseextent(0,rl);
+    ibset const difference = snapped - regions.AT(rl);
+    return (difference & baseextent).empty();
   }
   
   void snap_coarse::
diff --git a/Carpet/CarpetRegrid2/src/regrid.cc b/Carpet/CarpetRegrid2/src/regrid.cc
index d656224c2..752d600b7 100644
--- a/Carpet/CarpetRegrid2/src/regrid.cc
+++ b/Carpet/CarpetRegrid2/src/regrid.cc
@@ -367,6 +367,7 @@ namespace CarpetRegrid2 {
       // Enforce properties on this level
       for (int count=0;; ++count) {
         bool done_enforcing = true;
+        ibset const old_regions = regions.at(rl);
         for (vector<property*>::iterator
                pi = properties.begin(); pi != properties.end(); ++ pi)
         {
@@ -379,8 +380,11 @@ namespace CarpetRegrid2 {
         
         if (done_enforcing) break;
         
+        if (regions.at(rl) == old_regions) {
+          CCTK_WARN (CCTK_WARN_ABORT, "Could not enforce grid structure properties (not making any progress); giving up");
+        }
         if (count == 10) {
-          CCTK_WARN (CCTK_WARN_ABORT, "Could not enforce grid structure properties; giving up");
+          CCTK_WARN (CCTK_WARN_ABORT, "Could not enforce grid structure properties (maximum number of iterations reached); giving up");
         }
         if (count != 0) {
           // This may not be true. However, the previous version of