CarpetRegrid2: Ensure and test proper nesting

Introduce a parameter ensure_proper_nesting, which automatically
enlarges the levels of the grid hierarchy to ensure proper nesting.
(It does not enlarge the coarsest level.)  If this parameter is not
set, proper nesting is still tested, and a warning is emitted if it is
violated.

darcs-hash:20070425002649-dae7b-7431881af786bbae5048c5329807e90a4f7b2e29.gz

2 files changed, 66 insertions, 39 deletions

diff --git a/Carpet/CarpetRegrid2/param.ccl b/Carpet/CarpetRegrid2/param.ccl
index 97a42ce57..c97e71810 100644
--- a/Carpet/CarpetRegrid2/param.ccl
+++ b/Carpet/CarpetRegrid2/param.ccl
@@ -11,13 +11,17 @@ CCTK_INT min_distance "Minimum distance (in grid points) between coarse and fine
   0:* :: ""
 } 4
 
+BOOLEAN ensure_proper_nesting "Ensure proper nesting automatically" STEERABLE=always
+{
+} "yes"
+
+
+
 CCTK_REAL min_fraction "Minimum fraction of required refined points that need to be present in a refined region" STEERABLE=always
 {
   0:* :: ""
 } 0.9
 
-
-
 BOOLEAN snap_to_coarse "Ensure that the fine grid extent coincides with coarse grid points" STEERABLE=always
 {
 } "no"
diff --git a/Carpet/CarpetRegrid2/src/regrid.cc b/Carpet/CarpetRegrid2/src/regrid.cc
index ab178317e..5fe4ac8df 100644
--- a/Carpet/CarpetRegrid2/src/regrid.cc
+++ b/Carpet/CarpetRegrid2/src/regrid.cc
@@ -2,6 +2,7 @@
 #include <cmath>
 #include <cstdlib>
 #include <iostream>
+#include <sstream>
 #include <vector>
 
 #include <cctk.h>
@@ -265,7 +266,10 @@ namespace CarpetRegrid2 {
     vector <ibboxset> regions (1);
     
     // Set up coarsest level
-    regions.at(0) = ibboxset (hh.extent (0, 0, 0));
+    for (size_t c = 0; c < regss.at(0).size(); ++ c) {
+      regions.at(0) += regss.at(0).at(c).extent;
+    }
+    regions.at(0).normalize();
     
     // Loop over all centres
     for (int n = 0; n < num_centres; ++ n) {
@@ -301,8 +305,61 @@ namespace CarpetRegrid2 {
       
     regss.resize (regions.size());
     
+    assert (regions.size() > 0);
     for (size_t rl = regions.size() - 1; rl >= 1; -- rl) {
       
+      // Sanity check
+      assert (not regions.at(rl).empty());
+      
+      
+      
+      //
+      // Ensure that this grid contains all finer ones
+      //
+      if (rl < regions.size() - 1) {
+        bool is_properly_nested = true;
+        
+        ibbox const & coarse0 = * regions.at(rl).begin();
+        
+        i2vect const fdistance = dd.ghost_width;
+        i2vect const cdistance =
+          i2vect (min_distance + dd.prolongation_stencil_size());
+        
+        for (ibboxset::const_iterator ibb = regions.at(rl+1).begin();
+             ibb != regions.at(rl+1).end();
+             ++ ibb)
+        {
+          ibbox const & fbb = * ibb;
+          
+          bvect const lower_is_outer = fbb.lower() <= physical_ilower;
+          bvect const upper_is_outer = fbb.upper() >= physical_iupper;
+          b2vect const ob (lower_is_outer, upper_is_outer);
+          
+          ibbox const ebb = fbb.expand (i2vect (ob) * fdistance);
+          ibbox const cbb = ebb.expanded_for (coarse0);
+          ibbox const ecbb = cbb.expand (i2vect (ob) * cdistance);
+          
+          is_properly_nested = is_properly_nested and ecbb <= regions.at(rl);
+          
+          // Enlarge this level
+          if (ensure_proper_nesting) {
+            regions.at(rl) |= ecbb;
+          }
+        }
+        
+        if (ensure_proper_nesting) {
+          regions.at(rl).normalize();
+        } else {
+          if (not is_properly_nested) {
+            ostringstream msg;
+            msg << "Level " << rl << " of the refinement hierarchy is not properly nested.  It does not contain level " << (rl+1) << ".";
+            CCTK_WARN (CCTK_WARN_ALERT, msg.str().c_str());
+          }
+        }
+      }
+      
+      
+      
       //
       // Add buffer zones
       //
@@ -533,43 +590,9 @@ namespace CarpetRegrid2 {
       
       
       //
-      // Ensure that the coarser grids contain the finer ones
+      // Ensure that this grid is contained in the domain
       //
-      // TODO: move this to the top, and check the current grid
-      // instead of the next coarser one
-      if (rl > 0) {
-        assert (not regions.at(rl-1).empty());
-        ibbox const & coarse = * regions.at(rl-1).begin();
-        
-        i2vect const fdistance =
-          i2vect (0 * min_distance) + dd.ghost_width + dd.buffer_width;
-        i2vect const cdistance =
-          dd.buffer_width +
-          i2vect (min_distance + dd.prolongation_stencil_size());
-        
-        ibboxset coarsified;
-        for (ibboxset::const_iterator ibb = regions.at(rl).begin();
-             ibb != regions.at(rl).end();
-             ++ ibb)
-        {
-          ibbox const & fbb = * ibb;
-          ibbox const efbb = fbb.expand (fdistance[0], fdistance[1]);
-          ibbox const cbb = efbb.expanded_for(coarse);
-          ibbox const ecbb = cbb.expand (cdistance[0], cdistance[1]);
-          coarsified |= ecbb;
-        }
-        
-#if 0
-        // We cannot just enlarge the grid; all the nice properties
-        // that were ensured above are then broken
-        if (rl > 1) {
-          // Do not enlarge the coarsest grid
-          regions.at(rl-1) |= coarsified;
-          regions.at(rl-1).normalize();
-        }
-#endif
-        assert (coarsified <= regions.at(rl-1));
-      }
+      assert (regions.at(rl) <= hh.baseextents.at(0).at(rl));