Merge branch 'master' of carpetgit@carpetcode.dyndns.org:carpet

25 files changed, 670 insertions, 362 deletions

diff --git a/Carpet/Carpet/src/CarpetParamCheck.cc b/Carpet/Carpet/src/CarpetParamCheck.cc
index 30483d5db..fb9b40861 100644
--- a/Carpet/Carpet/src/CarpetParamCheck.cc
+++ b/Carpet/Carpet/src/CarpetParamCheck.cc
@@ -111,7 +111,7 @@ namespace Carpet {
         break;
       case all_timelevels:
         if (setup_method != init_all_levels) {
-          CCTK_PARAMWARN ("When you set neither Carpet::init_each_timelevel=yes nor Carpet::init_fill_timelevels=yes nor Carpet::init_3_timelevel=yes, then you must also use InitBase::initial_data_setup_method=\"init_all_levels\"");
+          CCTK_PARAMWARN ("When you set neither Carpet::init_each_timelevel=yes nor Carpet::init_fill_timelevels=yes nor Carpet::init_3_timelevels=yes, then you must also use InitBase::initial_data_setup_method=\"init_all_levels\"");
         }
         break;
       default:
diff --git a/Carpet/Carpet/src/Initialise.cc b/Carpet/Carpet/src/Initialise.cc
index 8d7cffff5..1b81dec11 100644
--- a/Carpet/Carpet/src/Initialise.cc
+++ b/Carpet/Carpet/src/Initialise.cc
@@ -271,6 +271,19 @@ namespace Carpet {
       CallRegridInitialMeta (cctkGH);
     }
     
+    // Poison early, since grid functions may be initialised in global
+    // loop-local mode, ane we must not overwrite them accidentally
+    for (int rl=0; rl<reflevels; ++rl) {
+      BEGIN_MGLEVEL_LOOP(cctkGH) {
+        ENTER_LEVEL_MODE (cctkGH, rl) {
+          
+          // Checking
+          Poison (cctkGH, alltimes, CCTK_GF);
+          
+        } LEAVE_LEVEL_MODE;
+      } END_MGLEVEL_LOOP;
+    } // for rl
+    
     for (int rl=0; rl<reflevels; ++rl) {
       BEGIN_MGLEVEL_LOOP(cctkGH) {
         ENTER_LEVEL_MODE (cctkGH, rl) {
@@ -284,9 +297,6 @@ namespace Carpet {
                     (do_global_mode ? " (global)" : ""),
                     (do_meta_mode ? " (meta)" : ""));
           
-          // Checking
-          Poison (cctkGH, alltimes, CCTK_GF);
-          
           // Timing statistics
           if (do_global_mode) {
             InitTimingStats (cctkGH);
@@ -327,6 +337,7 @@ namespace Carpet {
             
             // Set up the initial data
             ScheduleTraverse (where, "CCTK_INITIAL", cctkGH);
+            ScheduleTraverse (where, "CCTK_POSTINITIAL", cctkGH);
             
             if (init_fill_timelevels) {
               assert (tl==0);
@@ -347,7 +358,6 @@ namespace Carpet {
           
         } LEAVE_LEVEL_MODE;
       } END_MGLEVEL_LOOP;
-      
     } // for rl
     
     timer.stop();
@@ -396,7 +406,6 @@ namespace Carpet {
             ScheduleTraverse (where, "CCTK_POSTRESTRICTINITIAL", cctkGH);
           }
           
-          ScheduleTraverse (where, "CCTK_POSTINITIAL", cctkGH);
           ScheduleTraverse (where, "CCTK_POSTSTEP", cctkGH);
           
           PoisonCheck (cctkGH, alltimes);
@@ -667,7 +676,6 @@ namespace Carpet {
           
         } LEAVE_LEVEL_MODE;
       } END_MGLEVEL_LOOP;
-      
     } // for rl
     
     do_global_mode = old_do_global_mode;
diff --git a/Carpet/Carpet/src/Recompose.cc b/Carpet/Carpet/src/Recompose.cc
index 3027290de..9bd7d4066 100644
--- a/Carpet/Carpet/src/Recompose.cc
+++ b/Carpet/Carpet/src/Recompose.cc
@@ -91,13 +91,17 @@ namespace Carpet {
     }
     assert ((int)regsss.at(0).size() <= maxreflevels);
     for (int ml=0; ml<(int)regsss.size(); ++ml) {
+      int num_regions = 0;
       for (int rl=0; rl<(int)regsss.at(0).size(); ++rl) {
         // No empty levels
-        assert (regsss.at(ml).at(rl).size() > 0);
+        // (but allow some empty maps)
+        // assert (regsss.at(ml).at(rl).size() > 0);
+        num_regions += regsss.at(ml).at(rl).size();
         for (int c=0; c<(int)regsss.at(ml).at(rl).size(); ++c) {
           // Check sizes
-          // Do allow processors with zero grid points
-//           assert (all(regsss.at(rl).at(c).at(ml).extent.lower() <= regsss.at(rl).at(c).at(ml).extent.upper()));
+          // (but allow processors with zero grid points)
+          // assert (all(regsss.at(rl).at(c).at(ml).extent.lower() <=
+          //             regsss.at(rl).at(c).at(ml).extent.upper()));
           // Check strides
           const ivect str
             = (maxspacereflevelfact / spacereffacts.at(rl) * ipow(mgfact, ml));
@@ -107,6 +111,8 @@ namespace Carpet {
           assert (all(regsss.at(ml).at(rl).at(c).extent.upper() % str == 0));
         }
       }
+      // No empty levels
+      assert (num_regions > 0);
     }
   }
   
@@ -782,7 +788,7 @@ namespace Carpet {
           }
           for (int c=0; c<hh->components(rl); ++c) {
             ++ num_comps;
-            dh::dboxes const & b = dd->boxes.AT(m).AT(rl).AT(c);
+            dh::dboxes const & b = dd->boxes.AT(ml).AT(rl).AT(c);
             num_active_mem_points += num_gfs * b.active.size();
             num_owned_mem_points  += num_gfs * b.owned.size();
             num_total_mem_points  += num_gfs * b.exterior.size();
@@ -1137,7 +1143,6 @@ namespace Carpet {
     for (int n=0; n<nslices; ++n) {
       mynpoints.at(n) = int (floor (CCTK_REAL(1) * npoints_left * mynprocs.at(n)
                                     / nprocs_left + CCTK_REAL(0.5)));
-      assert (mynpoints.at(n) > 0);
       assert (mynprocs .at(n) > 0);
       npoints_left -= mynpoints.at(n);
       nprocs_left  -= mynprocs.at(n);
diff --git a/Carpet/Carpet/src/SetupGH.cc b/Carpet/Carpet/src/SetupGH.cc
index 109c3bcab..5e34fcd30 100644
--- a/Carpet/Carpet/src/SetupGH.cc
+++ b/Carpet/Carpet/src/SetupGH.cc
@@ -1221,6 +1221,8 @@ namespace Carpet {
       
     } else if (domain_from_coordbase) {
       
+      assert (not CCTK_IsFunctionAliased ("MultiPatch_GetDomainSpecification"));
+      
       // Ensure that CartGrid3D::type = "coordbase"
       ensure_CartGrid3D_type ();
       
@@ -1235,6 +1237,8 @@ namespace Carpet {
     } else {
       // Legacy code
       
+      assert (not CCTK_IsFunctionAliased ("MultiPatch_GetDomainSpecification"));
+      
       if (max_refinement_levels > 1) {
         // Ensure that CartGrid3D::avoid_origin = no
         ensure_CartGrid3D_avoid_origin ();
diff --git a/Carpet/CarpetIOASCII/src/ioascii.cc b/Carpet/CarpetIOASCII/src/ioascii.cc
index ade2299dd..f90770cde 100644
--- a/Carpet/CarpetIOASCII/src/ioascii.cc
+++ b/Carpet/CarpetIOASCII/src/ioascii.cc
@@ -358,7 +358,6 @@ namespace CarpetIOASCII {
   ::OutputVarAs (const cGH* const cctkGH,
 		 const char* const varname, const char* const alias)
   {
-    DECLARE_CCTK_ARGUMENTS;
     DECLARE_CCTK_PARAMETERS;
 
     int vindex = -1;
@@ -434,7 +433,7 @@ namespace CarpetIOASCII {
       for (int ml = 0; ml < mglevels; ++ml) {
         last_outputs[ml].resize (maxreflevels);
         for (int rl = 0; rl < maxreflevels; ++rl) {
-          last_outputs[ml][rl].resize (numelems, cctk_iteration - 1);
+          last_outputs[ml][rl].resize (numelems, cctkGH->cctk_iteration - 1);
         }
       }
       // TODO: this makes a copy of last_outputs, which is expensive
@@ -454,7 +453,7 @@ namespace CarpetIOASCII {
       elem = vindex;
     }
     int& last_output = thisfile->second.at(mglevel).at(reflevel).at(elem);
-    if (last_output == cctk_iteration) {
+    if (last_output == cctkGH->cctk_iteration) {
       // Has already been output during this iteration
       char * const fullname = CCTK_FullName (vindex);
       CCTK_VWarn (5, __LINE__, __FILE__, CCTK_THORNSTRING,
@@ -465,8 +464,8 @@ namespace CarpetIOASCII {
       free (fullname);
       return 0;
     }
-    assert (last_output < cctk_iteration);
-    last_output = cctk_iteration;
+    assert (last_output < cctkGH->cctk_iteration);
+    last_output = cctkGH->cctk_iteration;
 
     // Check for storage
     if (not CCTK_QueryGroupStorageI (cctkGH, group)) {
@@ -547,7 +546,9 @@ namespace CarpetIOASCII {
 	if (desired) {
 
 	  // Traverse all maps on this refinement and multigrid level
-	  BEGIN_MAP_LOOP(cctkGH, grouptype) {
+          int const m_min = 0;
+          int const m_max = grouptype == CCTK_GF ? Carpet::maps : 1;
+          for (int m = m_min; m < m_max; ++ m) {
 
             // Find the output offset.
             ivect offset(0);
@@ -555,17 +556,17 @@ namespace CarpetIOASCII {
               switch (outdim) {
               case 0:
                 offset[0] = GetGridOffset
-                  (cctkGH, 1,
+                  (cctkGH, m, 1,
                    "out%dD_point_xi", /*"out_point_xi"*/ NULL,
                    "out%dD_point_x",  /*"out_point_x"*/  NULL,
                    /*out_point_x*/ 0.0);
                 offset[1] = GetGridOffset
-                  (cctkGH, 2,
+                  (cctkGH, m, 2,
                    "out%dD_point_yi", /*"out_point_yi"*/ NULL,
                    "out%dD_point_y",  /*"out_point_y"*/  NULL,
                    /*out_point_y*/ 0.0);
                 offset[2] = GetGridOffset
-                  (cctkGH, 3,
+                  (cctkGH, m, 3,
                    "out%dD_point_zi", /*"out_point_zi"*/ NULL,
                    "out%dD_point_z",  /*"out_point_z"*/  NULL,
                    /*out_point_z*/ 0.0);
@@ -573,31 +574,31 @@ namespace CarpetIOASCII {
               case 1:
                 switch (dirs[0]) {
                 case 0:
-                  offset[1] = GetGridOffset (cctkGH, 2,
+                  offset[1] = GetGridOffset (cctkGH, m, 2,
                                              "out%dD_xline_yi", "out_xline_yi",
                                              "out%dD_xline_y",  "out_xline_y",
                                              out_xline_y);
-                  offset[2] = GetGridOffset (cctkGH, 3,
+                  offset[2] = GetGridOffset (cctkGH, m, 3,
                                              "out%dD_xline_zi", "out_xline_zi",
                                              "out%dD_xline_z",  "out_xline_z",
                                              out_xline_z);
                   break;
                 case 1:
-                  offset[0] = GetGridOffset (cctkGH, 1,
+                  offset[0] = GetGridOffset (cctkGH, m, 1,
                                              "out%dD_yline_xi", "out_yline_xi",
                                              "out%dD_yline_x",  "out_yline_x",
                                              out_yline_x);
-                  offset[2] = GetGridOffset (cctkGH, 3,
+                  offset[2] = GetGridOffset (cctkGH, m, 3,
                                              "out%dD_yline_zi", "out_yline_zi",
                                              "out%dD_yline_z",  "out_yline_z",
                                              out_yline_z);
                   break;
                 case 2:
-                  offset[0] = GetGridOffset (cctkGH, 1,
+                  offset[0] = GetGridOffset (cctkGH, m, 1,
                                              "out%dD_zline_xi", "out_zline_xi",
                                              "out%dD_zline_x",  "out_zline_x",
                                              out_zline_x);
-                  offset[1] = GetGridOffset (cctkGH, 2,
+                  offset[1] = GetGridOffset (cctkGH, m, 2,
                                              "out%dD_zline_yi", "out_zline_yi",
                                              "out%dD_zline_y",  "out_zline_y",
                                              out_zline_y);
@@ -615,19 +616,19 @@ namespace CarpetIOASCII {
               case 2:
                 if (dirs[0]==0 and dirs[1]==1) {
                   offset[2] = GetGridOffset
-                    (cctkGH, 3,
+                    (cctkGH, m, 3,
                      "out%dD_xyplane_zi", "out_xyplane_zi",
                      "out%dD_xyplane_z",  "out_xyplane_z",
                      out_xyplane_z);
                 } else if (dirs[0]==0 and dirs[1]==2) {
                   offset[1] = GetGridOffset
-                    (cctkGH, 2,
+                    (cctkGH, m, 2,
                      "out%dD_xzplane_yi", "out_xzplane_yi",
                      "out%dD_xzplane_y",  "out_xzplane_y",
                      out_xzplane_y);
                 } else if (dirs[0]==1 and dirs[1]==2) {
                   offset[0] = GetGridOffset
-                    (cctkGH, 1,
+                    (cctkGH, m, 1,
                      "out%dD_yzplane_xi", "out_yzplane_xi",
                      "out%dD_yzplane_x",  "out_yzplane_x",
                      out_yzplane_x);
@@ -653,7 +654,7 @@ namespace CarpetIOASCII {
               ostringstream filenamebuf;
               filenamebuf << my_out_dir << "/" << alias;
               if (maps > 1 and grouptype == CCTK_GF) {
-                filenamebuf << "." << Carpet::map;
+                filenamebuf << "." << m;
               }
               filenamebuf << ".";
               if (new_filename_scheme) {
@@ -766,7 +767,12 @@ namespace CarpetIOASCII {
 
             // Traverse and components on this multigrid and
             // refinement level and map
-            BEGIN_COMPONENT_LOOP(cctkGH, grouptype) {
+            int const c_min = 0;
+            int const c_max =
+              grouptype == CCTK_GF ?
+              vhh.at(m)->components(reflevel) :
+              CCTK_nProcs(cctkGH);
+            for (int c = c_min; c < c_max; ++ c) {
 
               cGroup groupdata;
               {
@@ -781,18 +787,14 @@ namespace CarpetIOASCII {
                 assert (not ierr);
               }
 
-              if (groupdata.disttype != CCTK_DISTRIB_CONSTANT
-                  or component == 0)
-              {
+              if (groupdata.disttype != CCTK_DISTRIB_CONSTANT or c == 0) {
 
-                const ggf* const ff
-                  = arrdata.at(group).at(Carpet::map).data.at(var);
+                const ggf* const ff = arrdata.at(group).at(m).data.at(var);
 
                 const int maxtl = output_all_timelevels ? num_tl : 1;
                 for (int tl=0; tl<maxtl; ++tl) {
 
-                  const gdata* const data
-                    = (*ff) (tl, my_reflevel, component, mglevel);
+                  const gdata* const data = (*ff) (tl, my_reflevel, c, mglevel);
                   ibbox ext = data->extent();
 
                   ivect lo = ext.lower();
@@ -855,13 +857,17 @@ namespace CarpetIOASCII {
                   rvect global_lower;
                   rvect coord_delta;
                   if (grouptype == CCTK_GF) {
+                    rvect const cctk_origin_space =
+                      origin_space.at(m).at(mglevel);
+                    rvect const cctk_delta_space =
+                      delta_space.at(m) * rvect (mglevelfact);
                     for (int d=0; d<dim; ++d) {
                       // lower boundary of Carpet's integer indexing
-                      global_lower[d] = cctkGH->cctk_origin_space[d];
+                      global_lower[d] = cctk_origin_space[d];
                       // grid spacing of Carpet's integer indexing
                       coord_delta[d] =
-                        cctkGH->cctk_delta_space[d] /
-                        vhh.at(Carpet::map)->baseextents.at(0).at(0).stride()[d];
+                        cctk_delta_space[d] /
+                        vhh.at(m)->baseextents.at(0).at(0).stride()[d];
                     }
                   } else {
                     for (int d=0; d<dim; ++d) {
@@ -877,7 +883,7 @@ namespace CarpetIOASCII {
                   ivect offset1;
                   if (grouptype == CCTK_GF) {
                     const ibbox& baseext
-                      = vhh.at(Carpet::map)->baseextents.at(mglevel).at(reflevel);
+                      = vhh.at(m)->baseextents.at(mglevel).at(reflevel);
                     offset1 = baseext.lower() + offset * ext.stride();
                   } else {
                     offset1 = offset * ext.stride();
@@ -893,9 +899,8 @@ namespace CarpetIOASCII {
                     int const numvars = CCTK_NumVarsInGroupI(group);
                     datas.resize (numvars);
                     for (int n=0; n<numvars; ++n) {
-                      const ggf* const ff1
-                        = arrdata.at(group).at(Carpet::map).data.at(n);
-                      datas.at(n) = (*ff1) (tl, my_reflevel, component, mglevel);
+                      const ggf* const ff1 = arrdata.at(group).at(m).data.at(n);
+                      datas.at(n) = (*ff1) (tl, my_reflevel, c, mglevel);
                     }
                   } else {
                     datas.resize (1);
@@ -903,7 +908,7 @@ namespace CarpetIOASCII {
                   }
                   WriteASCII (file, datas, ext, vindex, cctkGH->cctk_iteration,
                               offset1, dirs,
-                              my_reflevel, mglevel, Carpet::map, component, tl,
+                              my_reflevel, mglevel, m, c, tl,
                               coord_time, coord_lower, coord_upper);
 
                   // Append EOL after every component
@@ -917,9 +922,9 @@ namespace CarpetIOASCII {
 
                 } // for tl
 
-              } // if distrib!=CONSTANT or component==0
+              } // if distrib!=CONSTANT or c==0
 
-            } END_COMPONENT_LOOP;
+            } // for c
 
             // Append EOL after every complete set of components
             if (dist::rank()==ioproc) {
@@ -937,7 +942,7 @@ namespace CarpetIOASCII {
             CCTK_REAL const io_bytes = io_bytes_end - io_bytes_begin;
             EndTimingIO (cctkGH, io_files, io_bytes, false);
 
-          } END_MAP_LOOP;
+          } // for m
 
 	} // if (desired)
 
@@ -993,7 +998,7 @@ namespace CarpetIOASCII {
 
   template<int outdim>
   int IOASCII<outdim>
-  ::GetGridOffset (const cGH* const cctkGH, const int dir,
+  ::GetGridOffset (const cGH* const cctkGH, const int m, const int dir,
 		   const char* const itempl, const char* const iglobal,
 		   const char* const ctempl, const char* const cglobal,
 		   const CCTK_REAL cfallback)
@@ -1011,7 +1016,7 @@ namespace CarpetIOASCII {
       assert (pcoord);
       const CCTK_REAL coord = *pcoord;
       assert (ptype == PARAMETER_REAL);
-      return CoordToOffset (cctkGH, dir, coord, 0);
+      return CoordToOffset (cctkGH, m, dir, coord, 0);
     }
 
     // Second choice: explicit index
@@ -1042,7 +1047,7 @@ namespace CarpetIOASCII {
         assert (pcoord);
         const CCTK_REAL coord = *pcoord;
         assert (ptype == PARAMETER_REAL);
-        return CoordToOffset (cctkGH, dir, coord, 0);
+        return CoordToOffset (cctkGH, m, dir, coord, 0);
       }
     }
 
@@ -1062,22 +1067,30 @@ namespace CarpetIOASCII {
     }
 
     // Fifth choice: default coordinate
-    return CoordToOffset (cctkGH, dir, cfallback, 0);
+    return CoordToOffset (cctkGH, m, dir, cfallback, 0);
   }
 
 
 
   template<int outdim>
   int IOASCII<outdim>
-  ::CoordToOffset (const cGH* cctkGH, const int dir, const CCTK_REAL coord,
-		   const int ifallback)
+  ::CoordToOffset (const cGH* cctkGH, const int m, const int dir,
+                   const CCTK_REAL coord, const int ifallback)
   {
-    assert (dir>=1 and dir<=dim);
+    assert (m>=0 and m<Carpet::maps and dir>=1 and dir<=dim);
+
+    assert (mglevel!=-1 and reflevel!=-1 and Carpet::map==-1);
 
-    assert (mglevel!=-1 and reflevel!=-1 and Carpet::map!=-1);
+    rvect const cctk_origin_space = origin_space.at(m).at(mglevel);
+    rvect const cctk_delta_space  = delta_space.at(m) * rvect (mglevelfact);
+    ivect const cctk_levfac = spacereffacts.at (reflevel);
+    ibbox const & coarseext = vhh.at(m)->baseextents.at(mglevel).at(0       );
+    ibbox const & baseext   = vhh.at(m)->baseextents.at(mglevel).at(reflevel);
+    ivect const cctk_levoff      = baseext.lower() - coarseext.lower();
+    ivect const cctk_levoffdenom = baseext.stride();
 
-    const CCTK_REAL delta = cctkGH->cctk_delta_space[dir-1] / cctkGH->cctk_levfac[dir-1];
-    const CCTK_REAL lower = cctkGH->cctk_origin_space[dir-1] + cctkGH->cctk_delta_space[dir-1] / cctkGH->cctk_levfac[dir-1] * cctkGH->cctk_levoff[dir-1] / cctkGH->cctk_levoffdenom[dir-1];
+    const CCTK_REAL delta = cctk_delta_space[dir-1] / cctk_levfac[dir-1];
+    const CCTK_REAL lower = cctk_origin_space[dir-1] + cctk_delta_space[dir-1] / cctk_levfac[dir-1] * cctk_levoff[dir-1] / cctk_levoffdenom[dir-1];
 
     const CCTK_REAL rindex = (coord - lower) / delta;
     int cindex = (int)floor(rindex + 0.75);
@@ -1317,7 +1330,7 @@ namespace CarpetIOASCII {
 	  const vect<int,3> str = gfext.stride();
 	  const bbox<int,3> ext(lo,up,str);
           
-	  gh const & hh = *vhh.at(Carpet::map);
+	  gh const & hh = *vhh.at(m);
 	  ibbox const & base = hh.baseextents.at(mglevel).at(reflevel);
           
        	  assert (base.stride()[0] ==  base.stride()[1]
diff --git a/Carpet/CarpetIOASCII/src/ioascii.hh b/Carpet/CarpetIOASCII/src/ioascii.hh
index 921490af9..6f0bb74dc 100644
--- a/Carpet/CarpetIOASCII/src/ioascii.hh
+++ b/Carpet/CarpetIOASCII/src/ioascii.hh
@@ -46,12 +46,12 @@ namespace CarpetIOASCII {
     static int TimeToOutput (const cGH* cctkGH, int vindex);
     static int TriggerOutput (const cGH* cctkGH, int vindex);
 
-    static int GetGridOffset (const cGH* cctkGH, int dir,
+    static int GetGridOffset (const cGH* cctkGH, int m, int dir,
                               const char* itempl, const char* iglobal,
                               const char* ctempl, const char* cglobal,
                               CCTK_REAL cfallback);
-    static int CoordToOffset (const cGH* cctkGH, int dir, CCTK_REAL coord,
-                              int ifallback);
+    static int CoordToOffset (const cGH* cctkGH, int m, int dir,
+                              CCTK_REAL coord, int ifallback);
 
     static void CheckSteerableParameters (const cGH* cctkGH);
     static const char* GetStringParameter (const char* parametertemplate);
diff --git a/Carpet/CarpetInterp/src/interp.cc b/Carpet/CarpetInterp/src/interp.cc
index 3e6d2a801..2ff115029 100644
--- a/Carpet/CarpetInterp/src/interp.cc
+++ b/Carpet/CarpetInterp/src/interp.cc
@@ -3,6 +3,7 @@
 #include <cmath>
 #include <cstdlib>
 #include <cstring>
+#include <iostream>
 #include <vector>
 
 #include <mpi.h>
@@ -387,6 +388,7 @@ namespace CarpetInterp {
                 coords_list, coords_buffer,
                 dstprocs, N_points_to,
                 rlev, home, allhomecnts);
+    //cout << "CarpetInterp: N_points_to=" << N_points_to << endl;
 
     //////////////////////////////////////////////////////////////////////
     // Communicate the number of points each processor is going to communicate
@@ -402,6 +404,7 @@ namespace CarpetInterp {
     MPI_Alltoall (&N_points_to[0],   1, dist::datatype (N_points_to[0]),
                   &N_points_from[0], 1, dist::datatype (N_points_from[0]),
                   dist::comm());
+    //cout << "CarpetInterp: N_points_from=" << N_points_from << endl;
 
     //////////////////////////////////////////////////////////////////////
     // Communicate the interpolation coordinates
@@ -507,16 +510,38 @@ namespace CarpetInterp {
         recvdispl.AT(p) = recvdispl.AT(p-1) + recvcnt.AT(p-1);
       }
 
-      // Re-order source maps into a temporary buffer
-      vector<CCTK_INT> tmp (N_interp_points);
-      vector<int> tmpcnts (N_points_to.size());
-      for (int n = 0; n < N_interp_points; n++) {
-        int const p = dstprocs.AT(n);
-        int const offset = senddispl.AT(p) + tmpcnts.AT(p);
-        tmp.AT(offset) = source_map.AT(n);
-        tmpcnts.AT(p)++;
+      // Set up the per-component maps
+      // as offset into the single communication send buffer
+      vector<CCTK_INT> tmp (source_map.size());
+      vector<CCTK_INT*> maps (allhomecnts.size());
+      {
+        int offset = 0;
+        for (int c = 0; c < (int)allhomecnts.size(); c++) {
+          maps.AT(c) = &tmp.front() + offset;
+          offset += allhomecnts.AT(c);
+        }
+        assert (offset == N_interp_points);
+      }
+  
+      // Copy the input maps into the communication send buffer
+      {
+        // totalhomecnts is the accumulated number of points over all components
+        vector<int> totalhomecnts (allhomecnts.size());
+        for (int idx = 1; idx < (int)allhomecnts.size(); idx++) {
+          totalhomecnts.AT(idx) =
+            totalhomecnts.AT(idx-1) + allhomecnts.AT(idx-1);
+        }
+  
+        vector<int> tmpcnts (allhomecnts.size(), 0);
+        for (int n = 0; n < N_interp_points; n++) {
+          int const idx = component_idx
+            (dstprocs.AT(n), source_map.AT(n), rlev.AT(n), home.AT(n));
+          assert (tmpcnts.AT(idx) < allhomecnts.AT(idx));
+          maps.AT(idx)[tmpcnts.AT(idx)] = source_map.AT(n);
+          tmpcnts.AT(idx)++;
+        }
+        assert (tmpcnts == allhomecnts);
       }
-      assert (tmpcnts == N_points_to);
 
       // Transmit the source maps
       source_map.resize (N_points_local);
@@ -556,10 +581,7 @@ namespace CarpetInterp {
     //////////////////////////////////////////////////////////////////////
     // Map (local) interpolation points to components
     //////////////////////////////////////////////////////////////////////
-    rlev.resize (N_points_local);          // reflevel of (local) point n
-    home.resize (N_points_local);          // component of (local) point n
     vector<int> srcprocs (N_points_local); // which processor requested point n
-    vector<int> homecnts (allhomecnts.size());  // points per components
 
     // Remember from where point n came
     {
@@ -588,6 +610,9 @@ namespace CarpetInterp {
     // (One could argue though that the per-point status codes as returned
     // by the local interpolator could be used to determine the global
     // interpolator return value instead.)
+    rlev.resize (N_points_local);          // reflevel of (local) point n
+    home.resize (N_points_local);          // component of (local) point n
+    vector<int> homecnts (allhomecnts.size(), 0);  // points per components
     map_points (cctkGH, coord_system_handle, coord_group,
                 ml, minrl, maxrl, maxncomps, N_dims, N_points_local,
                 source_map,
@@ -912,6 +937,7 @@ namespace CarpetInterp {
         GetCoordRange (cctkGH, m, mglevel, dim,
                        & gsh[0],
                        & lower.AT(m)[0], & upper.AT(m)[0], & delta.AT(m)[0]);
+        //cout << "CarpetInterp distribute: m=" << m << " gsh=" << gsh << " lower=" << lower.AT(m) << " upper=" << upper.AT(m) << " delta=" << delta.AT(m) << endl;
         delta.AT(m) /= maxspacereflevelfact;
       }
       break;
@@ -931,6 +957,7 @@ namespace CarpetInterp {
         assert (iret == 0);
       }
       
+#warning "Why can the map number for grid arrays be larger than 0?"
       for (int m = 0; m < maps; ++m) {
         ibbox const & baseextent =
           arrdata.AT(coord_group).AT(m).hh->baseextents.AT(mglevel).AT(0);
@@ -946,6 +973,15 @@ namespace CarpetInterp {
       assert (0);
     }
 
+    for (int m=0; m<maps; ++m) {
+      gh const * const hh = arrdata.AT(coord_group).AT(m).hh;
+      for (int rl=minrl; rl<maxrl; ++rl) {
+        for (int c = 0; c < hh->components(rl); ++c) {
+          //cout << "CarpetInterp: gh: m=" << m << " rl=" << rl << " c=" << c << " ext=" << hh->extent(0,rl,c) << " ob=" << hh->outer_boundaries(rl,c) << " p=" << hh->processor(rl, c) << " local=" << hh->is_local(rl,c) << endl;
+        }
+      }
+    }
+
     // Assign interpolation points to processors/components
     for (int n = 0; n < npoints; ++n) {
 
@@ -965,6 +1001,7 @@ namespace CarpetInterp {
 
       // Find the component that this grid point belongs to
       int rl = -1, c = -1;
+      //cout << "CarpetInterp: assign: n=" << n << " m=" << m << " pos=" << pos << endl;
       if (all (pos >= lower.AT(m) and pos <= upper.AT(m))) {
         // Try finer levels first
         for (rl = maxrl-1; rl >= minrl; --rl) {
@@ -977,6 +1014,7 @@ namespace CarpetInterp {
 
           // TODO: use something faster than a linear search
           for (c = 0; c < hh->components(rl); ++c) {
+            //cout << "   rl=" << rl << " ipos=" << ipos << " c=" << c << " ext=" << hh->extent(ml,rl,c) << endl;
             if (hh->extent(ml,rl,c).contains(ipos)) {
               goto found;
             }
@@ -991,9 +1029,8 @@ namespace CarpetInterp {
       // Only warn once
       if (map_onto_processors) {
         CCTK_VWarn (CCTK_WARN_PICKY, __LINE__, __FILE__, CCTK_THORNSTRING,
-                    "Interpolation point #%d at [%g,%g,%g] is not on "
-                    "any grid patch",
-                    n, (double)pos[0], (double)pos[1], (double)pos[2]);
+                    "Interpolation point #%d at [%g,%g,%g] of patch #%d is not on any component",
+                    n, (double)pos[0], (double)pos[1], (double)pos[2], (int)m);
       }
     found:
       assert (rl >= minrl and rl < maxrl);
@@ -1006,6 +1043,7 @@ namespace CarpetInterp {
       rlev.AT(n) = rl;
       home.AT(n) = c;
       ++ homecnts.AT(component_idx (procs.AT(n), m, rl, c));
+      //cout << "   p=" << procs.AT(n) << endl;
 
     } // for n
   }
@@ -1097,6 +1135,7 @@ namespace CarpetInterp {
               for (int p = 0; p < dist::size(); p++) {
                 int const idx =
                   component_idx (p, Carpet::map,reflevel,component);
+                //cout << "CarpetInterp: interpolate_single_component: rl=" << reflevel << " map=" << (Carpet::map) << " component=" << component << " p=" << p << endl;
                 if (homecnts.AT(idx) > 0) {
                   interpolate_single_component
                     (cctkGH, coord_system_handle, coord_group,
@@ -1313,6 +1352,7 @@ namespace CarpetInterp {
       jvect gsh;
       GetCoordRange (cctkGH, Carpet::map, mglevel, dim,
                      & gsh[0], & lower[0], & upper[0], & delta[0]);
+      //cout << "CarpetInterp single: m=" << (Carpet::map) << " gsh=" << gsh << " lower=" << lower << " upper=" << upper << " delta=" << delta << endl;
       delta /= maxspacereflevelfact;
       break;
     }
diff --git a/Carpet/CarpetLib/src/bbox.cc b/Carpet/CarpetLib/src/bbox.cc
index 97a8641ef..214f26a66 100644
--- a/Carpet/CarpetLib/src/bbox.cc
+++ b/Carpet/CarpetLib/src/bbox.cc
@@ -2,6 +2,7 @@
 #include <cassert>
 #include <iostream>
 #include <limits>
+#include <typeinfo>
 
 #include "cctk.h"
 
@@ -27,8 +28,9 @@ void bbox<T,D>::assert_bbox_limits () const
         any (_upper >= numeric_limits<T>::max() / 2) or
         any (_upper <= numeric_limits<T>::min() / 2))
     {
-      CCTK_WARN (CCTK_WARN_ABORT,
-                 "Tried to create a very large bbox -- it is likely that this would lead to an integer overflow");
+      CCTK_VWarn (CCTK_WARN_ABORT, __LINE__, __FILE__, CCTK_THORNSTRING,
+                  "Tried to create a very large bbox of type %s -- it is likely that this would lead to an integer overflow",
+                  typeid(*this).name());
     }
   }
 }
@@ -46,8 +48,9 @@ typename bbox<T,D>::size_type bbox<T,D>::size () const {
   size_type sz = 1, max = numeric_limits<size_type>::max();
   for (int d=0; d<D; ++d) {
     if (sh[d] > max) {
-      CCTK_WARN (CCTK_WARN_ABORT,
-                 "bbox size is too large -- integer overflow");
+      CCTK_VWarn (CCTK_WARN_ABORT, __LINE__, __FILE__, CCTK_THORNSTRING,
+                  "size of bbox of type %s is too large -- integer overflow",
+                  typeid(*this).name());
     }
     sz *= sh[d];
     max /= sh[d];
diff --git a/Carpet/CarpetLib/src/copy_3d.cc b/Carpet/CarpetLib/src/copy_3d.cc
index 5c9950f50..36a48df40 100644
--- a/Carpet/CarpetLib/src/copy_3d.cc
+++ b/Carpet/CarpetLib/src/copy_3d.cc
@@ -103,9 +103,9 @@ namespace CarpetLib {
     
     // Loop over region
 #pragma omp parallel for
-    for (ptrdiff_t k=0; k<regkext; ++k) {
-      for (ptrdiff_t j=0; j<regjext; ++j) {
-        for (ptrdiff_t i=0; i<regiext; ++i) {
+    for (int k=0; k<regkext; ++k) {
+      for (int j=0; j<regjext; ++j) {
+        for (int i=0; i<regiext; ++i) {
           
           dst [DSTIND3(i, j, k)] = src [SRCIND3(i, j, k)];
           
diff --git a/Carpet/CarpetLib/src/data.cc b/Carpet/CarpetLib/src/data.cc
index bad80a00a..da1f8ab56 100644
--- a/Carpet/CarpetLib/src/data.cc
+++ b/Carpet/CarpetLib/src/data.cc
@@ -58,7 +58,8 @@ call_operator (void
 #  if ! defined (CARPET_OPTIMISE)
   ibset allregbboxes;
 #  endif
-  _Pragma ("omp parallel") {
+#pragma omp parallel
+  {
     int const num_threads = omp_get_num_threads();
     int const thread_num = omp_get_thread_num();
     // Parallelise in z direction
@@ -83,9 +84,8 @@ call_operator (void
     if (not myregbbox.empty()) {
       (* the_operator) (src, srcext, dst, dstext, srcbbox, dstbbox, myregbbox);
 #  if ! defined (NDEBUG) && ! defined (CARPET_OPTIMISE)
-      _Pragma ("omp critical") {
-        allregbboxes += myregbbox;
-      }
+#pragma omp critical
+      allregbboxes += myregbbox;
 #  endif
     }
   }
diff --git a/Carpet/CarpetLib/src/gh.cc b/Carpet/CarpetLib/src/gh.cc
index 6a7738048..fa9ff2125 100644
--- a/Carpet/CarpetLib/src/gh.cc
+++ b/Carpet/CarpetLib/src/gh.cc
@@ -108,7 +108,7 @@ regrid (mregs const & regs)
   // Check component consistency
   for (int ml=0; ml<mglevels(); ++ml) {
     for (int rl=0; rl<reflevels(); ++rl) {
-      assert (components(rl)>0);
+      assert (components(rl)>=0);
       for (int c=0; c<components(rl); ++c) {
         ibbox const & b  = extent(ml,rl,c);
         ibbox const & b0 = extent(ml,rl,0);
@@ -136,31 +136,34 @@ regrid (mregs const & regs)
     bool have_error = false;
     for (int ml=0; ml<mglevels(); ++ml) {
       for (int rl=1; rl<reflevels(); ++rl) {
-        assert (all (extent(ml,rl,0).stride() * reffacts.AT(rl) ==
-                     extent(ml,rl-1,0).stride() * reffacts.AT(rl-1)));
-        // Check contained-ness:
-        // first take all coarse grids
-        ibset coarse;
-        for (int c=0; c<components(rl-1); ++c) {
-          coarse += extent(ml,rl-1,c);
-        }
-        coarse.normalize();
-        // then check all fine grids
-        for (int c=0; c<components(rl); ++c) {
-          ibbox const & fine =
-            extent(ml,rl,c).contracted_for(extent(ml,rl-1,0));
-          if (not (fine <= coarse)) {
-            if (not have_error) {
-              cout << "The following components are not properly nested, i.e.," << endl
-                   << "they are not contained within the next coarser level's components:" << endl;
-              have_error = true;
-            }
-            cout << "   ml " << ml << " rl " << rl << " c " << c << ":   "
-                 << fine << endl;
+        if (components(rl)>0) {
+          assert (components(rl-1)>0);
+          assert (all (extent(ml,rl,0).stride() * reffacts.AT(rl) ==
+                       extent(ml,rl-1,0).stride() * reffacts.AT(rl-1)));
+          // Check contained-ness:
+          // first take all coarse grids
+          ibset coarse;
+          for (int c=0; c<components(rl-1); ++c) {
+            coarse += extent(ml,rl-1,c);
           }
-        } // for c
-      }   // for rl
-    }     // for ml
+          coarse.normalize();
+          // then check all fine grids
+          for (int c=0; c<components(rl); ++c) {
+            ibbox const & fine =
+              extent(ml,rl,c).contracted_for(extent(ml,rl-1,0));
+            if (not (fine <= coarse)) {
+              if (not have_error) {
+                cout << "The following components are not properly nested, i.e.," << endl
+                     << "they are not contained within the next coarser level's components:" << endl;
+                have_error = true;
+              }
+              cout << "   ml " << ml << " rl " << rl << " c " << c << ":   "
+                   << fine << endl;
+            }
+          } // for c
+        }   // if c
+      }     // for rl
+    }       // for ml
     if (have_error) {
       cout << "The grid hierarchy is:" << endl;
       for (int ml=0; ml<mglevels(); ++ml) {
diff --git a/Carpet/CarpetLib/src/interpolate_3d_2tl.cc b/Carpet/CarpetLib/src/interpolate_3d_2tl.cc
index 78170b1ee..0b984b142 100644
--- a/Carpet/CarpetLib/src/interpolate_3d_2tl.cc
+++ b/Carpet/CarpetLib/src/interpolate_3d_2tl.cc
@@ -121,9 +121,9 @@ namespace CarpetLib {
     
     // Loop over region
 #pragma omp parallel for
-    for (ptrdiff_t k=0; k<regkext; ++k) {
-      for (ptrdiff_t j=0; j<regjext; ++j) {
-        for (ptrdiff_t i=0; i<regiext; ++i) {
+    for (int k=0; k<regkext; ++k) {
+      for (int j=0; j<regjext; ++j) {
+        for (int i=0; i<regiext; ++i) {
           
           dst [DSTIND3(i, j, k)] =
             + s1fac * src1 [SRCIND3(i, j, k)]
diff --git a/Carpet/CarpetLib/src/interpolate_3d_3tl.cc b/Carpet/CarpetLib/src/interpolate_3d_3tl.cc
index 9a1d0e5d7..c0e8b44fd 100644
--- a/Carpet/CarpetLib/src/interpolate_3d_3tl.cc
+++ b/Carpet/CarpetLib/src/interpolate_3d_3tl.cc
@@ -125,9 +125,9 @@ namespace CarpetLib {
     
     // Loop over region
 #pragma omp parallel for
-    for (ptrdiff_t k=0; k<regkext; ++k) {
-      for (ptrdiff_t j=0; j<regjext; ++j) {
-        for (ptrdiff_t i=0; i<regiext; ++i) {
+    for (int k=0; k<regkext; ++k) {
+      for (int j=0; j<regjext; ++j) {
+        for (int i=0; i<regiext; ++i) {
           
           dst [DSTIND3(i, j, k)] =
             + s1fac * src1 [SRCIND3(i, j, k)]
diff --git a/Carpet/CarpetLib/src/interpolate_3d_4tl.cc b/Carpet/CarpetLib/src/interpolate_3d_4tl.cc
index 9892d0bbf..7d4c7fe26 100644
--- a/Carpet/CarpetLib/src/interpolate_3d_4tl.cc
+++ b/Carpet/CarpetLib/src/interpolate_3d_4tl.cc
@@ -132,9 +132,9 @@ namespace CarpetLib {
     
     // Loop over region
 #pragma omp parallel for
-    for (ptrdiff_t k=0; k<regkext; ++k) {
-      for (ptrdiff_t j=0; j<regjext; ++j) {
-        for (ptrdiff_t i=0; i<regiext; ++i) {
+    for (int k=0; k<regkext; ++k) {
+      for (int j=0; j<regjext; ++j) {
+        for (int i=0; i<regiext; ++i) {
           
           dst [DSTIND3(i, j, k)] =
             + s1fac * src1 [SRCIND3(i, j, k)]
diff --git a/Carpet/CarpetLib/src/interpolate_3d_5tl.cc b/Carpet/CarpetLib/src/interpolate_3d_5tl.cc
index abad807b1..f4204ea68 100644
--- a/Carpet/CarpetLib/src/interpolate_3d_5tl.cc
+++ b/Carpet/CarpetLib/src/interpolate_3d_5tl.cc
@@ -137,9 +137,9 @@ namespace CarpetLib {
     
     // Loop over region
 #pragma omp parallel for
-    for (ptrdiff_t k=0; k<regkext; ++k) {
-      for (ptrdiff_t j=0; j<regjext; ++j) {
-        for (ptrdiff_t i=0; i<regiext; ++i) {
+    for (int k=0; k<regkext; ++k) {
+      for (int j=0; j<regjext; ++j) {
+        for (int i=0; i<regiext; ++i) {
           
           dst [DSTIND3(i, j, k)] =
             + s1fac * src1 [SRCIND3(i, j, k)]
diff --git a/Carpet/CarpetLib/src/interpolate_eno_3d_3tl.cc b/Carpet/CarpetLib/src/interpolate_eno_3d_3tl.cc
index ef3a69053..e84059d10 100644
--- a/Carpet/CarpetLib/src/interpolate_eno_3d_3tl.cc
+++ b/Carpet/CarpetLib/src/interpolate_eno_3d_3tl.cc
@@ -155,9 +155,9 @@ namespace CarpetLib {
     
     // Loop over region
 #pragma omp parallel for
-    for (ptrdiff_t k=0; k<regkext; ++k) {
-      for (ptrdiff_t j=0; j<regjext; ++j) {
-        for (ptrdiff_t i=0; i<regiext; ++i) {
+    for (int k=0; k<regkext; ++k) {
+      for (int j=0; j<regjext; ++j) {
+        for (int i=0; i<regiext; ++i) {
           
           T const s1 = src1 [SRCIND3(i, j, k)];
           T const s2 = src2 [SRCIND3(i, j, k)];
diff --git a/Carpet/CarpetLib/src/restrict_3d_rf2.cc b/Carpet/CarpetLib/src/restrict_3d_rf2.cc
index e64c4f6e5..26031f304 100644
--- a/Carpet/CarpetLib/src/restrict_3d_rf2.cc
+++ b/Carpet/CarpetLib/src/restrict_3d_rf2.cc
@@ -102,9 +102,9 @@ namespace CarpetLib {
     
     // Loop over coarse region
 #pragma omp parallel for
-    for (ptrdiff_t k=0; k<regkext; ++k) {
-      for (ptrdiff_t j=0; j<regjext; ++j) {
-        for (ptrdiff_t i=0; i<regiext; ++i) {
+    for (int k=0; k<regkext; ++k) {
+      for (int j=0; j<regjext; ++j) {
+        for (int i=0; i<regiext; ++i) {
           
           dst [DSTIND3(i, j, k)] = src [SRCIND3(2*i, 2*j, 2*k)];
           
diff --git a/Carpet/CarpetRegrid2/interface.ccl b/Carpet/CarpetRegrid2/interface.ccl
index 7b55ed523..aa6f799da 100644
--- a/Carpet/CarpetRegrid2/interface.ccl
+++ b/Carpet/CarpetRegrid2/interface.ccl
@@ -47,6 +47,42 @@ CCTK_INT FUNCTION ConvertFromPhysicalBoundary \
      CCTK_REAL IN  ARRAY spacing)
 USES FUNCTION ConvertFromPhysicalBoundary
 
+# The location of the boundary points
+CCTK_INT FUNCTION MultiPatch_GetBoundarySpecification \
+    (CCTK_INT IN        map,                          \
+     CCTK_INT IN        size,                         \
+     CCTK_INT OUT ARRAY nboundaryzones,               \
+     CCTK_INT OUT ARRAY is_internal,                  \
+     CCTK_INT OUT ARRAY is_staggered,                 \
+     CCTK_INT OUT ARRAY shiftout)
+USES FUNCTION MultiPatch_GetBoundarySpecification
+
+# The overall size of the domain
+CCTK_INT FUNCTION MultiPatch_GetDomainSpecification \
+    (CCTK_INT  IN        map,                       \
+     CCTK_INT  IN        size,                      \
+     CCTK_REAL OUT ARRAY physical_min,              \
+     CCTK_REAL OUT ARRAY physical_max,              \
+     CCTK_REAL OUT ARRAY interior_min,              \
+     CCTK_REAL OUT ARRAY interior_max,              \
+     CCTK_REAL OUT ARRAY exterior_min,              \
+     CCTK_REAL OUT ARRAY exterior_max,              \
+     CCTK_REAL OUT ARRAY spacing)
+USES FUNCTION MultiPatch_GetDomainSpecification
+
+# Conversion between boundary types
+CCTK_INT FUNCTION MultiPatch_ConvertFromPhysicalBoundary \
+    (CCTK_INT  IN        map,                            \
+     CCTK_INT  IN        size,                           \
+     CCTK_REAL IN  ARRAY physical_min,                   \
+     CCTK_REAL IN  ARRAY physical_max,                   \
+     CCTK_REAL OUT ARRAY interior_min,                   \
+     CCTK_REAL OUT ARRAY interior_max,                   \
+     CCTK_REAL OUT ARRAY exterior_min,                   \
+     CCTK_REAL OUT ARRAY exterior_max,                   \
+     CCTK_REAL IN  ARRAY spacing)
+USES FUNCTION MultiPatch_ConvertFromPhysicalBoundary
+
 
 
 # The true prototype of the routine below:
diff --git a/Carpet/CarpetRegrid2/src/regrid.cc b/Carpet/CarpetRegrid2/src/regrid.cc
index 983f7acba..c08705d29 100644
--- a/Carpet/CarpetRegrid2/src/regrid.cc
+++ b/Carpet/CarpetRegrid2/src/regrid.cc
@@ -174,6 +174,68 @@ namespace CarpetRegrid2 {
   
   
   
+  void
+  get_physical_boundary (rvect & physical_lower,
+                         rvect & physical_upper,
+                         rvect & spacing)
+  {
+    rvect interior_lower, interior_upper;
+    rvect exterior_lower, exterior_upper;
+    if (CCTK_IsFunctionAliased ("MultiPatch_GetDomainSpecification")) {
+      assert (Carpet::map >= 0);
+      CCTK_INT const ierr = MultiPatch_GetDomainSpecification
+        (Carpet::map, dim,
+         & physical_lower[0], & physical_upper[0],
+         & interior_lower[0], & interior_upper[0],
+         & exterior_lower[0], & exterior_upper[0],
+         & spacing[0]);
+      assert (not ierr);
+    } else if (CCTK_IsFunctionAliased ("GetDomainSpecification")) {
+      CCTK_INT const ierr = GetDomainSpecification
+        (dim,
+         & physical_lower[0], & physical_upper[0],
+         & interior_lower[0], & interior_upper[0],
+         & exterior_lower[0], & exterior_upper[0],
+         & spacing[0]);
+      assert (not ierr);
+    } else {
+      assert (0);
+    }
+  }
+  
+  void
+  calculate_exterior_boundary (rvect const & physical_lower,
+                               rvect const & physical_upper,
+                               rvect       & exterior_lower,
+                               rvect       & exterior_upper,
+                               rvect const & spacing)
+  {
+    rvect interior_lower, interior_upper;
+    if (CCTK_IsFunctionAliased ("MultiPatch_ConvertFromPhysicalBoundary")) {
+      assert (Carpet::map >= 0);
+      CCTK_INT const ierr = MultiPatch_ConvertFromPhysicalBoundary
+        (Carpet::map, dim,
+         & physical_lower[0], & physical_upper[0],
+         & interior_lower[0], & interior_upper[0],
+         & exterior_lower[0], & exterior_upper[0],
+         & spacing[0]);
+      assert (not ierr);
+    } else if (CCTK_IsFunctionAliased ("ConvertFromPhysicalBoundary")) {
+      CCTK_INT const ierr =
+        ConvertFromPhysicalBoundary
+        (dim,
+         & physical_lower[0], & physical_upper[0],
+         & interior_lower[0], & interior_upper[0],
+         & exterior_lower[0], & exterior_upper[0],
+         & spacing[0]);
+      assert (not ierr);
+    } else {
+      assert (0);
+    }
+  }
+  
+  
+  
   extern "C" {
     CCTK_INT
     CarpetRegrid2_Regrid (CCTK_POINTER_TO_CONST const cctkGH_,
@@ -208,6 +270,7 @@ namespace CarpetRegrid2 {
 #warning "TODO: check this (for Carpet, and maybe also for CartGrid3D)"
 #warning "TODO: (the check that these two are consistent should be in Carpet)"
     
+#if 0
     typedef vect<vect<CCTK_INT,2>,3> jjvect;
     jjvect nboundaryzones;
     jjvect is_internal;
@@ -222,7 +285,9 @@ namespace CarpetRegrid2 {
          & shiftout[0][0]);
       assert (not ierr);
     }
+#endif
     
+#if 0
     rvect physical_lower, physical_upper;
     rvect interior_lower, interior_upper;
     rvect exterior_lower, exterior_upper;
@@ -236,10 +301,16 @@ namespace CarpetRegrid2 {
          & spacing[0]);
       assert (not ierr);
     }
+#endif
+    rvect physical_lower, physical_upper;
+    rvect spacing;
+    get_physical_boundary (physical_lower, physical_upper,
+                           spacing);
     
     // Adapt spacing for convergence level
     spacing *= ipow ((CCTK_REAL) Carpet::mgfact, Carpet::basemglevel);
     
+#if 0
     {
       CCTK_INT const ierr =
         ConvertFromPhysicalBoundary
@@ -250,6 +321,11 @@ namespace CarpetRegrid2 {
          & spacing[0]);
       assert (not ierr);
     }
+#endif
+    rvect exterior_lower, exterior_upper;
+    calculate_exterior_boundary (physical_lower, physical_upper,
+                                 exterior_lower, exterior_upper,
+                                 spacing);
     
     //
     // Calculate the union of the bounding boxes for all levels
@@ -274,38 +350,43 @@ namespace CarpetRegrid2 {
     }
     regions.at(0).normalize();
     
-    // Loop over all centres
-    for (int n = 0; n < num_centres; ++ n) {
-      centre_description centre (cctkGH, n);
-      if (centre.active) {
-        
-        // Loop over all levels for this centre
-        for (int rl = 1; rl < centre.num_levels; ++ rl) {
-          
-          // Calculate a bbox for this region
-          rvect const rmin = centre.position - centre.radius.at(rl);
-          rvect const rmax = centre.position + centre.radius.at(rl);
-          
-          // Convert to an integer bbox
-          ivect const istride = hh.baseextents.at(0).at(rl).stride();
-          ivect const imin =
-            rpos2ipos (rmin, origin, scale, hh, rl)
-            - boundary_shiftout * istride;
-          ivect const imax =
-            rpos2ipos1 (rmax, origin, scale, hh, rl)
-            + boundary_shiftout * istride;
-          
-          ibbox const region (imin, imax, istride);
+    // Refine only patch 0
+    if (Carpet::map == 0) {
+      // Loop over all centres
+      for (int n = 0; n < num_centres; ++ n) {
+        centre_description centre (cctkGH, n);
+        if (centre.active) {
           
-          // Add this region to the list of regions
-          if (static_cast <int> (regions.size()) < rl+1) regions.resize (rl+1);
-          regions.at(rl) |= region;
-          regions.at(rl).normalize();
+          // Loop over all levels for this centre
+          for (int rl = 1; rl < centre.num_levels; ++ rl) {
+            
+            // Calculate a bbox for this region
+            rvect const rmin = centre.position - centre.radius.at(rl);
+            rvect const rmax = centre.position + centre.radius.at(rl);
+            
+            // Convert to an integer bbox
+            ivect const istride = hh.baseextents.at(0).at(rl).stride();
+            ivect const imin =
+              rpos2ipos (rmin, origin, scale, hh, rl)
+              - boundary_shiftout * istride;
+            ivect const imax =
+              rpos2ipos1 (rmax, origin, scale, hh, rl)
+              + boundary_shiftout * istride;
+            
+            ibbox const region (imin, imax, istride);
+            
+            // Add this region to the list of regions
+            if (static_cast <int> (regions.size()) < rl+1) {
+              regions.resize (rl+1);
+            }
+            regions.at(rl) |= region;
+            regions.at(rl).normalize();
+            
+          } // for rl
           
-        } // for rl
-        
-      } // if centre is active
-    }   // for n
+        } // if centre is active
+      }   // for n
+    }     // if map==0
     
     
       
@@ -448,6 +529,7 @@ namespace CarpetRegrid2 {
       rvect const level_physical_lower = physical_lower;
       rvect const level_physical_upper = physical_upper;
       rvect const level_spacing = spacing / rvect (hh.reffacts.at(rl));
+#if 0
       rvect level_interior_lower, level_interior_upper;
       rvect level_exterior_lower, level_exterior_upper;
       {
@@ -460,6 +542,11 @@ namespace CarpetRegrid2 {
            & level_spacing[0]);
         assert (not ierr);
       }
+#endif
+      rvect level_exterior_lower, level_exterior_upper;
+      calculate_exterior_boundary (level_physical_lower, level_physical_upper,
+                                   level_exterior_lower, level_exterior_upper,
+                                   level_spacing);
       
       ivect const level_exterior_ilower =
         rpos2ipos (level_exterior_lower, origin, scale, hh, rl);
@@ -975,20 +1062,20 @@ namespace CarpetRegrid2 {
       for (int m = 0; m < maps; ++ m) {
         maxrl = max (maxrl, rls.at(m));
       }
+      // All maps must have the same number of levels
+      for (int m = 0; m < maps; ++ m) {
+        regsss.at(m).resize (maxrl);
+      }
       
       // Make multiprocessor aware
       for (int rl = 0; rl < maxrl; ++ rl) {
         vector <vector <region_t> > regss (maps);
         for (int m = 0; m < maps; ++ m) {
-          if (rl < rls.at(m)) {
-            regss.at(m) = regsss.at(m).at(rl);
-          }
+          regss.at(m) = regsss.at(m).at(rl);
         }
         Carpet::SplitRegionsMaps (cctkGH, regss);
         for (int m = 0; m < maps; ++ m) {
-          if (rl < rls.at(m)) {
-            regsss.at(m).at(rl) = regss.at(m);
-          }
+          regsss.at(m).at(rl) = regss.at(m);
         }
       } // for rl
       
diff --git a/Carpet/CarpetWeb/version-4.html b/Carpet/CarpetWeb/version-4.html
index ce6f9dc6c..ab6dbb2f0 100644
--- a/Carpet/CarpetWeb/version-4.html
+++ b/Carpet/CarpetWeb/version-4.html
@@ -44,16 +44,16 @@
       defines or updates the grid hierarchy, and reduces the chance of
       inconsistencies therein.</p></li>
 
-      <li></p>During checkpointing and recovery, the grid structure is
+      <li><p>During checkpointing and recovery, the grid structure is
       saved and restored by default.  This avoids accidental changes
       upon recovery.</p></li>
 
-      <li></p>The efficiency of I/O has been increased, especially for
+      <li><p>The efficiency of I/O has been increased, especially for
       HDF5 based binary I/O.  It is possible to combine several
       variables into one file to reduce the number of output
       files.</p></li>
 
-      <li></p>A new thorn LoopControl offers iterators over grid
+      <li><p>A new thorn LoopControl offers iterators over grid
       points, implemented as C-style macros.  These iterators allow
       additional important loop-level optimisations, such
       as <a href="http://en.wikipedia.org/wiki/Loop_tiling">loop
@@ -143,7 +143,7 @@
       reduces the number of output files.</p></li>
     </ul>
 
-    <h3><Communication, Mesh Refinement</h3>
+    <h3>Communication, Mesh Refinement</h3>
 
     <ul>
       <li><p>The algorithm determining the communication schedule is
diff --git a/Carpet/LoopControl/src/lc_hill.c b/Carpet/LoopControl/src/lc_hill.c
index f24b0efae..6ab44dba8 100644
--- a/Carpet/LoopControl/src/lc_hill.c
+++ b/Carpet/LoopControl/src/lc_hill.c
@@ -34,9 +34,9 @@ drand (void)
 
 static
 int
-irand (int const imax)
+irand (int const imaxval)
 {
-  return rand() / (RAND_MAX + 1.0) * imax;
+  return rand() / (RAND_MAX + 1.0) * imaxval;
 }
 
 
diff --git a/Carpet/LoopControl/src/loopcontrol.c b/Carpet/LoopControl/src/loopcontrol.c
index 6f420e0c1..c03cbf452 100644
--- a/Carpet/LoopControl/src/loopcontrol.c
+++ b/Carpet/LoopControl/src/loopcontrol.c
@@ -295,10 +295,8 @@ lc_stattime_init (lc_stattime_t * restrict const lt,
   
   
   /* Append to loop statistics list */
-/*   _Pragma ("omp critical") { */
-    lt->next = ls->stattime_list;
-    ls->stattime_list = lt;
-/*   } */
+  lt->next = ls->stattime_list;
+  ls->stattime_list = lt;
 }
 
 lc_stattime_t *
@@ -450,10 +448,8 @@ lc_statset_init (lc_statset_t * restrict const ls,
   ls->time_calc_sum2  = 0.0;
   
   /* Append to loop statistics list */
-/*   _Pragma ("omp critical") { */
-    ls->next = lm->statset_list;
-    lm->statset_list = ls;
-/*   } */
+  ls->next = lm->statset_list;
+  lm->statset_list = ls;
 }
 
 lc_statset_t *
@@ -510,21 +506,34 @@ lc_statset_find_create (lc_statmap_t * restrict const lm,
 
 /* Initialise loop statistics */
 void
-lc_statmap_init (lc_statmap_t * restrict const lm,
+lc_statmap_init (int * restrict const initialised,
+                 lc_statmap_t * restrict const lm,
                  char const * restrict const name)
 {
   /* Check arguments */
+  assert (initialised);
   assert (lm);
   
-  /* Set name */
-  lm->name = strdup (name);
-  
-  /* Initialise list */
-  lm->statset_list = NULL;
+#pragma omp single
+  {
+    
+    /* Set name */
+    lm->name = strdup (name);
+    
+    /* Initialise list */
+    lm->statset_list = NULL;
+    
+    /* Append to loop statistics list */
+    lm->next = lc_statmap_list;
+    lc_statmap_list = lm;
+    
+  }
   
-  /* Append to loop statistics list */
-  lm->next = lc_statmap_list;
-  lc_statmap_list = lm;
+#pragma omp single
+  {
+    /* Set this flag only after initialising */
+    * initialised = 1;
+  }
 }
 
 
@@ -563,7 +572,8 @@ lc_control_init (lc_control_t * restrict const lc,
   
   
   lc_statset_t * restrict ls;
-  _Pragma ("omp single copyprivate (ls)") {
+#pragma omp single copyprivate (ls)
+  {
     /* Get number of threads */
     int const num_threads = omp_get_num_threads();
     
@@ -579,7 +589,8 @@ lc_control_init (lc_control_t * restrict const lc,
   
   
   lc_stattime_t * restrict lt;
-  _Pragma ("omp single copyprivate (lt)") {
+#pragma omp single copyprivate (lt)
+  {
     
     lc_state_t state;
     
@@ -778,7 +789,8 @@ lc_control_finish (lc_control_t * restrict const lc)
   lc_statset_t * restrict const ls = lc->statset;
   
   int ignore_iteration;
-  _Pragma ("omp single copyprivate (ignore_iteration)") {
+#pragma omp single copyprivate (ignore_iteration)
+  {
     DECLARE_CCTK_PARAMETERS;
     ignore_iteration = ignore_initial_overhead && lt->time_count == 0.0;
   }
@@ -798,7 +810,8 @@ lc_control_finish (lc_control_t * restrict const lc)
   double const time_calc_sum2 = pow (time_calc_sum, 2);
 
   /* Update statistics */
-  _Pragma ("omp critical") {
+#pragma omp critical
+  {
     lt->time_count += 1.0;
     
     lt->time_setup_sum  += time_setup_sum;
@@ -816,21 +829,24 @@ lc_control_finish (lc_control_t * restrict const lc)
     ls->time_calc_sum2 += time_calc_sum2;
   }
   
-  _Pragma ("omp master") {
+#pragma omp master
+  {
     lt->last_updated = time_calc_end;
   }
   
-  _Pragma ("omp barrier");
+#pragma omp barrier
   
   {
     DECLARE_CCTK_PARAMETERS;
     if (use_simulated_annealing) {
-      _Pragma ("omp single") {
+#pragma omp single
+      {
         lc_auto_finish (ls, lt);
       }
     }
     if (use_random_restart_hill_climbing) {
-      _Pragma ("omp single") {
+#pragma omp single
+      {
         lc_hill_finish (ls, lt);
       }
     }
@@ -912,11 +928,12 @@ lc_printstats (CCTK_ARGUMENTS)
 
 CCTK_FCALL
 void
-CCTK_FNAME (lc_statmap_init) (lc_statmap_t * restrict const lm,
+CCTK_FNAME (lc_statmap_init) (int * restrict const initialised,
+                              lc_statmap_t * restrict const lm,
                               ONE_FORTSTRING_ARG)
 {
   ONE_FORTSTRING_CREATE (name);
-  lc_statmap_init (lm, name);
+  lc_statmap_init (initialised, lm, name);
   free (name);
 }
 
diff --git a/Carpet/LoopControl/src/loopcontrol.h b/Carpet/LoopControl/src/loopcontrol.h
index 5d640c3b6..17eb8a8cc 100644
--- a/Carpet/LoopControl/src/loopcontrol.h
+++ b/Carpet/LoopControl/src/loopcontrol.h
@@ -10,6 +10,10 @@
 
 #include <cctk_Arguments.h>
 
+#ifdef __cplusplus
+extern "C" {
+#endif
+
 
 
 #if 0
@@ -247,7 +251,8 @@ lc_max (int const i, int const j)
 
 
 void
-lc_statmap_init (lc_statmap_t * restrict ls,
+lc_statmap_init (int * restrict initialised,
+                 lc_statmap_t * restrict ls,
                  char const * restrict name);
 
 void
@@ -264,16 +269,10 @@ lc_control_finish (lc_control_t * restrict lc);
 
 #define LC_LOOP3(name, i,j,k, imin,jmin,kmin, imax,jmax,kmax, ilsh,jlsh,klsh) \
   do {                                                                  \
-    static lc_statmap_t lc_lm;                                          \
     static int lc_initialised = 0;                                      \
+    static lc_statmap_t lc_lm;                                          \
     if (! lc_initialised) {                                             \
-      _Pragma ("omp single") {                                          \
-        lc_statmap_init (& lc_lm, #name);                               \
-      }                                                                 \
-      _Pragma ("omp single") {                                          \
-        /* Set this flag only after initialising */                     \
-        lc_initialised = 1;                                             \
-      }                                                                 \
+      lc_statmap_init (& lc_initialised, & lc_lm, #name);               \
     }                                                                   \
     lc_control_t lc_lc;                                                 \
     lc_control_init (& lc_lc, & lc_lm,                                  \
@@ -315,6 +314,10 @@ lc_control_finish (lc_control_t * restrict lc);
 void
 lc_printstats (CCTK_ARGUMENTS);
 
+#ifdef __cplusplus
+}
+#endif
+
 #endif
 
 
diff --git a/Carpet/doc/domain-decomposition.tex b/Carpet/doc/domain-decomposition.tex
index 97252e87f..f494c03ca 100644
--- a/Carpet/doc/domain-decomposition.tex
+++ b/Carpet/doc/domain-decomposition.tex
@@ -1,86 +1,128 @@
-\documentclass[nofootinbib, twocolumn]{revtex4}
+\documentclass[oneside]{amsart}
 
+\usepackage{amsmath}
 \usepackage{amssymb}
 \usepackage{color}
+\usepackage{graphicx}
+\usepackage[latin9]{inputenc}
 \usepackage{mathpazo}
+\usepackage[hyphens]{url}
+
+% Put this package last
+\usepackage[bookmarks, bookmarksopen, bookmarksnumbered]{hyperref}
+% Put this package after hyperref
+\usepackage[all]{hypcap}
+% Don't use tt font for urls
+\urlstyle{rm}
+
+
 
 % Make a comment stand out visually
 \newcommand{\todo}[1]{{\color{blue}$\blacksquare$~\textsf{[TODO: #1]}}}
 
+% Name of a code
+\newcommand{\code}[1]{\texttt{#1}}
+
+
+
+\hyphenation{Cac-tus-Ein-stein Schwarz-schild South-amp-ton}
+
+\sloppypar
+
+
+
 \begin{document}
 
-\title{Domain Decomposition in Carpet:\\Regridding and Determination
-  of the Communication Schedule}
+\title[Domain Decomposition in Carpet]{Domain Decomposition in
+  Carpet:\\Regridding and Determination of the Communication Schedule}
 
-\author{Erik Schnetter}
+\author{Erik Schnetter$^{(1,2)}$}
 
-\email{schnetter@cct.lsu.edu}
+% \email
+\thanks{\emph{Email}:~\url{mailto:schnetter@cct.lsu.edu}}
 
-\affiliation{Center for Computation \& Technology, 216 Johnston Hall,
-  Louisiana State University, Baton Rouge, LA 70803, USA}
+% \urladdr
+\thanks{\emph{Web}:~\url{http://www.cct.lsu.edu/}}
 
-\homepage{http://www.cct.lsu.edu/about/focus/numerical}
+% \address
+\thanks{{}$^{(1)}$ Center for Computation \& Technology,
+  Louisiana State University, Baton Rouge, LA, USA}
 
-\affiliation{Department of Physics and Astronomy, 202 Nicholson Hall,
-  Louisiana State University, Baton Rouge, LA 70803, USA}
+\thanks{{}$^{(2)}$ Department of Physics and Astronomy,
+  Louisiana State University, Baton Rouge LA, USA}
 
-\homepage{http://relativity.phys.lsu.edu/}
 
-\date{March 24, 2007}
+
+\date{March 9, 2008}
 
 
 
 \begin{abstract}
-  This paper describes the algorithms that Carpet
+  This text describes the algorithms that Carpet
   \cite{Schnetter-etal-03b, carpetweb} uses for regridding, domain
-  decomposition, and to set up the communication schedule.
+  decomposition, and setting up the communication schedule.  It is
+  intended for readers interested more details than a casual user
+  would need.  This text explains the concepts that Carpet uses
+  internally to describe the grid hierarchy and to ensure its
+  consistency.
 \end{abstract}
 
 \maketitle
 
-\todo{Submit this paper to Tom}
-
 
 
 \section{Introduction}
 
-Regridding in Carpet is handled by three entities: A \emph{regridding
-  thorn} which is responsible for deciding on a grid hierarchy and for
-the domain decomposition, \emph{Carpet} itself which decides the
-extent of the domain and the type of outer boundary conditions, and
-\emph{CarpetLib} which handles the details.
+Setting up a grid hierarch (``regridding'') is in Carpet handled by
+three different entities: \emph{Carpet} itself decides the extent of
+the domain, the type of outer boundary conditions, and distributes the
+domain onto processors; a \emph{regridding thorn} is responsible for
+deciding the shape of the grid hierarchy, and \emph{CarpetLib} handles
+the details and actually manages the data.  (Technically speaking, it
+is the regridding thorn which has to do the domain decomposition;
+however, it can simply call a convenient helper routine in Carpet for
+this task.)
 
-In the following we concentrated on a single patch which contains a
-mesh refinement hierarchy.  If there are multiple patches, then each
-patch is conceptually handled independently.  Certain conditions may
-have to be satisfied if a refined mesh touches an inter-patch
-boundary, but these is not discussed here.
+This separation leaves the decision on the shape of the grid hiearchy
+to a thorn which can be replaced or augmented as necessary.  All
+handling of data happens in CarpetLib, which is thus the only entity
+which needs to be optimised for speed.  Finally, Carpet retains the
+overview over the regridding process.
 
-We assume that the domain has $D$ spatial dimensions.  Usually $D=3$.
+In the following we assume that there is a single patch (block) which
+contains a mesh refinement hierarchy.  If there are multiple patches,
+then each patch is conceptually handled independently.  Certain
+conditions may have to be satisfied if a refined mesh touches an
+inter-patch boundary, but these are not discussed here.
+
+We assume that the domain has $D$ spatial dimensions, usually $D=3$.
 
 
 
 \section{Domain description}
 
+We assume that boundary location and boundary discretisation are set
+up via \emph{CoordBase}.  This is necessary since other methods do not
+allow specifying sufficient details to handle e.g.\ refined regions
+intersecting mesh refinement boundaries.  Below we give an overview
+over the information that needs to be specified to describe a
+boundary.  See the CoordBase thorn guide for details.
+
 The extent of the overall simulation domain is described in terms or
 real-valued coordinates.  The domain is cuboidal, i.e., it can be
 described in terms of two vectors $x^i_\mathrm{min}$ and
 $x^i_\mathrm{max}$.
 
-The thorn \textrm{CoordBase} has facilities to describe the domain
-extent in this manner.  Thorn \textrm{CartGrid3D} can also be used
-instead.
-
 The type of boundary conditions for each of the $2D$ faces is
 described by three quantities:
 \begin{itemize}
-  \item the total number of boundary points,
-  \item how many of these are outside of the domain boundary,
-  \item whether the boundary is staggered with respect to the domain
-    boundary, or whether one of the boundary points is located at the
-    domain boundary.
+\item the total number of boundary points,
+\item how many of these points are outside of the domain boundary,
+\item whether the boundary is staggered with respect to the domain
+  boundary, or whether one of the boundary points is located on the
+  domain boundary.
 \end{itemize}
-\todo{Show figure.}
+\todo{Show the relevant figures from the CoordBase thorn guide.}
 
 The extent of the domain $L^i := x^i_\mathrm{max} - x^i_\mathrm{min}$
 must be commensurate with the coarse grid spacing $h^i$.  This means
@@ -88,28 +130,36 @@ that $L^i$ must be a multiple of $h^i$, modulo the fact that the
 boundary may be staggered.
 
 An \emph{outer boundary condition} can either be a \emph{physical}
-boundary condition, such as e.g.\ a Sommerfeld condition, or a
+boundary condition, such as e.g.\ a Direchlet condition, or a
 \emph{symmetry} boundary condition, such as e.g.\ a reflection
 symmetry.  This distinction is irrelevant for Carpet.  Both kinds of
 outer boundary conditions are applied by thorns which are scheduled in
 the appropriate way.
 
-Note that Carpet does currently not support evolving outer boundary
-points in time if there is a mesh refinement boundary touching the
-refined region.  Carpet cannot fill the overlap of the outer boundary
-and the mesh refinement boundary through interpolation, because this
-requires off-centred interpolation stencils, which are not implemented
-at the moment.  Therefore Carpet does not fill any outer boundary
-points through interpolation.
-  
-Carpet also does not fill outer boundary points through
-synchronisation.  The outer boundary condition is supposed to handle
-these.  This requires that the outer boundary condition is apply after
-synchronisation.  This is usally automatically the case when using the
-Cactus boundary condition selection mechanism, i.e., when the group
-\texttt{ApplyBCs} is scheduled.  Synchronising outer boundary points
-would be possible, but is not done so that refinement boundaries and
-inter-processor boundaries are treated in the same way.
+The main distinction between an \emph{outer boundary point} and an
+\emph{interior point} from Carpet's point of view is that an outer
+boundary point is not evolved in time.  Instead, the value of boundary
+points must be completely determined by the value of interior points.
+(This is clearly the case for Direchlet or symmetry boundary
+conditions.)  The reason for this distinction is that Carpet cannot
+fill outer boundary points on refined grids via interpolation, because
+this requires off-centred interpolation stencils which are not
+implemented at the moment.  Therefore, Carpet does not fill any outer
+boundary points through interpolation.
+
+\todo{We should introduce mesh refinement and parallelisation before
+  this paragraph.}
+%
+If the refined region abuts the outer boundary, then outer boundary
+points can also be refinement boundary points.  Carpet does not fill
+these outer boundary points through synchronisation (prolongation)
+either.  This requires that the outer boundary condition is applied
+after every synchronisation.  This is usally automatically the case
+when using the Cactus boundary condition selection mechanism, i.e.,
+when the group \texttt{ApplyBCs} is scheduled.  Synchronising outer
+boundary points would be possible, but this is not done so that
+refinement boundaries and inter-processor boundaries are treated in
+the same way.
 
 
 
@@ -120,41 +170,46 @@ A regridding thorn, such as e.g.\ \texttt{CarpetRegrid} or
 hierarchy consists of several \emph{refinement levels}, and each
 refinement level consists of several \emph{refined regions}.  A
 refined region is a cuboid set of grid points which is assigned to a
-particular processor.  The refined regions on one level may not
-overlap.  (If they do not touch each other or the outer boundary, then
-they usually have to keep a certain minimum distance, as described
-below.)
-
-Refined regions are described by the data structure
-\texttt{region\_t}, which is declared in
-\texttt{CarpetLib/src/region.hh}.
-
-The semantics of a grid hierarchy is independent of the number of
-refined regions, or of the processors which are assigned to them.  The
-only relevant quantity is the set of refined grid points, i.e., the
-conjunction of all refined regions.  For example, when running on more
-processors, then regions will be split up so that there are more and
-smaller regions.  There can be more than one region per processor for
-each level.
+particular processor.  (The ``refined regions'' which are specified
+e.g.\ in a parameter file are broken up during domain is decomposition
+into a set of regions, so that each region is assigned to exactly one
+processor.)  There can be zero or multiple regions per processor, and
+differnet processors may own different numbers of regions.  Each face
+of a region is either an outer or an internal boundary.  Refined
+regions are described by the data structure \code{region\_t}, which is
+declared in \code{CarpetLib/src/region.hh}.
+
+The refined regions on one level may not overlap.  (If they do not
+touch each other or the outer boundary, then they usually even have to
+keep a certain minimum distance, as described below.)
+
+The ``semantics'' (the result obtained in a simulation) of a grid
+hierarchy is independent of the number of refined regions, or of the
+processors which are assigned to them.  The only relevant quantity is
+the set of refined grid points, i.e., the conjunction of all refined
+regions.  When running on more processors, then regions will be split
+up so that there are more and smaller regions.
 
 The assignment of regions to processors is handled during regridding
 because it affects performance.  When there are only few processors
 available, then it may be more efficient to use fewer regions.
 Combining regions may require some fill-in.  \todo{Show figure.}  No
-current regridding thorn in Carpet uses this at the moment; currently,
-the set of refined points is independent of the number of processors.
-
-Each face of a refined region can be an outer boundary face.  This
-means that Carpet will add no ghost or buffer zones there, and will
-mark this face as outer boundary when calling thorns.  An outer
-boundary face must extend up to the outermost outer boundary point,
-since otherwise thorns will apply the outer boundary conditions
-incorrectly.  The regridding thorn has to ensure this.  Faces which
-are not outer boundary faces must be sufficiently far away from outer
-boundaries, so that any ghost or buffer zones which are added to that
-face do not intersect the outer boundary.  \todo{Show figure.}
-
-Each face which is not an outer boundary face can be either an
+current regridding thorn in Carpet offers this at the moment;
+currently, the set of refined points is independent of the number of
+processors.
+
+Each face of a refined region is either an interior or an outer
+boundary face.  Carpet adds no ghost or buffer zones to outer boundary
+faces, and marks this face as outer boundary when calling thorns.  An
+outer boundary face must extend up to the outermost outer boundary
+point, since otherwise thorns will apply the outer boundary conditions
+incorrectly.  The regridding thorn has to ensure this property.  Faces
+which are not outer boundary faces must be sufficiently far away from
+outer boundaries, so that any ghost or buffer zones which are added to
+that face do not intersect the outer boundary.  \todo{Show figure.}
+The regridding thorn also has to ensure this property.
+
+Each face which is not an outer boundary face is either an
 inter-processor or a refinement boundary face.  Inter-processor faces
 have no buffer zones, and the ghost zones can be filled in completely
 from other refined regions on the same level.  Refinement boundary
@@ -163,22 +218,27 @@ be filled via prolongation from the next coarser grid.  It is possible
 to have faces which are partly an inter-processor boundary and partly
 a refinement boundary.  These are counted as and handled as refinement
 boundaries, although some of the ghost zones may still be filled via
-synchronisation.  This is explained in more detail below.
+synchronisation.  (Any grid points which can be filled via
+synchronisation are always also filled via synchronisation.)  This is
+explained in more detail below.
 
-As described below, ghost zones are added at the outside of grids by
+As described below, ghost zones are added at the outside of regions by
 Carpet.  Buffer zones need to be taken into account by the regridding
 thorn, most likely by extending the regined regions appropriately.
-(In terms of previous terminology: buffer zones are ``inner'' buffer
-zones; ``outer'' buffer zones do not exist any more.  However, since
-they are added by the regridding thorn, they do not need to be taken
-into account when specifying the refinement hierarchy.)  Carpet places
-buffer zones at all faces which are marked as refinement boundaries.
+(In terms of previous terminology: buffer zones are always ``inner''
+buffer zones; ``outer'' buffer zones do not exist any more.  However,
+if they are added by the regridding thorn, they do not need to be
+taken into account when specifying the refinement hierarchy in a
+parameter file -- this depends on the regridding thorn.
+\todo{CarpetRegrid2 does this, CarpetRegrid does this not.})  Carpet
+assumes buffer zones at all faces which are marked as refinement
+boundaries.
 
 It is the responsibility of the regridding thorn to mark outer
 boundaries and refinement boundaries appropriately.  If the outer
-boundary mark is chosen wrong, then the simulation is inconsistent,
+boundary mark is chosen wrongly, then the simulation is inconsistent,
 since the outer boundary condition may be applied at a the wrong
-location, or the outer boundary may be filled by interpolation, which
+location, or the outer boundary may be filled by interpolation which
 is less accurate.  Depending on the number of boundary points and
 stencil sizes, this may or may not be detected later by
 self-consistency checks.
@@ -187,18 +247,31 @@ self-consistency checks.
 
 \section{Zoning}
 
-Consider a single refinement level.
+This section defines the algorithm which Carpet uses to define which
+grid points are defined by what action.  This algorithm is codified in
+\code{dh.cc}.  Since the code in \code{dh.cc} is tested, it should be
+assumed to be correct where it differs with this description.  This
+algorithm is applied to every refinement level.
+
+\emph{Grid arrays} are handled in the same way as grid functions,
+except that there are no refined regions.  This may seem like overkill
+at first, but in fact it greatly simplifies the implementation since
+grid arrays have (almost) only a subset of the capabilities of grid
+functions.
 
 \subsection{Domain}
 
 Let $dom$ be the simulation domain.  Let $domact$ be the set of grid
-points which are evolved in time, which is required to be cuboidal and
-not empty.  (An empty region is also counted as cuboidal.)
+points which are evolved in time, and which is required to be cuboidal
+and not empty.\footnote{While grid functions cannot be empty, $domact$
+  for grid arrays can be empty.}  (An empty region is also counted as
+cuboidal.)
 
 Each face has a layer of outer boundary points.  Let $domob$ be the
 set of these layers.  We require that all values in $domob$ can be
-calculated from the values in $domact$.  Let $domext :=domact \cup
-domob$.  It follows that $domext$ is cuboidal and not empty.
+calculated from the values in $domact$.  Let $domext := domact \cup
+domob$.  It follows that $domext$ is cuboidal and not
+empty.
 
 \todo{Show figure.}
 
@@ -206,11 +279,11 @@ The following properties hold for $dom$:
 \begin{itemize}
 \item $domact$ is cuboidal and not empty
 \item $domob$ is a possibly empty layer around $domact$
-\item $domob$ has a width of \texttt{nboundaryzones} as obtained from
-  \texttt{GetBoundarySpecification}
+\item $domob$ has a width of \code{nboundaryzones} as obtained from
+  \code{GetBoundarySpecification}
 \item $domact \cap domob = \emptyset$
 \item $domext = domact \cup domob$ is cuboidal and not empty
-\item $domext$ has the size \texttt{cctk\_gsh}
+\item $domext$ has the size \code{cctk\_gsh}
 \end{itemize}
 For completeness, one can define $dombnd = dombuf = \emptyset$, and
 $domint = domown = domact$.  Then the same relations hold for $dom$ as
@@ -226,69 +299,80 @@ There may be other refined regions $reg'$ on the same level.  We
 require that $\forall_{reg' \ne reg}: regint \cap regint' =
 \emptyset$.
 
-Each face of $regown$ which is not an outer boundary face has a layer
-of ghost zones.  Let $regghost$ be the set of these layers.  We
-require that $regghost \subseteq domact$.  Let $regext := regint \cup
-regghost$.  It follows that $regext$ is cuboidal and not empty.  We
-require that $regext \subseteq domext$.
+Each face of $regint$ which is not an outer boundary face has a layer
+of \code{cctk\_nghostzones} ghost zones added to it.  (These are what
+Cactus calls ``ghost zones'', but note that not all ghost zones are
+filled via synchronisation.)  Let $regghost$ be the set of these
+layers.  We require that $regghost \subseteq domact$.  Let $regext :=
+regint \cup regghost$ be the interior with the ghost zones added.  It
+follows that $regext$ is cuboidal and not empty.  We require that
+$regext \subseteq domext$.
 
 \todo{Show figure.}
 
-Let $regown := regint - regbnd$.  Then $regown$ is cuboidal, and we
-require it to be not empty.
-
-Each face of $regown$ which is not an outer boundary face has a layer
-of ghost zones.  Let $regbnd$ be the set of these layers.  Let
-$regouter := regint \cup regbnd$.  It follows that $regouter$ is
-cuboidal and not empty.  We require that $regouter \subseteq domact$.
-It follows that $regouter \subseteq regext$.
-
-Let $regob := regext - regouter$.  It follows that $regob \subseteq
-domob$.
-
-Wherever a face which is not an outer boundary face does not touch
-another refined region, buffer zones exist in $regown$.  Let $allown
-:= \bigcup regown$.  Let $allbuf$ be the layer of outermost grid
-points of $allown$ on those faces which do not touch the outer
-boundary.  Then $regbuf := allbuf \cup regown$.  It follows that
-$regbuf \subseteq domact$.  Let $regact := regown - rebuf$.  It
-follows that $regact$ is cuboidal.  (Should we required that $regact$
-be not empty?)
+On each face of $regext$, the outermost layer of \code{nboundaryzones}
+points are not communicated.  Let $regcomm$ be the set of communicated
+points.  The $regcomm \subseteq regext$, we require it to be not
+empty, and we require that $recomm \subseteq domact$.
+
+Let $regob := regext - regcomm$ be the set of outer boundary points.
+It follows that $regob \subseteq domob$.
+
+The owned region $regown$ is the interior region, extended up to the
+boundary, i.e., $regint$ with a layer of \code{nboundaryzones} points
+added to it at all outer boundary faces.  Thus $regown$ is cuboidal,
+and we require it to be not empty.  It is also $regown \subseteq
+domact$, and $regown \subseteq regext$.  It is also $\forall_{reg' \ne
+  reg}: regown \cap regown = \emptyset$.
+
+Let $regbnd := regcomm - regown$.  Then $regbnd \subseteq domact$.
+\todo{Does this follow, or is this a requirement?}  These points are
+filled via synchronisation \todo{Is this so?  Is this always via
+  inter-processor communiocation, or also via prolongation?}, and they
+cannot be too close to the outer boundary.
+
+Let $regbuf$ be the set of grid points which have a distance less than
+\code{buffer\_width} from the boundary of the conjunction of all
+active grid points.  Then $regbuf \subseteq regown$.  Let $regact :=
+regown - regbuf$.  Then also $regact \subseteq regown$.  In general,
+$regact$ is not cuboidal.  $regact$ can be empty.  \todo{This may be
+  controversial.}
 
 \todo{Show figure.}
 
 Let $regsync := regbnd \cap allown$ be the boundary points which can
 be filled via synchronisation.  Note that $regbnd \cap regown =
 \emptyset$ by construction, which means that a region cannot
-synchronise from itself.
+synchronise from itself.  Note that outer boundary points are
+synchronised for backward compatibility.
 
 Let $regref := (regbnd - regsync) \cup regbuf$ be the set of all
-points which need to be filled via prolongation.  Note that now
-$regref \cap regsync = \emptyset$, which means that no point is both
-synchronised and prolongated.  Note also that $regsync \cap regob =
-\emptyset$ and $regref \cap regob = \emptyset$, which means that no
-outer boundary point is synchronised or prolongated.
+points which need to be filled via prolongation.  Note that $regref
+\cap regsync = \emptyset$, which means that no point is both
+synchronised and prolongated.  Note also that $regref \cap regob =
+\emptyset$, which means that no outer boundary point is prolongated.
 
 \todo{Show figure.}
 
 The following properties hold for $reg$:
 \begin{itemize}
-\item $regact$ is cuboidal
-\item $regbuf$ is a layer around $regact$ on those faces which are
-  marked ``refinement boundary''
-\item $regown = regact \cup regbuf$ is cuboidal and not empty
+\item $regint \subseteq domext$ is cuboidal and not empty
+\item $regghost$ is a layer on the outside of $regint$ on those faces
+  which are not marked ``outer boundary''
+\item $regext = regint \cup regghost \subseteq domext$ is cuboidal and
+  not empty
+\item $recgomm \subseteq domact$
+\item $regob \subseteq domob$
+\item $regown \subseteq domact$ is cuboidal and not empty
 \item $\forall_{reg' \ne reg}: regown \cup regown' = \emptyset$
 \item $regbnd$ is a layer around $regown$ on those faces which are not
   marked ``outer boundary''
-\item $regouter = regown \cup regbnd$ is cuboidal and not empty
-\item $regouter \subseteq domact$
-\item $regob \subseteq domob$
-\item $regext = regint \cup regob$ is cuboidal and not empty
-\item $regext \subseteq domext$
-
-\item $regghost$ is a layer on the outside of $regext$ on those faces
-  which are not marked ``outer boundary''
-\item $regint = regext - regghost$ is cuboidal and not empty
+\item $regbuf$ is a layer around $regact$ on those faces which are
+  marked ``refinement boundary''
+\item $regsync$ are those boundary points which can be filled via
+  synchronisation
+\item $regref$ are those boundary or buffer points which are filled
+  via prolongation
 \end{itemize}
 $regint$ is not important for Carpet's working, but only for Cactus.
 Since Cactus has no notion of outer boundaries, we introduce $regint$,
@@ -345,7 +429,7 @@ number of buffer zones.  Then $allbuf := allown - notact$.
 
 
 % With titles and urls
-\bibliographystyle{bibtex/apsrev-titles}
+\bibliographystyle{bibtex/amsplain-url}
 \bibliography{bibtex/references}
 
 \end{document}
diff --git a/Carpet/doc/scheduling.tex b/Carpet/doc/scheduling.tex
index e148fbfea..5d298c0e3 100644
--- a/Carpet/doc/scheduling.tex
+++ b/Carpet/doc/scheduling.tex
@@ -639,20 +639,25 @@ PARAMCHECK
 \fbox{\begin{minipage}[t]{\textwidth}
 Recover? (yes/no)
 \\
+\fbox{\begin{minipage}[t]{0.47\textwidth}
+PREREGRIDINITIAL\\
+Regrid\\
+POSTREGRIDINITIAL
+\end{minipage}}
+\fbox{\begin{minipage}[t]{0.47\textwidth}
+Recover grid structure
+\end{minipage}}\\
 \fbox{\begin{minipage}[t]{0.1\textwidth}
 initial\\
 loop
 \end{minipage}
 \fbox{\begin{minipage}[t]{0.345\textwidth}
-PREREGRIDINITIAL\\
-Regrid\\
-POSTREGRIDINITIAL\\
 BASEGRID\\
 INITIAL\\
+POSTINITIAL\\
 $\longrightarrow$ Recurse\\
 Restrict\\
 POSTRESTRICTINITIAL\\
-POSTINITIAL\\
 POSTSTEP
 \end{minipage}}}
 \fbox{\begin{minipage}[t]{0.1\textwidth}
@@ -784,6 +789,7 @@ SHUTDOWN
         \begin{minipage}[c]{0.6\textwidth}
          BASEGRID\\
          INITIAL
+         POSTINITIAL\\
         \end{minipage}
        }
       }
@@ -846,7 +852,6 @@ SHUTDOWN
      \fbox{
       \begin{minipage}[c]{0.49\textwidth}
        POSTRESTRICTINITIAL\\
-       POSTINITIAL\\
        POSTSTEP
       \end{minipage}
      }