Enhance support for multi-patch simulations

Carpet:

Ensure that at most one of GetDomainSpecificatio or
MultiPatch_GetDomainSpecification is defined.

Allow the boundary case of having zero regions on a refinement level.

CarpetLib:

Allow the boundary case of having zero components on a patch, if there
are still more than zero components overall.

CarpetIOASCII:

Loop over the components using explicit light-weight for loops instead of
Carpet's looping macros.  This is faster, since less state information
needs to be updated.

Correct an inconsistency in converting between integer indices and
coordinates in multi-patch simulations.

CarpetRegrid2:

Take multi-patch systems into account.

7 files changed, 266 insertions, 117 deletions

diff --git a/Carpet/Carpet/src/Recompose.cc b/Carpet/Carpet/src/Recompose.cc
index 3027290de..4da7afb73 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();
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/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/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