Merge /Users/eschnett/Cbeta/carpet

35 files changed, 1704 insertions, 547 deletions

diff --git a/Carpet/Carpet/interface.ccl b/Carpet/Carpet/interface.ccl
index 734020bae..cee9c1075 100644
--- a/Carpet/Carpet/interface.ccl
+++ b/Carpet/Carpet/interface.ccl
@@ -12,10 +12,12 @@ uses include header: dist.hh
 
 uses include header: bbox.hh
 uses include header: bboxset.hh
+uses include header: fulltree.hh
 uses include header: region.hh
 uses include header: vect.hh
 
 uses include header: data.hh
+
 uses include header: dh.hh
 uses include header: gf.hh
 uses include header: ggf.hh
@@ -205,18 +207,22 @@ PROVIDES FUNCTION GetRefinementLevels \
 
 # The true prototype of the routine below:
 #   int Carpet_Regrid (const cGH * cctkGH,
+#                      gh::rregs * superregss,
 #                      gh::mregs * regsss,
 #                      int         force);
 CCTK_INT FUNCTION Carpet_Regrid (CCTK_POINTER_TO_CONST IN cctkGH, \
+                                 CCTK_POINTER IN superregss,      \
                                  CCTK_POINTER IN regsss,          \
                                  CCTK_INT IN force)
 USES FUNCTION Carpet_Regrid
 
 # The true prototype of the routine below:
 #   int Carpet_Regrid (const cGH * cctkGH,
+#                      vector<gh::rregs> * superregsss,
 #                      vector<gh::mregs> * regssss,
 #                      int                  force);
 CCTK_INT FUNCTION Carpet_RegridMaps (CCTK_POINTER_TO_CONST IN cctkGH, \
+                                     CCTK_POINTER IN superregsss,     \
                                      CCTK_POINTER IN regssss,         \
                                      CCTK_INT IN force)
 USES FUNCTION Carpet_RegridMaps
diff --git a/Carpet/Carpet/src/Recompose.cc b/Carpet/Carpet/src/Recompose.cc
index 9bd7d4066..696adf7ca 100644
--- a/Carpet/Carpet/src/Recompose.cc
+++ b/Carpet/Carpet/src/Recompose.cc
@@ -60,12 +60,14 @@ namespace Carpet {
   
   
   static void
-  SplitRegionsMaps_Automatic_Recursively (bvect    const & dims,
-                                          int      const   nprocs,
-                                          region_t const & reg,
+  SplitRegionsMaps_Automatic_Recursively (bvect    const   & dims,
+                                          int      const     firstproc,
+                                          int      const     nprocs,
+                                          region_t         & superreg,
                                           vector<region_t> & newregs);
   static void
   SplitRegions_AsSpecified (cGH const * cctkGH,
+                            vector<region_t> & superregs,
                             vector<region_t> & regs);
   
   
@@ -124,6 +126,8 @@ namespace Carpet {
   {
     DECLARE_CCTK_PARAMETERS;
     
+    Checkpoint ("Regridding level %d...", reflevel);
+    
     assert (is_level_mode());
     
     bool const have_regrid     = CCTK_IsFunctionAliased ("Carpet_Regrid");
@@ -155,37 +159,40 @@ namespace Carpet {
       
       BEGIN_MAP_LOOP(cctkGH, CCTK_GF) {
         
-        gh::mregs regsss = vhh.at(map)->regions;
+        gh::rregs superregss = vhh.at(map)->superregions;
+        gh::mregs regsss;
         
         // Check whether to recompose
         CCTK_INT const do_recompose =
-          Carpet_Regrid (cctkGH, & regsss, force_recompose);
+          Carpet_Regrid (cctkGH, & superregss, & regsss, force_recompose);
         assert (do_recompose >= 0);
         did_change = did_change or do_recompose;
         
         if (do_recompose) {
-          RegridMap (cctkGH, map, regsss);
+          RegridMap (cctkGH, map, superregss, regsss);
         }
         
       } END_MAP_LOOP;
       
     } else {
       
+      vector<gh::rregs> superregsss (maps);
       vector<gh::mregs> regssss (maps);
       BEGIN_MAP_LOOP(cctkGH, CCTK_GF) {
-        regssss.at(map) = vhh.at(map)->regions;
+        superregsss.at(map) = vhh.at(map)->superregions;
       } END_MAP_LOOP;
       
       // Check whether to recompose
       CCTK_INT const do_recompose =
-        Carpet_RegridMaps (cctkGH, & regssss, force_recompose);
+        Carpet_RegridMaps (cctkGH, & superregsss, & regssss, force_recompose);
       assert (do_recompose >= 0);
       did_change = did_change or do_recompose;
       
       if (do_recompose) {
         BEGIN_MAP_LOOP(cctkGH, CCTK_GF) {
+          gh::rregs const & superregss = superregsss.at(map);
           gh::mregs const & regsss = regssss.at(map);
-          RegridMap (cctkGH, map, regsss);
+          RegridMap (cctkGH, map, superregss, regsss);
         } END_MAP_LOOP;
       }
       
@@ -209,6 +216,7 @@ namespace Carpet {
   void
   RegridMap (cGH const * const cctkGH,
              int const m,
+             gh::rregs const & superregss,
              gh::mregs const & regsss)
   {
     DECLARE_CCTK_PARAMETERS;
@@ -221,13 +229,15 @@ namespace Carpet {
     // not change
     
     // Regrid
-    vhh.at(m)->regrid (regsss);
+    vhh.at(m)->regrid (superregss, regsss);
     
     // Write grid structure to file
     OutputGridStructure (cctkGH, m, regsss);
     OutputGridCoordinates (cctkGH, m, regsss);
     
-    if (verbose) OutputGrids (cctkGH, m, * vhh.at(m), * vdd.at(m));
+    if (verbose or veryverbose) {
+      OutputGrids (cctkGH, m, * vhh.at(m), * vdd.at(m));
+    }
     
     Waypoint ("Done regridding map %d.", m);
   }
@@ -717,18 +727,20 @@ namespace Carpet {
   // SplitRegions_AlongZ for grid arrays)
   void
   SplitRegions (cGH const * const cctkGH,
+                vector<region_t> & superregs,
                 vector<region_t> & regs)
   {
     DECLARE_CCTK_PARAMETERS;
     
+    assert (regs.empty());
     if (CCTK_EQUALS (processor_topology, "along-z")) {
-      SplitRegions_AlongZ (cctkGH, regs);
+      SplitRegions_AlongZ (cctkGH, superregs, regs);
     } else if (CCTK_EQUALS (processor_topology, "along-dir")) {
-      SplitRegions_AlongDir (cctkGH, regs, split_direction);
+      SplitRegions_AlongDir (cctkGH, superregs, regs, split_direction);
     } else if (CCTK_EQUALS (processor_topology, "automatic")) {
-      SplitRegions_Automatic (cctkGH, regs);
+      SplitRegions_Automatic (cctkGH, superregs, regs);
     } else if (CCTK_EQUALS (processor_topology, "manual")) {
-      SplitRegions_AsSpecified (cctkGH, regs);
+      SplitRegions_AsSpecified (cctkGH, superregs, regs);
     } else {
       assert (0);
     }
@@ -835,29 +847,38 @@ namespace Carpet {
   
   void
   SplitRegions_AlongZ (cGH const * const cctkGH,
+                       vector<region_t> & superregs,
                        vector<region_t> & regs)
   {
-    SplitRegions_AlongDir (cctkGH, regs, dim - 1);
+    assert (regs.empty());
+    SplitRegions_AlongDir (cctkGH, superregs, regs, dim - 1);
   }
   
   
   
   void
   SplitRegions_AlongDir (cGH const * const cctkGH,
+                         vector<region_t> & superregs,
                          vector<region_t> & regs,
                          int const dir)
   {
+    assert (regs.empty());
+    
     // Something to do?
-    if (regs.size() == 0) {
+    if (superregs.size() == 0) {
       return;
     }
     
     const int nprocs = CCTK_nProcs (cctkGH);
     
-    assert (regs.size() == 1);
+    assert (superregs.size() == 1);
+    regs = superregs;
     
     if (nprocs == 1) {
       regs.at(0).processor = 0;
+      pseudoregion_t const preg (regs.at(0).extent, regs.at(0).processor);
+      assert (superregs.at(0).processors == NULL);
+      superregs.at(0).processors = new ipfulltree (preg);
       return;
     }
     
@@ -870,6 +891,8 @@ namespace Carpet {
     const b2vect obnd0 = reg0.outer_boundaries;
     
     regs.resize (nprocs);
+    vector<int> bounds (nprocs+1);
+    vector<ipfulltree *> subtrees (nprocs);
     for (int c=0; c<nprocs; ++c) {
       ivect cstr = rstr0;
       ivect clb  = rlb0;
@@ -892,7 +915,14 @@ namespace Carpet {
       obnd[0][dir] &= clb[dir] == rlb0[dir];
       obnd[1][dir] &= cub[dir] == rub0[dir];
       proc = c;
+      pseudoregion_t const preg (reg.extent, c);
+      bounds.at(c) = reg.extent.lower()[dir];
+      subtrees.at(c) = new ipfulltree (preg);
     }
+    bounds.at(nprocs) = rub0[dir] + rstr0[dir];
+    
+    assert (superregs.at(0).processors == NULL);
+    superregs.at(0).processors = new ipfulltree (dir, bounds, subtrees);
     
     for (int c=0; c<(int)regs.size(); ++c) {
       assert (regs.at(c).processor == c);
@@ -903,10 +933,15 @@ namespace Carpet {
   
   void
   SplitRegions_Automatic (cGH const * const cctkGH,
+                          vector<region_t> & superregs,
                           vector<region_t> & regs)
   {
-    vector<vector<region_t> > regss (1, regs);
-    SplitRegionsMaps_Automatic (cctkGH, regss);
+    assert (regs.empty());
+    vector<vector<region_t> > superregss (1, superregs);
+    vector<vector<region_t> > regss (1);
+    SplitRegionsMaps_Automatic (cctkGH, superregss, regss);
+    assert (superregss.size() == 1);
+    superregs = regss.at(0);
     assert (regss.size() == 1);
     regs = regss.at(0);
   }
@@ -915,20 +950,23 @@ namespace Carpet {
   
   static void
   SplitRegions_AsSpecified (cGH const * const cctkGH,
+                            vector<region_t> & superregs,
                             vector<region_t> & regs)
   {
     DECLARE_CCTK_PARAMETERS;
     
+    assert (regs.empty());
+    
     // Something to do?
-    if (regs.size() == 0) {
+    if (superregs.size() == 0) {
       return;
     }
     
     const int nprocs = CCTK_nProcs (cctkGH);
     
-    assert (regs.size() == 1);
+    assert (superregs.size() == 1);
     
-    region_t const & reg0 = regs.at(0);
+    region_t const & reg0 = superregs.at(0);
     const ivect rstr0  = reg0.extent.stride();
     const ivect rlb0   = reg0.extent.lower();
     const ivect rub0   = reg0.extent.upper() + rstr0;
@@ -956,9 +994,17 @@ namespace Carpet {
     const ivect rem = glonp % nprocs_dir;
     const ivect step = locnp * cstr;
     assert (dim==3);
+    
+    vector<int> boundsz(nprocs_dir[2]);
+    vector<ipfulltree *> subtreesz(nprocs_dir[2]+1);
     for (int k=0; k<nprocs_dir[2]; ++k) {
+      vector<int> boundsy(nprocs_dir[2]);
+      vector<ipfulltree *> subtreesy(nprocs_dir[2]+1);
       for (int j=0; j<nprocs_dir[1]; ++j) {
+        vector<int> boundsx(nprocs_dir[2]);
+        vector<ipfulltree *> subtreesx(nprocs_dir[2]+1);
 	for (int i=0; i<nprocs_dir[0]; ++i) {
+          
 	  const int c = i + nprocs_dir[0] * (j + nprocs_dir[1] * k);
 	  const ivect ipos (i, j, k);
 	  ivect clb = rlb0 + step *  ipos;
@@ -988,9 +1034,20 @@ namespace Carpet {
           obnd[0] &= clb == rlb0;
           obnd[1] &= cub == rub0;
           proc = c;
+          
+          pseudoregion_t preg (reg.extent, c);
+          subtreesx.at(i) = new ipfulltree (preg);
 	}
+        boundsx.at(nprocs_dir[0]) = rub0[0] + rstr0[0];
+        subtreesy.at(j) = new ipfulltree (0, boundsx, subtreesx);
       }
+      boundsy.at(nprocs_dir[1]) = rub0[1] + rstr0[1];
+      subtreesz.at(k) = new ipfulltree (1, boundsy, subtreesy);
     }
+    boundsz.at(nprocs_dir[2]) = rub0[2] + rstr0[2];
+    
+    assert (superregs.at(0).processors == NULL);
+    superregs.at(0).processors = new ipfulltree (2, boundsz, subtreesz);
     
     for (int c=0; c<(int)regs.size(); ++c) {
       assert (regs.at(c).processor == c);
@@ -1003,21 +1060,26 @@ namespace Carpet {
   // SplitRegions_AlongZ for grid arrays)
   void
   SplitRegionsMaps (cGH const * const cctkGH,
+                    vector<vector<region_t> > & superregss,
                     vector<vector<region_t> > & regss)
   {
     DECLARE_CCTK_PARAMETERS;
     
+    assert (regss.size() == superregss.size());
+    for (size_t m=0; m<regss.size(); ++m) {
+      assert (regss.AT(m).empty());
+    }
     if (CCTK_EQUALS (processor_topology, "along-z")) {
       assert (0);
-//       SplitRegionsMaps_AlongZ (cctkGH, regss);
+//       SplitRegionsMaps_AlongZ (cctkGH, superregss, regss);
     } else if (CCTK_EQUALS (processor_topology, "along-dir")) {
       assert (0);
-//       SplitRegionsMaps_AlongDir (cctkGH, regss, split_direction);
+//       SplitRegionsMaps_AlongDir (cctkGH, superregss, regss, split_direction);
     } else if (CCTK_EQUALS (processor_topology, "automatic")) {
-      SplitRegionsMaps_Automatic (cctkGH, regss);
+      SplitRegionsMaps_Automatic (cctkGH, superregss, regss);
     } else if (CCTK_EQUALS (processor_topology, "manual")) {
       assert (0);
-//       SplitRegionsMaps_AsSpecified (cctkGH, regss);
+//       SplitRegionsMaps_AsSpecified (cctkGH, superregss, regss);
     } else {
       assert (0);
     }
@@ -1026,9 +1088,10 @@ namespace Carpet {
   
   
   static void
-  SplitRegionsMaps_Automatic_Recursively (bvect    const & dims,
-                                          int      const   nprocs,
-                                          region_t const & reg,
+  SplitRegionsMaps_Automatic_Recursively (bvect    const   & dims,
+                                          int      const     firstproc,
+                                          int      const     nprocs,
+                                          region_t         & superreg,
                                           vector<region_t> & newregs)
   {
     if (DEBUG) cout << "SRMAR enter" << endl;
@@ -1046,7 +1109,11 @@ namespace Carpet {
       assert (nprocs == 1);
       
       // Return argument
-      newregs.push_back (reg);
+      region_t newreg = superreg;
+      newreg.processor = firstproc;
+      newregs.push_back (newreg);
+      pseudoregion_t const preg (newreg.extent, newreg.processor);
+      superreg.processors = new ipfulltree (preg);
       
       // Check postcondition
       assert (newregs.size() == oldsize + nprocs);
@@ -1056,19 +1123,21 @@ namespace Carpet {
     }
     
     // Special case
-    if (reg.extent.empty()) {
+    if (superreg.extent.empty()) {
       if (DEBUG) cout << "SRMAR empty" << endl;
       
       // Create a new region
-      region_t newreg (reg);
-      newreg.outer_boundaries      = b2vect(false);
+      region_t newreg (superreg);
+      newreg.outer_boundaries = b2vect(false);
       if (DEBUG) cout << "SRMAR newreg " << newreg << endl;
       
       // Store
       for (int p=0; p<nprocs; ++p) {
-        newreg.processor = reg.processor + p;
+        newreg.processor = firstproc + p;
         newregs.push_back (newreg);
       }
+      // Do not set any superreg.processors information
+      // TODO: Maybe create an empty tree for this?
       
       // Check postcondition
       assert (newregs.size() == oldsize + nprocs);
@@ -1078,10 +1147,10 @@ namespace Carpet {
     }
     
     // Calculate cost factors
-    rvect rcost = cost (reg);
+    rvect rcost = cost (superreg);
     CCTK_REAL const rfact = pow (nprocs / prod(rcost), CCTK_REAL(1)/dim);
     rcost *= rfact;
-    assert (abs (prod (rcost) - nprocs) < 1.0e-6);
+    assert (abs (prod (rcost) - nprocs) <= 1.0e-6);
     if (DEBUG) cout << "SRMA shapes " << rcost << endl;
     
     // Choose a direction
@@ -1134,7 +1203,8 @@ namespace Carpet {
     
     // Split the region
     vector<int> mynpoints(nslices);
-    int const npoints = (reg.extent.shape() / reg.extent.stride())[mydim];
+    int const npoints =
+      (superreg.extent.shape() / superreg.extent.stride())[mydim];
     
     // Keep track of how many points and processors we have left to
     // distribute
@@ -1155,16 +1225,18 @@ namespace Carpet {
     
     // Create the regions and recurse
     if (DEBUG) cout << "SRMAR " << mydim << ": create bboxes" << endl;
-    ivect lo = reg.extent.lower();
-    ivect up = reg.extent.upper();
-    ivect const str = reg.extent.stride();
-    int p = reg.processor;
+    ivect lo = superreg.extent.lower();
+    ivect up = superreg.extent.upper();
+    ivect const str = superreg.extent.stride();
+    int p = firstproc;
+    vector<int> bounds (nslices+1);
+    vector<ipfulltree *> subtrees (nslices);
     for (int n=0; n<nslices; ++n) {
       if (DEBUG) cout << "SRMAR " << mydim << " n " << n << endl;
       
       // Create a new region
       up[mydim] = lo[mydim] + (mynpoints.at(n) - 1) * str[mydim];
-      b2vect newob (reg.outer_boundaries);
+      b2vect newob (superreg.outer_boundaries);
       if (n > 0) {
         newob[0][mydim] = false;
       }
@@ -1172,23 +1244,33 @@ namespace Carpet {
         up[mydim] = lo[mydim] + (mynpoints.at(n) - 1) * str[mydim];
         newob[1][mydim] = false;
       }
-      region_t newreg (reg);
+      region_t newreg (superreg);
       newreg.extent = ibbox(lo, up, str);
-      newreg.outer_boundaries      = newob;
-      newreg.processor = p;
+      newreg.outer_boundaries = newob;
       if (DEBUG) cout << "SRMAR " << mydim << " newreg " << newreg << endl;
       
       // Recurse
+      bounds.at(n) = lo[mydim];
       SplitRegionsMaps_Automatic_Recursively
-        (newdims, mynprocs.at(n), newreg, newregs);
+        (newdims, p, mynprocs.at(n), newreg, newregs);
       if (DEBUG) cout << "SRMAR " << mydim << " newregs " << newregs << endl;
       
+      assert (newreg.processors != NULL);
+      subtrees.at(n) = newreg.processors;
+      newreg.processors = NULL;
+      newreg.processor = p;
+      
       // Next slice
       lo[mydim] = up[mydim] + str[mydim];
       p += mynprocs.at(n);
     }
-    assert (up[mydim] == reg.extent.upper()[mydim]);
-    assert (p == reg.processor + nprocs);
+    assert (up[mydim] == superreg.extent.upper()[mydim]);
+    assert (p == firstproc + nprocs);
+    bounds.at(nslices) = up[mydim] + str[mydim];
+    
+    // Create tree
+    assert (superreg.processors == NULL);
+    superreg.processors = new ipfulltree (mydim, bounds, subtrees);
     
     // Check postcondition
     assert (newregs.size() == oldsize + nprocs);
@@ -1200,14 +1282,21 @@ namespace Carpet {
   
   void
   SplitRegionsMaps_Automatic (cGH const * const cctkGH,
+                              vector<vector<region_t> > & superregss,
                               vector<vector<region_t> > & regss)
   {
     if (DEBUG) cout << "SRMA enter" << endl;
     
-    int const nmaps = regss.size();
+    int const nmaps = superregss.size();
+    int minmap = 1000000000;
+    for (int m=0; m<nmaps; ++m) {
+      for (int r=0; r<int(superregss.at(m).size()); ++r) {
+        minmap = min (minmap, superregss.at(m).at(r).map);
+      }
+    }
     int nregs = 0;
     for (int m=0; m<nmaps; ++m) {
-      nregs += regss.at(m).size();
+      nregs += superregss.at(m).size();
     }
     if (DEBUG) cout << "SRMA nregs " << nregs << endl;
     
@@ -1217,84 +1306,12 @@ namespace Carpet {
     }
     
     // Collect slices
-    vector<region_t> regs;
+    vector<region_t> superregs;
     {
-#if 0
-      regs.resize (nregs);
-      int r=0;
       for (int m=0; m<nmaps; ++m) {
-        for (int c=0; c<int(regss.at(m).size()); ++c, ++r) {
-          regs.at(r) = regss.at(m).at(c);
-          regs.at(r).map       = m;
-          regs.at(r).processor = -1;
-        }
+        combine_regions (superregss.at(m), superregs);
       }
-      assert (r == nregs);
-#endif
-      for (int m=0; m<nmaps; ++m) {
-        ibset comps;
-        ibset cobnds[2][dim];
-        for (size_t c=0; c<regss.at(m).size(); ++c) {
-          region_t const & reg = regss.at(m).at(c);
-          comps += reg.extent;
-          for (int f = 0; f < 2; ++ f) {
-            for (int d = 0; d < dim; ++ d) {
-              if (reg.outer_boundaries[f][d]) {
-                ibbox bnd = reg.extent;
-                ivect lo = bnd.lower();
-                ivect up = bnd.upper();
-                if (f==0) {
-                  up[d] = lo[d];
-                } else {
-                  lo[d] = up[d];
-                }
-                bnd = ibbox (lo, up, bnd.stride());
-                cobnds[f][d] += bnd;
-              }
-            }
-          }
-        }
-        comps.normalize();
-        for (int f = 0; f < 2; ++ f) {
-          for (int d = 0; d < dim; ++ d) {
-            cobnds[f][d].normalize();
-          }
-        }
-        size_t const needsize = regs.size() + comps.setsize();
-        if (regs.capacity() < needsize) {
-          regs.reserve (1000 + 2 * needsize);
-        }
-        for (ibset::const_iterator ci = comps.begin(); ci != comps.end(); ++ ci)
-        {
-          ibbox const & c = * ci;
-          b2vect obnds;
-          for (int f = 0; f < 2; ++ f) {
-            for (int d = 0; d < dim; ++ d) {
-              obnds[f][d] = cobnds[f][d].intersects (c);
-              if (obnds[f][d]) {
-                ivect lo = c.lower();
-                ivect up = c.upper();
-                if (f) lo[d]=up[d]; else up[d]=lo[d];
-                ibbox const cbnds (lo, up, c.stride());
-                if (not ((cobnds[f][d] & ibset(c)) == ibset(cbnds))) {
-                  cout << "cobnds[f][d] = " << cobnds[f][d] << endl
-                       << "ibset(c) = " << ibset(c) << endl
-                       << "(cobnds[f][d] & ibset(c)) = " << (cobnds[f][d] & ibset(c)) << endl
-                       << "ibset(cbnds) = " << ibset(cbnds) << endl;
-                }
-                assert ((cobnds[f][d] & ibset(c)) == ibset(cbnds));
-              }
-            }
-          }
-          region_t reg;
-          reg.extent           = c;
-          reg.outer_boundaries = obnds;
-          reg.map              = m;
-          reg.processor        = -1;
-          regs.push_back (reg);
-        }
-      } // for m
-      nregs = regs.size();
+      nregs = superregs.size();
       
       // If the last region was removed, add a new empty region again.
       // A set of regions (corresponding to a refinement level or a
@@ -1307,8 +1324,8 @@ namespace Carpet {
         reg.outer_boundaries = b2vect (bvect (true), bvect (true));
         reg.map              = 0;
         reg.processor        = -1;
-        regs.push_back (reg);
-        nregs = regs.size();
+        superregs.push_back (reg);
+        nregs = superregs.size();
       }
     }
     
@@ -1326,7 +1343,7 @@ namespace Carpet {
     if (DEBUG) cout << "SRMA: distributing processors to regions" << endl;
     vector<CCTK_REAL> mycosts(nregs);
     for (int r=0; r<nregs; ++r) {
-      mycosts.at(r) = prod (cost (regs.at(r)));
+      mycosts.at(r) = prod (cost (superregs.at(r)));
     }
     int nregs_left = newnregs;
     vector<int> mynprocs(nregs);
@@ -1361,44 +1378,77 @@ namespace Carpet {
     vector<region_t> newregs;
     newregs.reserve (newnregs);
     for (int r=0, p=0; r<nregs; p+=mynprocs.at(r), ++r) {
-      regs.at(r).processor = p;
-      if (DEBUG) cout << "SRMA reg[" << r << "] " << regs.at(r) << endl;
+      if (DEBUG) cout << "SRMA superreg[" << r << "] " << superregs.at(r) << endl;
       bvect const dims = false;
       SplitRegionsMaps_Automatic_Recursively
-        (dims, mynprocs.at(r), regs.at(r), newregs);
+        (dims, p, mynprocs.at(r), superregs.at(r), newregs);
     } // for r
     if (DEBUG) cout << "SRMA newregs " << newregs << endl;
     assert (int(newregs.size()) == newnregs);
     
+    // Count components per map
+    vector<int> myncomps(nmaps, 0);
+    for (int r=0; r<newnregs; ++r) {
+      int const m = newregs.at(r).map - minmap;
+      assert (m>=0 and m<nmaps);
+      ++ myncomps.at(m);
+    }
+    vector<int> mynregs(nmaps, 0);
+    for (int r=0; r<nregs; ++r) {
+      int const m = superregs.at(r).map - minmap;
+      assert (m>=0 and m<nmaps);
+      ++ mynregs.at(m);
+    }
+    
     // Convert virtual to real processors
     for (int r=0; r<newnregs; ++r) {
       newregs.at(r).processor /= ncomps;
       assert (newregs.at(r).processor >= 0 and
               newregs.at(r).processor < nprocs);
     }
+    {
+      vector<int> tmpncomps(nmaps, 0);
+      for (int r=0; r<nregs; ++r) {
+        int const m = superregs.at(r).map - minmap;
+        ipfulltree * const regf = superregs.at(r).processors;
+        assert (regf != NULL);
+        for (ipfulltree::iterator
+               fti = regf->begin(); fti != regf->end(); ++ fti)
+        {
+          pseudoregion_t & preg = (* fti).payload();
+          preg.component = tmpncomps.at(m)++;
+          // preg.processor /= ncomps;
+        }
+      }
+      for (int m=0; m<nmaps; ++m) {
+        assert (tmpncomps.at(m) == myncomps.at(m));
+      }
+    }
     
     // Distribute regions
-    // Count components per map
-    vector<int> myncomps(nmaps);
-    for (int r=0; r<newnregs; ++r) {
-      int const m = newregs.at(r).map;
-      assert (m>=0 and m<nmaps);
-      ++ myncomps.at(m);
-    }
     // Allocate regions
+    assert ((int)regss.size() == nmaps);
     for (int m=0; m<nmaps; ++m) {
-      regss.at(m).resize (0);
+      assert (regss.at(m).empty());
       regss.at(m).reserve (myncomps.at(m));
+      superregss.at(m).clear();
+      superregss.at(m).reserve (mynregs.at(m));
     }
     // Assign regions
     for (int r=0; r<newnregs; ++r) {
-      int const m = newregs.at(r).map;
+      int const m = newregs.at(r).map - minmap;
       assert (m>=0 and m<nmaps);
       regss.at(m).push_back (newregs.at(r));
     }
+    for (int r=0; r<nregs; ++r) {
+      int const m = superregs.at(r).map - minmap;
+      assert (m>=0 and m<nmaps);
+      superregss.at(m).push_back (superregs.at(r));
+    }
     // Check sizes
     for (int m=0; m<nmaps; ++m) {
       assert (int(regss.at(m).size()) == myncomps.at(m));
+      assert (int(superregss.at(m).size()) == mynregs.at(m));
     }
     
     if (DEBUG) cout << "SRMA exit" << endl;
@@ -1477,7 +1527,7 @@ namespace Carpet {
   void
   MakeMultigridBoxes (cGH const * const cctkGH,
                       int const m,
-                      vector<vector<region_t> >  const & regss,
+                      vector<vector<region_t> > const & regss,
                       vector<vector<vector<region_t> > > & regsss)
   {
     regsss.resize (mglevels);
diff --git a/Carpet/Carpet/src/SetupGH.cc b/Carpet/Carpet/src/SetupGH.cc
index 5e34fcd30..b2bf94b70 100644
--- a/Carpet/Carpet/src/SetupGH.cc
+++ b/Carpet/Carpet/src/SetupGH.cc
@@ -128,7 +128,7 @@ namespace Carpet {
   static void
   find_processor_decomposition
   (cGH const * cctkGH,
-   vector<vector<vector<region_t> > >          & regsss,
+   vector<vector<vector<region_t> > >          & superregsss,
    vector<vector<vector<vector<region_t> > > > & regssss);
   
   static void
@@ -204,7 +204,7 @@ namespace Carpet {
       // Do not call Util_GetHostName.  Certain InfiniBand libraries
       // do not allow calling fork or exec, and getting the host name
       // seems also not allowed.  It leads to random crashes.
-      if (verbose) {
+      if (verbose or veryverbose) {
         // Collect host names
         char hostnamebuf[1000];
         Util_GetHostName (hostnamebuf, sizeof hostnamebuf);
@@ -572,13 +572,13 @@ namespace Carpet {
   {
     DECLARE_CCTK_PARAMETERS;
     
-    vector<vector<vector<region_t> > > regsss (maps);
+    vector<vector<vector<region_t> > > superregsss (maps);
     for (int m=0; m<maps; ++m) {
-      set_base_extent (m, regsss.at(m));
+      set_base_extent (m, superregsss.at(m));
     }
     
     vector<vector<vector<vector<region_t> > > > regssss;
-    find_processor_decomposition (cctkGH, regsss, regssss);
+    find_processor_decomposition (cctkGH, superregsss, regssss);
     
     for (int m=0; m<maps; ++m) {
       
@@ -586,7 +586,7 @@ namespace Carpet {
       CheckRegions (regssss.at(m));
       
       // Recompose grid hierarchy
-      vhh.at(m)->regrid (regssss.at(m));
+      vhh.at(m)->regrid (superregsss.at(m), regssss.at(m));
       int const rl = 0;
       vhh.at(m)->recompose (rl, false);
       
@@ -827,21 +827,25 @@ namespace Carpet {
                               ivect const & convoffsets)
   {
     // Set refinement structure for scalars and arrays
-    vector<region_t> regs(1);
-    int const m=0;
-    regs.at(0).extent = arrdata.at(group).at(m).hh->baseextents.at(0).at(0);
-    regs.at(0).outer_boundaries = b2vect(true);
-    regs.at(m).map = m;
+    vector<region_t> superregs(1);
+    int const m = 0;
+    int const rl = 0;
+    int const c = 0;
+    superregs.at(c).extent =
+      arrdata.at(group).at(m).hh->baseextents.at(rl).at(c);
+    superregs.at(c).outer_boundaries = b2vect(true);
+    superregs.at(c).map = m;
+    vector<region_t> regs;
     
     // Split it into components, one for each processor
     switch (gdata.disttype) {
     case CCTK_DISTRIB_DEFAULT: {
-      SplitRegions_Automatic (cctkGH, regs);
+      SplitRegions_Automatic (cctkGH, superregs, regs);
       break;
     }
     case CCTK_DISTRIB_CONSTANT: {
       int const d = gdata.dim==0 ? 0 : gdata.dim-1;
-      SplitRegions_AlongDir (cctkGH, regs, d);
+      SplitRegions_AlongDir (cctkGH, superregs, regs, d);
       break;
     }
     default:
@@ -849,8 +853,10 @@ namespace Carpet {
     }
     
     // Only one refinement level
+    vector<vector<region_t> > superregss(1);
+    superregss.at(rl) = superregs;
     vector<vector<region_t> > regss(1);
-    regss.at(0) = regs;
+    regss.at(rl) = regs;
     
     // Create all multigrid levels
     vector<vector<vector<region_t> > > regsss (mglevels);
@@ -895,7 +901,7 @@ namespace Carpet {
       char * const groupname = CCTK_GroupName (group);
       assert (groupname);
       Checkpoint ("Recomposing grid array group \"%s\"...", groupname);
-      arrdata.at(group).at(0).hh->regrid (regsss);
+      arrdata.at(group).at(0).hh->regrid (superregss, regsss);
       arrdata.at(group).at(0).hh->recompose (0, false);
       Checkpoint ("Done recomposing grid array group \"%s\".", groupname);
       free (groupname);
@@ -1445,7 +1451,7 @@ namespace Carpet {
   void
   find_processor_decomposition
   (cGH const * const cctkGH,
-   vector<vector<vector<region_t> > >          & regsss,
+   vector<vector<vector<region_t> > >          & superregsss,
    vector<vector<vector<vector<region_t> > > > & regssss)
   {
     DECLARE_CCTK_PARAMETERS;
@@ -1459,11 +1465,13 @@ namespace Carpet {
       for (int m=0; m<maps; ++m) {
         int const rl=0;
         
+        vector<vector<region_t> > regss(1);
+        
         // Distribute onto the processors
-        SplitRegions (cctkGH, regsss.at(m).at(rl));
+        SplitRegions (cctkGH, superregsss.at(m).at(rl), regss.at(rl));
         
         // Create all multigrid levels
-        MakeMultigridBoxes (cctkGH, m, regsss.at(m), regssss.at(m));
+        MakeMultigridBoxes (cctkGH, m, regss, regssss.at(m));
       } // for m
       
     } else {
@@ -1471,14 +1479,18 @@ namespace Carpet {
       
       int const rl=0;
       
-      vector<vector<region_t> >  regss(maps);
+      vector<vector<region_t> > superregss(maps);
       for (int m=0; m<maps; ++m) {
-        regss.at(m) = regsss.at(m).at(rl);
+        superregss.at(m) = superregsss.at(m).at(rl);
       }
       
-      SplitRegionsMaps (cctkGH, regss);
+      vector<vector<region_t> > regss(maps);
+      SplitRegionsMaps (cctkGH, superregss, regss);
       
+      vector<vector<vector<region_t> > > regsss(maps);
       for (int m=0; m<maps; ++m) {
+        superregsss.at(m).at(rl) = superregss.at(m);
+        regsss.at(m).resize(1);
         regsss.at(m).at(rl) = regss.at(m);
       }
       
diff --git a/Carpet/Carpet/src/functions.hh b/Carpet/Carpet/src/functions.hh
index bbe2ff724..ef86ca3c9 100644
--- a/Carpet/Carpet/src/functions.hh
+++ b/Carpet/Carpet/src/functions.hh
@@ -22,6 +22,7 @@
 namespace Carpet {
   
   using namespace std;
+  using namespace CarpetLib;
   
   int SyncGroupsByDirI (const cGH* cctkGH, int num_groups,
                         const int* groups, const int* directions);
@@ -102,6 +103,7 @@ namespace Carpet {
   void
   RegridMap (cGH const * cctkGH,
              int m,
+             gh::rregs const & supeerregss,
              gh::mregs const & regsss);
   void
   PostRegrid (cGH const * cctkGH);
@@ -137,28 +139,30 @@ namespace Carpet {
   // Functions for recomposing the grid hierarchy
   void
   SplitRegions (cGH const * cctkGH,
+                vector<region_t> & superregs,
                 vector<region_t> & regs);
   void
   SplitRegions_AlongZ (cGH const * cctkGH,
+                       vector<region_t> & superregs,
                        vector<region_t> & regs);
   void
   SplitRegions_AlongDir (cGH const * cctkGH,
+                         vector<region_t> & superregs,
                          vector<region_t> & regs,
-                         int const dir);
+                         int dir);
   void
   SplitRegions_Automatic (cGH const * cctkGH,
+                          vector<region_t> & superregs,
                           vector<region_t> & regs);
-  void
-  SplitRegions (cGH const * cctkGH,
-                vector<region_t> & regss);
-  
-  void
-  SplitRegionsMaps_Automatic (cGH const * cctkGH,
-                              vector<vector<region_t> > & regss);
   
   void
   SplitRegionsMaps (cGH const * cctkGH,
+                    vector<vector<region_t> > & superregss,
                     vector<vector<region_t> > & regss);
+  void
+  SplitRegionsMaps_Automatic (cGH const * cctkGH,
+                              vector<vector<region_t> > & superregss,
+                              vector<vector<region_t> > & regss);
   
   void
   MakeMultigridBoxes (cGH const * cctkGH,
diff --git a/Carpet/CarpetIOHDF5/src/Input.cc b/Carpet/CarpetIOHDF5/src/Input.cc
index 5242daedc..7abdd1232 100644
--- a/Carpet/CarpetIOHDF5/src/Input.cc
+++ b/Carpet/CarpetIOHDF5/src/Input.cc
@@ -108,7 +108,7 @@ void CarpetIOHDF5_RecoverGridStructure (CCTK_ARGUMENTS)
   // checkpoint file, or if the number of maps is wrong
   assert (int(fileset.grid_structure.size()) == maps);
   
-  vector<vector<vector<region_t> > > regsss = fileset.grid_structure;
+  vector<vector<vector<region_t> > > superregsss = fileset.grid_structure;
   vector<vector<vector<vector<region_t> > > > regssss (maps);
   
   int type;
@@ -122,26 +122,28 @@ void CarpetIOHDF5_RecoverGridStructure (CCTK_ARGUMENTS)
     // Distribute each map independently
     
     for (int m = 0; m < maps; ++ m) {
-      
-      vector <vector <region_t> > & regss = regsss.at(m);
+      vector <vector <region_t> > & superregss = superregsss.at(m);
       
       // Make multiprocessor aware
-      for (size_t rl = 0; rl < regss.size(); ++ rl) {
-        Carpet::SplitRegions (cctkGH, regss.at(rl));
+      vector <vector <region_t> > regss(superregss.size());
+      for (size_t rl = 0; rl < superregss.size(); ++ rl) {
+        Carpet::SplitRegions (cctkGH, superregss.at(rl), regss.at(rl));
       } // for rl
       
       // Make multigrid aware
-      Carpet::MakeMultigridBoxes (cctkGH, Carpet::map, regss, regssss.at(m));
+      Carpet::MakeMultigridBoxes (cctkGH, m, regss, regssss.at(m));
       
     } // for m
     
-  } else {
+  } else {                     // if regrid_in_level_mode
     // Distribute all maps at the same time
     
+    vector<vector<vector<region_t> > > regsss(maps);
+    
     // Count levels
     vector <int> rls (maps);
     for (int m = 0; m < maps; ++ m) {
-      rls.at(m) = regsss.at(m).size();
+      rls.at(m) = superregsss.at(m).size();
     }
     int maxrl = 0;
     for (int m = 0; m < maps; ++ m) {
@@ -149,30 +151,33 @@ void CarpetIOHDF5_RecoverGridStructure (CCTK_ARGUMENTS)
     }
     // All maps must have the same number of levels
     for (int m = 0; m < maps; ++ m) {
+      superregsss.at(m).resize (maxrl);
       regsss.at(m).resize (maxrl);
     }
     
     // Make multiprocessor aware
     for (int rl = 0; rl < maxrl; ++ rl) {
-      vector <vector <region_t> > regss (maps);
+      vector <vector <region_t> > superregss (maps);
       for (int m = 0; m < maps; ++ m) {
-        regss.at(m) = regsss.at(m).at(rl);
+        superregss.at(m) = superregsss.at(m).at(rl);
       }
-      Carpet::SplitRegionsMaps (cctkGH, regss);
+      vector <vector <region_t> > regss (maps);
+      Carpet::SplitRegionsMaps (cctkGH, superregss, regss);
       for (int m = 0; m < maps; ++ m) {
+        superregsss.at(m).at(rl) = superregss.at(m);
         regsss.at(m).at(rl) = regss.at(m);
       }
     } // for rl
     
-      // Make multigrid aware
+    // Make multigrid aware
     Carpet::MakeMultigridBoxesMaps (cctkGH, regsss, regssss);
     
-  } // if
+  } // if regrid_in_level_mode
   
   for (int m = 0; m < maps; ++ m) {
     
     // Regrid
-    RegridMap (cctkGH, m, regssss.at(m));
+    RegridMap (cctkGH, m, superregsss.at(m), regssss.at(m));
     
     // Set time hierarchy correctly after RegridMap created it
     for (int ml = 0; ml < mglevels; ++ ml) {
diff --git a/Carpet/CarpetInterp/param.ccl b/Carpet/CarpetInterp/param.ccl
index 3ba04386d..d75bad138 100644
--- a/Carpet/CarpetInterp/param.ccl
+++ b/Carpet/CarpetInterp/param.ccl
@@ -14,6 +14,14 @@ CCTK_REAL ipoison "Integer poison value" STEERABLE=always
   *:* :: ""
 } -420042
 
+BOOLEAN tree_search "Use a tree search to find the source processor" STEERABLE=always
+{
+} "no"
+
+BOOLEAN check_tree_search "Cross-check the result of the tree search" STEERABLE=always
+{
+} "yes"
+
 SHARES: Cactus
 
 USES CCTK_REAL cctk_initial_time
diff --git a/Carpet/CarpetInterp/src/interp.cc b/Carpet/CarpetInterp/src/interp.cc
index 9ebb438f2..3fcadc9ac 100644
--- a/Carpet/CarpetInterp/src/interp.cc
+++ b/Carpet/CarpetInterp/src/interp.cc
@@ -65,7 +65,7 @@ namespace CarpetInterp {
                                    int const N_input_arrays,
                                    int const N_output_arrays,
                                    bool& have_source_map,
-                                   int & num_time_derivs,
+                                   vector<int> & num_time_derivs,
                                    int & prolongation_order_time,
                                    CCTK_REAL & current_time,
                                    CCTK_REAL & delta_time,
@@ -109,7 +109,7 @@ namespace CarpetInterp {
                           CCTK_INT* const per_proc_retvals,
                           vector<CCTK_INT> & operand_indices,
                           vector<CCTK_INT> & time_deriv_order,
-                          int const num_time_derivs,
+                          vector<int> const & num_time_derivs,
                           CCTK_INT const local_interp_handle,
                           CCTK_INT const param_table_handle,
                           CCTK_REAL const current_time,
@@ -134,8 +134,8 @@ namespace CarpetInterp {
                                 vector<CCTK_INT> & interp_num_time_levels,
                                 CCTK_INT const local_interp_handle,
                                 CCTK_INT const param_table_handle,
-                                int const num_tl,
-                                bool const need_time_interp,
+                                vector<int> const & num_tl,
+                                vector<bool> const & need_time_interp,
                                 CCTK_REAL const current_time,
                                 CCTK_REAL const delta_time,
                                 int const N_input_arrays,
@@ -353,7 +353,7 @@ namespace CarpetInterp {
     vector<CCTK_INT> operand_indices (N_output_arrays);
     vector<CCTK_INT> time_deriv_order (N_output_arrays);
     bool have_source_map;
-    int num_time_derivs;
+    vector<int> num_time_derivs;
     CCTK_REAL current_time, delta_time;
     int prolongation_order_time;
 
@@ -398,8 +398,8 @@ namespace CarpetInterp {
     // that this processor is to receive from others
     vector<int> N_points_from (dist::size());
     {
-      assert (N_points_from.size() == dist::size());
-      assert (N_points_to.size() == dist::size());
+      assert ((int)N_points_from.size() == dist::size());
+      assert ((int)N_points_to.size() == dist::size());
     }
     MPI_Alltoall (&N_points_to[0],   1, dist::datatype (N_points_to[0]),
                   &N_points_from[0], 1, dist::datatype (N_points_from[0]),
@@ -474,10 +474,10 @@ namespace CarpetInterp {
       vector<CCTK_REAL> tmp (N_dims * N_points_local);
       MPI_Datatype const datatype = dist::datatype (tmp[0]);
       {
-        assert (sendcnt.size() == dist::size());
-        assert (recvcnt.size() == dist::size());
-        assert (senddispl.size() == dist::size());
-        assert (recvdispl.size() == dist::size());
+        assert ((int)sendcnt.size() == dist::size());
+        assert ((int)recvcnt.size() == dist::size());
+        assert ((int)senddispl.size() == dist::size());
+        assert ((int)recvdispl.size() == dist::size());
         int const sendbufsize = (int)coords_buffer.size();
         int const recvbufsize = (int)tmp.size();
         for (int n=0; n<dist::size(); ++n) {
@@ -491,7 +491,7 @@ namespace CarpetInterp {
       }
 #ifndef _NDEBUG
 #pragma omp parallel for
-      for (int i=0; i<tmp.size(); ++i) {
+      for (int i=0; i<(int)tmp.size(); ++i) {
         tmp.AT(i) = poison;
       }
 #endif
@@ -501,7 +501,7 @@ namespace CarpetInterp {
 #ifndef _NDEBUG
       {
         vector<bool> filled(tmp.size(), false);
-        for (size_t n=0; n<dist::size(); ++n) {
+        for (int n=0; n<(int)dist::size(); ++n) {
           //#pragma omp parallel for
           for (int i=0; i<recvcnt.AT(n); ++i) {
             assert (not filled.AT(recvdispl.AT(n)+i));
@@ -509,24 +509,24 @@ namespace CarpetInterp {
           }
         }
         bool error = false;
-        for (int i=0; i<filled.size(); ++i) {
+        for (int i=0; i<(int)filled.size(); ++i) {
           error = error or not (filled.AT(i));
         }
         if (error) {
-          cerr << "error" << endl;
-          cerr << "recvdispl: " << recvdispl << endl;
-          cerr << "recvcnt: " << recvcnt << endl;
-          cerr << "filled: " << filled << endl;
+          cout << "error" << endl;
+          cout << "recvdispl: " << recvdispl << endl;
+          cout << "recvcnt: " << recvcnt << endl;
+          cout << "filled: " << filled << endl;
         }
 #pragma omp parallel for
-        for (int i=0; i<filled.size(); ++i) {
+        for (int i=0; i<(int)filled.size(); ++i) {
           assert (filled.AT(i));
         }
       }
 #endif
 #ifndef _NDEBUG
 #pragma omp parallel for
-      for (int i=0; i<tmp.size(); ++i) {
+      for (int i=0; i<(int)tmp.size(); ++i) {
         assert (tmp.AT(i) != poison);
       }
 #endif
@@ -585,10 +585,10 @@ namespace CarpetInterp {
       source_map.resize (N_points_local);
       MPI_Datatype const datatype = dist::datatype (tmp[0]);
       {
-        assert (sendcnt.size() == dist::size());
-        assert (recvcnt.size() == dist::size());
-        assert (senddispl.size() == dist::size());
-        assert (recvdispl.size() == dist::size());
+        assert ((int)sendcnt.size() == dist::size());
+        assert ((int)recvcnt.size() == dist::size());
+        assert ((int)senddispl.size() == dist::size());
+        assert ((int)recvdispl.size() == dist::size());
         int const sendbufsize = (int)tmp.size();
         int const recvbufsize = (int)source_map.size();
         for (int n=0; n<dist::size(); ++n) {
@@ -602,7 +602,7 @@ namespace CarpetInterp {
       }
 #ifndef _NDEBUG
 #pragma omp parallel for
-      for (int i=0; i<source_map.size(); ++i) {
+      for (int i=0; i<(int)source_map.size(); ++i) {
         source_map[i] = ipoison;
       }
 #endif
@@ -612,7 +612,7 @@ namespace CarpetInterp {
 #ifndef _NDEBUG
       {
         vector<bool> filled(source_map.size(), false);
-        for (size_t n=0; n<dist::size(); ++n) {
+        for (int n=0; n<(int)dist::size(); ++n) {
           //#pragma omp parallel for
           for (int i=0; i<recvcnt.AT(n); ++i) {
             assert (not filled.AT(recvdispl.AT(n)+i));
@@ -620,14 +620,14 @@ namespace CarpetInterp {
           }
         }
 #pragma omp parallel for
-        for (int i=0; i<filled.size(); ++i) {
+        for (int i=0; i<(int)filled.size(); ++i) {
           assert (filled.AT(i));
         }
       }
 #endif
 #ifndef _NDEBUG
 #pragma omp parallel for
-      for (int i=0; i<source_map.size(); ++i) {
+      for (int i=0; i<(int)source_map.size(); ++i) {
         assert (source_map[i] != poison);
       }
 #endif
@@ -780,10 +780,10 @@ namespace CarpetInterp {
 
       vector<char> tmp (N_interp_points * N_output_arrays * vtypesize);
       {
-        assert (sendcnt.size() == dist::size());
-        assert (recvcnt.size() == dist::size());
-        assert (senddispl.size() == dist::size());
-        assert (recvdispl.size() == dist::size());
+        assert ((int)sendcnt.size() == dist::size());
+        assert ((int)recvcnt.size() == dist::size());
+        assert ((int)senddispl.size() == dist::size());
+        assert ((int)recvdispl.size() == dist::size());
         int const sendbufsize = (int)outputs_buffer.size();
         int const recvbufsize = (int)tmp.size() / vtypesize;
         for (int n=0; n<dist::size(); ++n) {
@@ -832,7 +832,7 @@ namespace CarpetInterp {
 
     for (int d = 0; d < N_output_arrays; d++) {
       char* output_array = static_cast<char*>(output_arrays[d]);
-      size_t offset = 0;
+      int offset = 0;
       for (int c = 0, i = 0; c < (int)allhomecnts.size(); c++) {
         assert ((int) (allhomecnts.AT(c)*(d+1) + offset) <=
                 N_output_arrays*N_interp_points);
@@ -854,7 +854,7 @@ namespace CarpetInterp {
         i += allhomecnts.AT(c);
         offset += N_output_arrays * allhomecnts.AT(c);
       }
-      assert ((int) offset == N_output_arrays * N_interp_points);
+      assert (offset == N_output_arrays * N_interp_points);
     }
 
     // set this processor's overall local interpolator status code
@@ -876,7 +876,7 @@ namespace CarpetInterp {
                                    int const N_input_arrays,
                                    int const N_output_arrays,
                                    bool& have_source_map,
-                                   int& num_time_derivs,
+                                   vector<int>& num_time_derivs,
                                    int& prolongation_order_time,
                                    CCTK_REAL& current_time,
                                    CCTK_REAL& delta_time,
@@ -937,12 +937,12 @@ namespace CarpetInterp {
                                  &time_deriv_order.front(), "time_deriv_order");
     if (iret == UTIL_ERROR_TABLE_NO_SUCH_KEY) {
       time_deriv_order.assign (time_deriv_order.size(), 0);
-      num_time_derivs = 0;
+      num_time_derivs.assign (N_output_arrays, 0);
     } else {
       assert (iret == N_output_arrays);
-      num_time_derivs = 0;
+      num_time_derivs.resize (N_output_arrays);
       for (int m = 0; m < N_output_arrays; ++m) {
-        num_time_derivs = max (num_time_derivs, (int)time_deriv_order[m]);
+        num_time_derivs.AT(m) = time_deriv_order[m];
       }
     }
 
@@ -960,8 +960,139 @@ namespace CarpetInterp {
 
     return 0;
   }
+  
+  
 
+  // Find the component and integer index to which a grid point
+  // belongs.  This uses a linear search over all components, which
+  // does NOT scale with the number of components.
+  static
+  void
+  find_location_linear (gh const * restrict const hh,
+                        rvect const & restrict pos,
+                        rvect const & restrict lower,
+                        rvect const & restrict upper,
+                        rvect const & restrict delta,
+                        int const ml,
+                        int const minrl, int const maxrl,
+                        int & restrict rl,
+                        int & restrict c)
+  {
+    // cout << "CarpetInterp: assign: m=" << m << " pos=" << pos << endl;
+    
+    assert (ml>=0 and ml<mglevels);
+    assert (minrl>=0 and minrl<maxrl and maxrl<=reflevels);
+    
+    CCTK_REAL const rone  = 1.0;
+    CCTK_REAL const rhalf = rone / 2;
+    
+    if (all (pos >= lower and pos <= upper)) {
+      // The point is within the domain
+      
+      // Try finer levels first
+      for (rl = maxrl-1; rl >= minrl; --rl) {
+        
+        ivect const fact =
+          maxspacereflevelfact / spacereffacts.AT(rl) * ipow(mgfact, ml);
+        ivect const ipos =
+          ivect(floor((pos - lower) / (delta * rvect(fact)) + rhalf)) * fact;
+        
+        ivect const & stride = hh->baseextent(ml,rl).stride();
+        assert (all (ipos % stride == 0));
+        
+        gh::cregs const & regs = hh->regions.AT(ml).AT(rl);
+        
+        // Search all components linearly
+        for (c = 0; c < int(regs.size()); ++c) {
+          region_t const & reg = regs.AT(c);
+          if (reg.extent.contains(ipos)) {
+            // We found the refinement level, component, and index to
+            // which this grid point belongs
+            return;
+          }
+        }
+      }
+    }
+    
+    // The point does not belong to any component.  This should happen
+    // only rarely.
+    rl = -1;
+    c = -1;
+  }
+  
+  
 
+  // Find the component and integer index to which a grid point
+  // belongs.  This uses a tree search over the superregions in the
+  // grid struction, which should scale reasonably (O(n log n)) better
+  // with the number of componets components.
+  static
+  void
+  find_location_tree (gh const * restrict const hh,
+                      rvect const & restrict pos,
+                      rvect const & restrict lower,
+                      rvect const & restrict upper,
+                      rvect const & restrict delta,
+                      int const ml,
+                      int const minrl, int const maxrl,
+                      int & restrict rl,
+                      int & restrict c)
+  {
+    // cout << "CarpetInterp: assign: m=" << m << " pos=" << pos << endl;
+    
+    assert (ml>=0 and ml<mglevels);
+    assert (minrl>=0 and minrl<maxrl and maxrl<=reflevels);
+    
+    CCTK_REAL const rone  = 1.0;
+    CCTK_REAL const rhalf = rone / 2;
+    
+    if (all (pos >= lower and pos <= upper)) {
+      // The point is within the domain
+      
+      // Try finer levels first
+      for (rl = maxrl-1; rl >= minrl; --rl) {
+        
+        ivect const fact =
+          maxspacereflevelfact / spacereffacts.AT(rl) * ipow(mgfact, ml);
+        ivect const ipos =
+          ivect(floor((pos - lower) / (delta * rvect(fact)) + rhalf)) * fact;
+        
+        ivect const & stride = hh->baseextent(ml,rl).stride();
+        assert (all (ipos % stride == 0));
+        
+        gh::cregs const & regs = hh->superregions.AT(rl);
+        
+        // Search all superregions linearly.  Each superregion
+        // corresponds to a "refined region", and the number of
+        // superregions is thus presumably independent of the number
+        // of processors.
+        for (size_t r = 0; r < regs.size(); ++r) {
+          region_t const & reg = regs.AT(r);
+          if (reg.extent.contains(ipos)) {
+            // We found the superregion to which this grid point
+            // belongs
+            
+            // Search the superregion hierarchically
+            pseudoregion_t const * const preg = reg.processors->search(ipos);
+            assert (preg);
+            
+            // We now know the refinement level, component, and index
+            // to which this grid point belongs
+            c = preg->component;
+            return;
+          }
+        }
+      }
+    }
+    
+    // The point does not belong to any component.  This should happen
+    // only rarely.
+    rl = -1;
+    c = -1;
+  }
+  
+  
+  
   static void
   map_points (cGH const* const cctkGH,
               int const coord_system_handle,
@@ -981,6 +1112,8 @@ namespace CarpetInterp {
               vector<int>& home,
               vector<int>& homecnts)
   {
+    DECLARE_CCTK_PARAMETERS;
+
     bool const map_onto_processors = coords_list != NULL;
 
     if (not map_onto_processors) {
@@ -993,11 +1126,11 @@ namespace CarpetInterp {
     assert ((int)source_map.size() == npoints);
 
 
-    // Find out about the coordinates: origin and delta
-    // for the Carpet grid indices
+    // Find out about the coordinates: origin and delta for the Carpet
+    // grid indices
     vector<rvect> lower (maps);
-    vector<rvect> delta (maps); // spacing on finest possible grid
     vector<rvect> upper (maps);
+    vector<rvect> delta (maps); // spacing on finest possible grid
     
     int const grouptype = CCTK_GroupTypeI (coord_group);
     switch (grouptype) {
@@ -1008,7 +1141,6 @@ 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;
@@ -1028,15 +1160,14 @@ 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);
-        lower.AT(m) = coord_lower;
-        upper.AT(m) = coord_upper;
-        delta.AT(m) = ((coord_upper - coord_lower) /
-                       rvect (baseextent.upper() - baseextent.lower()));
-      }
+      assert (arrdata.AT(coord_group).size() == 1);
+      int const m = 0;
+      ibbox const & baseextent =
+        arrdata.AT(coord_group).AT(m).hh->baseextents.AT(mglevel).AT(0);
+      lower.AT(m) = coord_lower;
+      upper.AT(m) = coord_upper;
+      delta.AT(m) = ((coord_upper - coord_lower) /
+                     rvect (baseextent.upper() - baseextent.lower()));
       break;
     }
       
@@ -1044,15 +1175,6 @@ 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
 #pragma omp parallel for
     for (int n = 0; n < npoints; ++n) {
@@ -1070,48 +1192,73 @@ namespace CarpetInterp {
           pos[d] = coords[N_dims*n + d];
         }
       }
-
+      
       // 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) {
-
-          ivect const fact = maxspacereflevelfact / spacereffacts.AT(rl) *
-                             ipow(mgfact, mglevel);
-          ivect const ipos = ivect(floor((pos - lower.AT(m)) / (delta.AT(m) *
-                             rvect(fact)) + (CCTK_REAL) 0.5)) * fact;
-          assert (all (ipos % hh->baseextents.AT(ml).AT(rl).stride() == 0));
-
-          // 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;
-            }
+      
+      // Calculate rl, c, and proc
+      int rl, c;
+      if (not tree_search) {
+        find_location_linear
+          (hh, pos, lower.AT(m), upper.AT(m), delta.AT(m), ml, minrl, maxrl,
+           rl, c);
+      } else {
+        find_location_tree
+          (hh, pos, lower.AT(m), upper.AT(m), delta.AT(m), ml, minrl, maxrl,
+           rl, c);
+        if (check_tree_search) {
+          int rl2, c2;
+          find_location_linear
+            (hh, pos, lower.AT(m), upper.AT(m), delta.AT(m), ml, minrl, maxrl,
+             rl2, c2);
+          if (rl2!=rl or c2!=c) {
+#pragma omp critical
+            CCTK_VWarn (CCTK_WARN_ABORT, __LINE__, __FILE__, CCTK_THORNSTRING,
+                      "Inconsistent search result from find_location_tree for 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);
           }
         }
       }
-      // Point could not be mapped onto any processor's domain.
-      // Map the point (arbitrarily) to the first component of the
-      // coarsest grid
-      rl = minrl;
-      c = 0;
-      // Only warn once
-      if (map_onto_processors) {
+      
+      if (c == -1) {
+        // The point could not be mapped onto any component
+        
+        // Warn only once, namely when mapping points onto processors.
+        // (This routine is called twice; first to map points onto
+        // processors, then to map points onto components.)
+        if (map_onto_processors) {
 #pragma omp critical
-        CCTK_VWarn (CCTK_WARN_PICKY, __LINE__, __FILE__, CCTK_THORNSTRING,
-                    "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);
+          CCTK_VWarn (CCTK_WARN_PICKY, __LINE__, __FILE__, CCTK_THORNSTRING,
+                      "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);
+        }
+        
+        // Map the point (arbitrarily) to the first component of the
+        // coarsest grid
+        // TODO: Handle these points explicitly later on
+        rl = minrl;
+        c = 0;
+        
+        // Does this refinement level actually have a component?  (It
+        // may not, if this is a multi-patch simulation.)
+        if (not (hh->components(rl) > c)) {
+          cout << "CI: m=" << m << " pos=" << pos << endl;
+          CCTK_WARN (CCTK_WARN_ABORT, "Cannot interpolate a point from a patch which does not have any components at this refinement level.  Very likely your multi-patch grid structure is inconsistent, e.g. with refinement boundaries too close to a multi-patch boundary.");
+        } 
+        assert (hh->components(rl) > c);
+      }
+      
+      if (not (rl >= minrl and rl < maxrl) or
+          not (c >= 0 and c < hh->components(rl)))
+      {
+        cout << "CI: m=" << m << " rl=" << rl << " c=" << c << " ext=" << hh->extent(ml,rl,c) << endl;
       }
-    found:
       assert (rl >= minrl and rl < maxrl);
       assert (c >= 0 and c < hh->components(rl));
-
+      
       if (map_onto_processors) {
-        procs.AT(n) = hh->processor(rl, c);
-        int & this_N_points_to = N_points_to.AT(procs.AT(n));
+        int const proc = hh->processor(rl,c);
+        procs.AT(n) = proc;
+        int & this_N_points_to = N_points_to.AT(proc);
 #pragma omp atomic
         ++ this_N_points_to;
       }
@@ -1120,8 +1267,7 @@ namespace CarpetInterp {
       int & this_homecnts = homecnts.AT(component_idx (procs.AT(n), m, rl, c));
 #pragma omp atomic
       ++ this_homecnts;
-      //cout << "   p=" << procs.AT(n) << endl;
-
+      
     } // for n
   }
 
@@ -1143,7 +1289,7 @@ namespace CarpetInterp {
                           CCTK_INT* const per_proc_retvals,
                           vector<CCTK_INT>& operand_indices,
                           vector<CCTK_INT>& time_deriv_order,
-                          int const num_time_derivs,
+                          vector<int> const & num_time_derivs,
                           CCTK_INT const local_interp_handle,
                           CCTK_INT const param_table_handle,
                           CCTK_REAL const current_time,
@@ -1195,18 +1341,25 @@ namespace CarpetInterp {
           // Number of neccessary time levels
           CCTK_REAL const level_time =
             cctkGH->cctk_time / cctkGH->cctk_delta_time;
-          bool const need_time_interp
-            = (num_time_derivs > 0
-               or (fabs(current_time - level_time)
-                   > 1e-12 * (fabs(level_time) + fabs(current_time)
-                              + fabs(cctkGH->cctk_delta_time))));
-          assert (not (not want_global_mode
-                       and num_time_derivs == 0
-                       and need_time_interp));
-          int const num_tl
-            = (need_time_interp
-               ? max (num_time_derivs + 1, prolongation_order_time + 1)
-               : 1);
+          vector<int> num_tl (N_input_arrays, 0);
+          vector<bool> need_time_interp (N_output_arrays);
+          for (int m=0; m<N_output_arrays; ++m) {
+            need_time_interp.AT(m)
+              = (num_time_derivs.AT(m) > 0
+                 or (fabs(current_time - level_time)
+                     > 1e-12 * (fabs(level_time) + fabs(current_time)
+                                + fabs(cctkGH->cctk_delta_time))));
+            assert (not (not want_global_mode
+                         and num_time_derivs.AT(m) == 0
+                         and need_time_interp.AT(m)));
+            int const n = operand_indices.AT(m);
+            num_tl.AT(n)
+              = max (num_tl.AT(n),
+                     (need_time_interp.AT(m)
+                      ? max (num_time_derivs.AT(m) + 1,
+                             prolongation_order_time + 1)
+                      : 1));
+          }
 
           BEGIN_MAP_LOOP(cctkGH, CCTK_GF) {
             BEGIN_LOCAL_COMPONENT_LOOP(cctkGH, CCTK_GF) {
@@ -1400,8 +1553,8 @@ namespace CarpetInterp {
                                 vector<CCTK_INT>& interp_num_time_levels,
                                 CCTK_INT const local_interp_handle,
                                 CCTK_INT const param_table_handle,
-                                int const num_tl,
-                                bool const need_time_interp,
+                                vector<int> const & num_tl,
+                                vector<bool> const & need_time_interp,
                                 CCTK_REAL const current_time,
                                 CCTK_REAL const delta_time,
                                 int const N_input_arrays,
@@ -1490,9 +1643,13 @@ namespace CarpetInterp {
       tmp_coords[d] = coords + d*npoints;
     }
 
-    vector<vector<void *> > tmp_output_arrays (num_tl);
+    int max_num_tl = 0;
+    for (int m=0; m<N_input_arrays; ++m) {
+      max_num_tl = max(max_num_tl, num_tl.AT(m));
+    }
+    vector<vector<void *> > tmp_output_arrays (max_num_tl);
 
-    for (int tl=0; tl<num_tl; ++tl) {
+    for (int tl=0; tl<max_num_tl; ++tl) {
 
       for (int n=0; n<N_input_arrays; ++n) {
 
@@ -1502,13 +1659,14 @@ namespace CarpetInterp {
         if (vi == -1) {
           input_arrays[n] = NULL;
         } else {
-          int const interp_num_tl = interp_num_time_levels[n] > 0 ?
-                                    interp_num_time_levels[n] : num_tl;
+          int const interp_num_tl =
+            interp_num_time_levels.AT(n) > 0 ?
+            interp_num_time_levels.AT(n) : num_tl.AT(n);
 
           // Do a dummy interpolation from a later timelevel
           // if the desired timelevel does not exist
           int const my_tl = tl >= interp_num_tl ? 0 : tl;
-          assert (my_tl < num_tl);
+          assert (my_tl < num_tl.AT(n));
           
           // Are there enough time levels?
           int const active_tl = CCTK_ActiveTimeLevelsVI (cctkGH, vi);
@@ -1534,7 +1692,7 @@ namespace CarpetInterp {
 
       tmp_output_arrays[tl].resize (N_output_arrays);
       for (int j=0; j<N_output_arrays; ++j) {
-        if (need_time_interp) {
+        if (need_time_interp.AT(j)) {
           if (output_array_type_codes[j] != CCTK_VARIABLE_REAL) {
             CCTK_WARN (CCTK_WARN_ABORT,
                        "time interpolation into output arrays of datatype "
@@ -1580,16 +1738,16 @@ namespace CarpetInterp {
     } // for tl
 
     // Interpolate in time, if neccessary
-    if (need_time_interp) {
-
-      for (int j=0; j<N_output_arrays; ++j) {
+    for (int j=0; j<N_output_arrays; ++j) {
+      if (need_time_interp.AT(j)) {
 
         // Find input array for this output array
         int const n = operand_indices.AT(j);
         CCTK_INT const deriv_order = time_deriv_order.AT(j);
 
-        int const interp_num_tl = interp_num_time_levels.AT(n) > 0 ?
-                                  interp_num_time_levels.AT(n) : num_tl;
+        int const interp_num_tl =
+          interp_num_time_levels.AT(n) > 0 ?
+          interp_num_time_levels.AT(n) : num_tl.AT(n);
         const InterpolationTimes times (interp_num_tl);
         const InterpolationWeights tfacs (deriv_order, times, current_time,
                                           delta_time);
@@ -1606,16 +1764,14 @@ namespace CarpetInterp {
             dest += tfacs[tl] * src;
           }
         }
-
-      } // for j
-
-      for (int tl=0; tl<num_tl; ++tl) {
-        for (int j=0; j<N_output_arrays; ++j) {
+        
+        for (int tl=0; tl<max_num_tl; ++tl) {
           delete [] (CCTK_REAL *)tmp_output_arrays[tl][j];
         }
-      }
 
-    } // if need_time_interp
+      } // if need_time_interp
+    } // for j
+
   }
 
 } // namespace CarpetInterp
diff --git a/Carpet/CarpetLib/interface.ccl b/Carpet/CarpetLib/interface.ccl
index 98a160bf9..44b7729d5 100644
--- a/Carpet/CarpetLib/interface.ccl
+++ b/Carpet/CarpetLib/interface.ccl
@@ -7,6 +7,7 @@ includes header: dist.hh in dist.hh
 
 includes header: bbox.hh in bbox.hh
 includes header: bboxset.hh in bboxset.hh
+includes header: fulltree.hh in fulltree.hh
 includes header: region.hh in region.hh
 includes header: vect.hh in vect.hh
 
diff --git a/Carpet/CarpetLib/src/defs.cc b/Carpet/CarpetLib/src/defs.cc
index 13ef01301..357061ee2 100644
--- a/Carpet/CarpetLib/src/defs.cc
+++ b/Carpet/CarpetLib/src/defs.cc
@@ -224,6 +224,7 @@ ostream& output (ostream& os, const vector<T>& v) {
 #include "bbox.hh"
 #include "bboxset.hh"
 #include "dh.hh"
+#include "fulltree.hh"
 #include "gdata.hh"
 #include "ggf.hh"
 #include "region.hh"
@@ -245,6 +246,7 @@ template size_t memoryof (vector<int> const & v);
 template size_t memoryof (vector<CCTK_REAL> const & v);
 template size_t memoryof (vector<bbox<int,3> > const & v);
 template size_t memoryof (vector<vect<int,3> > const & v);
+template size_t memoryof (vector<fulltree <int,3,pseudoregion_t> *> const & f);
 template size_t memoryof (vector<pseudoregion_t> const & v);
 template size_t memoryof (vector<region_t> const & v);
 template size_t memoryof (vector<sendrecv_pseudoregion_t> const & v);
diff --git a/Carpet/CarpetLib/src/defs.hh b/Carpet/CarpetLib/src/defs.hh
index 59d22f4c8..1152c154a 100644
--- a/Carpet/CarpetLib/src/defs.hh
+++ b/Carpet/CarpetLib/src/defs.hh
@@ -57,14 +57,19 @@ const int dim = 3;
 
 
 // Some shortcuts for type names
+template<typename T, int D> class vect;
 template<typename T, int D> class bbox;
 template<typename T, int D> class bboxset;
-template<typename T, int D> class vect;
+template<typename T, int D, typename P> class fulltree;
+
+struct pseudoregion_t;
+struct region_t;
 
 typedef vect<bool,dim>   bvect;
 typedef vect<int,dim>    ivect;
 typedef bbox<int,dim>    ibbox;
 typedef bboxset<int,dim> ibset;
+typedef fulltree<int,dim,pseudoregion_t> ipfulltree;
 
 // (Try to replace these by b2vect and i2vect)
 typedef vect<vect<bool,2>,dim> bbvect;
diff --git a/Carpet/CarpetLib/src/fulltree.cc b/Carpet/CarpetLib/src/fulltree.cc
new file mode 100644
index 000000000..b114dee19
--- /dev/null
+++ b/Carpet/CarpetLib/src/fulltree.cc
@@ -0,0 +1,492 @@
+#include <algorithm>
+#include <cmath>
+#include <iostream>
+
+#include "defs.hh"
+#include "region.hh"
+
+#include "fulltree.hh"
+
+
+
+#if 0
+// Create an empty tree
+template <typename T, int D, typename P>
+fulltree<T,D,P>::fulltree ()
+  : type (type_empty)
+{
+  assert (invariant());
+}
+#endif
+
+
+
+// Create a tree branch from a list of bounds and subtrees
+template <typename T, int D, typename P>
+fulltree<T,D,P>::fulltree (int dir_, vector <T> const & bounds_,
+                           vector <fulltree *> const & subtrees_)
+  : type (type_branch), dir (dir_), bounds (bounds_), subtrees (subtrees_)
+{
+  assert (dir>=0 and dir<D);
+  assert (bounds.size() == subtrees.size() + 1);
+  assert (invariant());
+}
+
+
+
+// Create a tree leaf from a payload
+template <typename T, int D, typename P>
+fulltree<T,D,P>::fulltree (P const & p_)
+  : type (type_leaf), p (p_)
+{
+  assert (invariant());
+}
+  
+
+
+// Create a tree as copy from another tree
+template <typename T, int D, typename P>
+fulltree<T,D,P>::fulltree (fulltree const & t)
+  : type (t.type)
+{
+  switch (type) {
+  case type_empty:
+    // do nothing
+    break;
+  case type_branch:
+    dir = t.dir;
+    bounds = t.bounds;
+    subtrees.resize (t.subtrees.size());
+    for (size_t i=0; i<subtrees.size(); ++i) {
+      subtrees.AT(i) = new fulltree (*t.subtrees.AT(i));
+    }
+    break;
+  case type_leaf:
+    p = t.p;
+    break;
+  default:
+    assert (0);
+  }
+  assert (invariant());
+}
+  
+
+
+// Copy a tree
+template <typename T, int D, typename P>
+fulltree<T,D,P> &
+fulltree<T,D,P>::operator= (fulltree const & t)
+{
+  assert (invariant());
+  if (is_branch()) {
+    for (size_t i=0; i<subtrees.size(); ++i) {
+      delete subtrees.AT(i);
+    }
+  }
+  type = t.type;
+  switch (type) {
+  case type_empty:
+    // do nothing
+    break;
+  case type_branch:
+    dir = t.dir;
+    bounds = t.bounds;
+    subtrees.resize (t.subtrees.size());
+    for (size_t i=0; i<subtrees.size(); ++i) {
+      subtrees.AT(i) = new fulltree (*t.subtrees.AT(i));
+    }
+    break;
+  case type_leaf:
+    p = t.p;
+    break;
+  default:
+    assert (0);
+  }
+  assert (invariant());
+  return *this;
+}
+
+
+
+// Delete a tree (and its subtrees)
+template <typename T, int D, typename P>
+fulltree<T,D,P>::~fulltree ()
+{
+  assert (invariant());
+  if (is_branch()) {
+    for (size_t i=0; i<subtrees.size(); ++i) {
+      delete subtrees.AT(i);
+    }
+  }
+}
+
+
+
+// Compare trees
+template <typename T, int D, typename P>
+bool
+fulltree<T,D,P>::operator== (fulltree const & t) const
+{
+  assert (invariant());
+  assert (t.invariant());
+  if (type != t.type) return false;
+  switch (type) {
+  case type_empty:
+    break;
+  case type_branch:
+    if (dir != t.dir) return false;
+    if (bounds != t.bounds) return false;
+    if (subtrees != t.subtrees) return false;
+    break;
+  case type_leaf:
+    return p == t.p;
+  default:
+    assert (0);
+  }
+  return true;
+}
+
+
+
+// Invariant
+template <typename T, int D, typename P>
+bool
+fulltree<T,D,P>::invariant () const
+{
+  return empty() or is_branch() or is_leaf();
+}
+
+
+
+// Find the leaf payload corresponding to a position
+template <typename T, int D, typename P>
+P const *
+fulltree<T,D,P>::search (tvect const & ipos) const
+{
+  assert (not empty());
+  // cout << "fulltree::search ipos=" << ipos << endl;
+  if (is_leaf()) return & p;
+  int const i = asearch (ipos[dir], bounds);
+  // cout << "fulltree::search i=" << i << " size=" << subtrees.size() << endl;
+  if (i<0 or i>=int(subtrees.size())) return NULL; // not found
+  return subtrees.AT(i)->search(ipos);
+}
+
+template <typename T, int D, typename P>
+P *
+fulltree<T,D,P>::search (tvect const & ipos)
+{
+  assert (not empty());
+  // cout << "fulltree::search ipos=" << ipos << endl;
+  if (is_leaf()) return & p;
+  int const i = asearch (ipos[dir], bounds);
+  // cout << "fulltree::search i=" << i << " size=" << subtrees.size() << endl;
+  if (i<0 or i>=int(subtrees.size())) return NULL; // not found
+  return subtrees.AT(i)->search(ipos);
+}
+
+
+
+// Const iterator
+template <typename T, int D, typename P>
+fulltree<T,D,P>::const_iterator::const_iterator (fulltree const & f_)
+  : f(f_), i(0), it(0)
+{
+  if (f.is_branch()) {
+    assert (f.subtrees.size() > 0);
+    it = new const_iterator (* f.subtrees.at(i));
+  }
+}
+
+template <typename T, int D, typename P>
+fulltree<T,D,P>::const_iterator::const_iterator (fulltree const & f_, int)
+  : f(f_), it(0)
+{
+  if (f.empty()) {
+    i = 0;
+  } else if (f.is_leaf()) {
+    i = 1;
+  } else {
+    i = f.subtrees.size();
+  }
+}
+
+template <typename T, int D, typename P>
+fulltree<T,D,P>::const_iterator::~const_iterator ()
+{
+  if (it) {
+    delete it;
+  }
+}
+
+template <typename T, int D, typename P>
+fulltree<T,D,P> const & 
+fulltree<T,D,P>::const_iterator::operator* () const
+{
+  assert (not f.empty());
+  if (f.is_leaf()) {
+    return f;
+  } else {
+    return **it;
+  }
+}
+
+template <typename T, int D, typename P>
+bool
+fulltree<T,D,P>::const_iterator::operator== (const_iterator const & it2) const
+{
+  // assert (f == it2.f);
+  assert (&f == &it2.f);
+  if (i != it2.i) return false;
+  if (it == 0 and it2.it == 0) return true;
+  if (it == 0 or it2.it == 0) return false;
+  return *it == *it2.it;
+}
+
+template <typename T, int D, typename P>
+typename fulltree<T,D,P>::const_iterator &
+fulltree<T,D,P>::const_iterator::operator++ ()
+{
+  assert (not done());
+  assert (not f.empty());
+  if (f.is_leaf()) {
+    ++ i;
+  } else {
+    ++ *it;
+    if ((*it).done()) {
+      delete it;
+      it = 0;
+      ++ i;
+      if (not done()) {
+        // to do: use a new function "reset iterator" instead
+        it = new const_iterator (* f.subtrees.at(i));
+      }
+    }
+  }
+  return *this;
+}
+
+template <typename T, int D, typename P>
+bool
+fulltree<T,D,P>::const_iterator::done () const
+{
+  if (f.empty()) {
+    return true;
+  } else if (f.is_leaf()) {
+    return i > 0;
+  } else {
+    return i == f.subtrees.size();
+  }
+}
+
+
+
+// Non-const iterator
+template <typename T, int D, typename P>
+fulltree<T,D,P>::iterator::iterator (fulltree & f_)
+  : f(f_), i(0), it(0)
+{
+  if (f.is_branch()) {
+    assert (f.subtrees.size() > 0);
+    it = new iterator (* f.subtrees.at(i));
+  }
+}
+
+template <typename T, int D, typename P>
+fulltree<T,D,P>::iterator::iterator (fulltree & f_, int)
+  : f(f_), it(0)
+{
+  if (f.empty()) {
+    i = 0;
+  } else if (f.is_leaf()) {
+    i = 1;
+  } else {
+    i = f.subtrees.size();
+  }
+}
+
+template <typename T, int D, typename P>
+fulltree<T,D,P>::iterator::~iterator ()
+{
+  if (it) {
+    delete it;
+  }
+}
+
+template <typename T, int D, typename P>
+fulltree<T,D,P> & 
+fulltree<T,D,P>::iterator::operator* ()
+{
+  assert (not f.empty());
+  if (f.is_leaf()) {
+    return f;
+  } else {
+    return **it;
+  }
+}
+
+template <typename T, int D, typename P>
+bool
+fulltree<T,D,P>::iterator::operator== (iterator const & it2) const
+{
+  // assert (f == it2.f);
+  assert (&f == &it2.f);
+  if (i != it2.i) return false;
+  if (it == 0 and it2.it == 0) return true;
+  if (it == 0 or it2.it == 0) return false;
+  return *it == *it2.it;
+}
+
+template <typename T, int D, typename P>
+typename fulltree<T,D,P>::iterator &
+fulltree<T,D,P>::iterator::operator++ ()
+{
+  assert (not done());
+  assert (not f.empty());
+  if (f.is_leaf()) {
+    ++ i;
+  } else {
+    ++ *it;
+    if ((*it).done()) {
+      delete it;
+      it = 0;
+      ++ i;
+      if (not done()) {
+        // to do: use a new function "reset iterator" instead
+        it = new iterator (* f.subtrees.at(i));
+      }
+    }
+  }
+  return *this;
+}
+
+template <typename T, int D, typename P>
+bool
+fulltree<T,D,P>::iterator::done () const
+{
+  if (f.empty()) {
+    return true;
+  } else if (f.is_leaf()) {
+    return i > 0;
+  } else {
+    return i == f.subtrees.size();
+  }
+}
+
+
+
+// Memory usage
+template <typename T, int D, typename P>
+size_t
+fulltree<T,D,P>::memory () const
+{
+  size_t size = sizeof *this;
+  if (is_branch()) {
+    size += memoryof (bounds) + memoryof (subtrees);
+    for (typename vector <fulltree *>::const_iterator
+           i = subtrees.begin(); i != subtrees.end(); ++ i)
+    {
+      size += memoryof(**i);
+    }
+  }
+  return size;
+}
+
+
+
+// Output helper
+template <typename T, int D, typename P>
+void
+fulltree<T,D,P>::output (ostream & os) const
+{
+  os << "fulltree{";
+  if (empty()) {
+    os << "empty";
+  } else if (is_branch()) {
+    os << "branch:"
+       << "dir=" << dir << ","
+       << "subtrees=[";
+    for (size_t i=0; i<subtrees.size(); ++i) {
+      os << bounds.at(i) << ":[" << i << "]=" << subtrees.at(i) << ":";
+    }
+    os << bounds.at(subtrees.size()) << "]";
+  } else {
+    os << "leaf:"
+       << "payload=" << p;
+  }
+  os << "}";
+}
+
+
+
+// Generic arithmetic search
+template <typename T>
+static
+int asearch (T const t, vector <T> const & ts)
+{
+  // cout << "fulltree::asearch t=" << t << " ts=" << ts << endl;
+  int imin = 0;
+  int imax = int(ts.size()) - 1;
+  if (imax <= imin) {
+    // cout << "fulltree::asearch: no values" << endl;
+    return imin;
+  }
+  T tmin = ts.AT(imin);
+  // cout << "fulltree::asearch: imin=" << imin << " tmin=" << tmin << endl;
+  if (t < tmin) {
+    // cout << "fulltree::asearch: below minimum" << endl;
+    return -1;
+  }
+  T tmax = ts.AT(imax);
+  // cout << "fulltree::asearch: imax=" << imax << " tmax=" << tmax << endl;
+  if (t >= tmax) {
+    // cout << "fulltree::asearch: above maximum" << endl;
+    return imax;
+  }
+  int isize = imax - imin;
+  for (;;) {
+    // check loop invariants
+    assert (imin < imax);
+    assert (t >= tmin);
+    assert (t < tmax);
+    // cout << "fulltree::asearch t=" << t << " imin=" << imin << " imax=" << imax << endl;
+    if (imax == imin+1) {
+      // cout << "fulltree::asearch: found value" << endl;
+      return imin;
+    }
+    assert (tmax > tmin);       // require that ts is strictly ordered
+    CCTK_REAL const rguess =
+      (imax - imin) * CCTK_REAL(t - tmin) / (tmax - tmin);
+    int const iguess = imin + max (1, int(floor(rguess)));
+    // handle round-off errors
+    if (iguess == imax) {
+      // cout << "fulltree::asearch: accidental hit after roundoff error" << endl;
+      return imax - 1;
+    }
+    assert (iguess > imin and iguess < imax);
+    T const tguess = ts.AT(iguess);
+    // cout << "fulltree::asearch: intersecting at iguess=" << iguess << " tguess=" << tguess << endl;
+    if (t < tguess) {
+      imax = iguess;
+      tmax = tguess;
+      // cout << "fulltree::asearch: new imax=" << imax << " tmax=" << tmax << endl;
+    } else {
+      imin = iguess;
+      tmin = tguess;
+      // cout << "fulltree::asearch: new imin=" << imin << " tmin=" << tmin << endl;
+    }
+    // check loop variant
+    int const newisize = imax - imin;
+    assert (newisize < isize);
+    isize = newisize;
+  }
+}
+
+
+
+template class fulltree <int, dim, pseudoregion_t>;
+
+template size_t memoryof (fulltree <int, dim, pseudoregion_t> const & f);
+template size_t memoryof (fulltree <int, dim, pseudoregion_t> * const & f);
+
+template ostream & operator<< (ostream & os, fulltree <int, dim, pseudoregion_t> const & f);
diff --git a/Carpet/CarpetLib/src/fulltree.hh b/Carpet/CarpetLib/src/fulltree.hh
new file mode 100644
index 000000000..6ffca209e
--- /dev/null
+++ b/Carpet/CarpetLib/src/fulltree.hh
@@ -0,0 +1,181 @@
+#ifndef FULLTREE_HH
+#define FULLTREE_HH
+
+#include <cassert>
+#include <cmath>
+#include <cstdlib>
+#include <iostream>
+#include <vector>
+
+#include <vect.hh>
+
+using namespace std;
+
+
+
+// This is a "full tree" data structure, i.e., a tree data structure
+// which decomposes a cuboid domain into a set of non-overlapping
+// cuboid subdomains.  It is an n-ary tree, i.e., each tree node can
+// have arbitrarily many subtrees.  Each node splits a domain in
+// exactly one direction.  Subdomains cannot be empty.
+
+// All intervals are closed at the lower and open at the upper
+// boundary.  This makes it possible to combine adjacent such
+// intervals, obtaining again an interval with this property.  (In
+// particular, if bboxes are stored, the upper bound of bboxes must be
+// increased by the stride to arrive at such intervals.)
+
+
+// Generic arithmetic search
+template <typename T>
+static
+int asearch (T t, vector <T> const & ts);
+
+
+
+// T: index space (usually integer, or maybe real)
+// D: number of dimensions (rank)
+// P: payload (e.g. region_t)
+template <typename T, int D, typename P>
+class fulltree {
+  
+public:
+  
+  // Short name for a small vector
+  typedef vect <T, D> tvect;
+  
+private:
+  
+  enum tree_type_t { type_empty, type_branch, type_leaf } type;
+  
+  // Direction in which the node is split (0 <= dir < D)
+  int dir;
+  
+  // If this is a branch:
+  // n+1 bounds, splitting the domain
+  vector <T> bounds;            // [n+1]
+  // n pointers to subtrees
+  vector <fulltree *> subtrees; // [n]
+  
+  // If this is a leaf:
+  // just the payload
+  P p;
+  
+public:
+  
+#if 0
+  // Create an empty tree
+  fulltree ();
+#endif
+  
+  // Create a tree branch from a list of bounds and subtrees
+  fulltree (int dir_, vector <T> const & bounds_,
+            vector <fulltree *> const & subtrees_);
+  
+  // Create a tree leaf from a payload
+  fulltree (P const & p_);
+  
+  // Create a tree as copy from another tree
+  fulltree (fulltree const & t);
+  
+  // Copy a tree
+  fulltree & operator= (fulltree const & t);
+  
+  // Delete a tree (and its subtrees)
+  ~fulltree ();
+  
+  // Inquire tree properties
+  bool empty() const { return type == type_empty; }
+  bool is_branch() const { return type == type_branch; }
+  bool is_leaf() const { return type == type_leaf; }
+  
+  // Compare trees
+  bool operator== (fulltree const & t) const;
+  bool operator!= (fulltree const & t) const
+  { return not (*this == t); }
+  
+  // Invariant
+  bool invariant () const;
+  
+  // Access the payload
+  P const & payload () const { assert (is_leaf()); return p; }
+  P & payload () { assert (is_leaf()); return p; }
+  
+  // Find the leaf payload corresponding to a position
+  P const * search (tvect const & ipos) const;
+  P * search (tvect const & ipos);
+  
+  class const_iterator {
+    fulltree const & f;
+    size_t i;
+    const_iterator * it;
+  public:
+    const_iterator (fulltree const & f_);
+    const_iterator (fulltree const & f_, int);
+    ~const_iterator ();
+    fulltree const & operator* () const;
+    bool operator== (const_iterator const & it2) const;
+    bool operator!= (const_iterator const & it2) const
+    { return not (*this == it2); }
+    const_iterator & operator++ ();
+    bool done () const;
+  };
+  
+  class iterator {
+    fulltree & f;
+    size_t i;
+    iterator * it;
+  public:
+    iterator (fulltree & f_);
+    iterator (fulltree & f_, int);
+    ~iterator ();
+    fulltree & operator* ();
+    bool operator== (iterator const & it2) const;
+    bool operator!= (iterator const & it2) const
+    { return not (*this == it2); }
+    iterator & operator++ ();
+    bool done () const;
+  };
+  
+  const_iterator begin() const
+  { return const_iterator (*this); }
+  
+  const_iterator end() const
+  { return const_iterator (*this, 0); }
+  
+  iterator begin()
+  { return iterator (*this); }
+  
+  iterator end()
+  { return iterator (*this, 0); }
+  
+  // Memory usage
+  size_t memory () const;
+  
+  // Output helper
+  void output (ostream & os) const;
+};
+
+
+
+// Memory usage
+template <typename T, int D, typename P>
+inline size_t memoryof (fulltree<T,D,P> const & f) { return f.memory(); }
+
+template <typename T, int D, typename P>
+inline size_t memoryof (fulltree<T,D,P> * const & fp) { return sizeof fp; }
+
+
+
+// Output
+template <typename T, int D, typename P>
+ostream &
+operator<< (ostream & os, fulltree<T,D,P> const & f)
+{
+  f.output (os);
+  return os;
+}
+
+
+
+#endif  // #ifndef FULLTREE_HH
diff --git a/Carpet/CarpetLib/src/ggf.cc b/Carpet/CarpetLib/src/ggf.cc
index b6eefa499..ae42668d8 100644
--- a/Carpet/CarpetLib/src/ggf.cc
+++ b/Carpet/CarpetLib/src/ggf.cc
@@ -512,8 +512,8 @@ transfer_from_all (comm_state & state,
     pseudoregion_t const & precv = (* ipsendrecv).recv;
     ibbox const & send = psend.extent;
     ibbox const & recv = precv.extent;
-    int const c2 = psend.processor;
-    int const c1 = precv.processor;
+    int const c2 = psend.component;
+    int const c1 = precv.component;
     
     // Source and destination data
     gdata * const dst = storage.AT(ml1).AT(rl1).AT(c1).AT(tl1);
diff --git a/Carpet/CarpetLib/src/ggf.hh b/Carpet/CarpetLib/src/ggf.hh
index e3c351b9d..b2c86b8db 100644
--- a/Carpet/CarpetLib/src/ggf.hh
+++ b/Carpet/CarpetLib/src/ggf.hh
@@ -62,7 +62,7 @@ protected:
 public:
   const int vectorlength;       // vector length
   const int vectorindex;        // index of *this
-  const ggf* vectorleader;      // first vector element
+  ggf* const vectorleader;      // first vector element
   
 private:
   mdata oldstorage;             // temporary storage
diff --git a/Carpet/CarpetLib/src/gh.cc b/Carpet/CarpetLib/src/gh.cc
index fa9ff2125..9b4a5f443 100644
--- a/Carpet/CarpetLib/src/gh.cc
+++ b/Carpet/CarpetLib/src/gh.cc
@@ -15,6 +15,8 @@
 
 using namespace std;
 
+using namespace CarpetLib;
+
 
 
   // Constructors
@@ -76,10 +78,12 @@ gh::
 // Modifiers
 void
 gh::
-regrid (mregs const & regs)
+regrid (rregs const & superregs, mregs const & regs)
 {
   DECLARE_CCTK_PARAMETERS;
   
+  superregions = superregs;
+  
   // Save the grid hierarchy
   oldregions.clear ();
   swap (oldregions, regions);
diff --git a/Carpet/CarpetLib/src/gh.hh b/Carpet/CarpetLib/src/gh.hh
index 0a283b4cf..ffac38804 100644
--- a/Carpet/CarpetLib/src/gh.hh
+++ b/Carpet/CarpetLib/src/gh.hh
@@ -47,6 +47,10 @@ public:				// should be readonly
   vector<vector<ibbox> > baseextents; // [ml][rl]
   const i2vect boundary_width;
   
+  // Extents of the regions before distributing them over the
+  // processors
+  rregs superregions;
+  
   mregs regions;                // extents and properties of all grids
   mregs oldregions;             // a copy, used during regridding
   
@@ -65,7 +69,7 @@ public:
   ~gh ();
   
   // Modifiers
-  void regrid (mregs const & regs);
+  void regrid (rregs const & superregs, mregs const & regs);
   bool recompose (int rl, bool do_prolongate);
   
 private:
diff --git a/Carpet/CarpetLib/src/make.code.defn b/Carpet/CarpetLib/src/make.code.defn
index aa9b399c0..88f6261ce 100644
--- a/Carpet/CarpetLib/src/make.code.defn
+++ b/Carpet/CarpetLib/src/make.code.defn
@@ -8,6 +8,7 @@ SRCS =	bbox.cc					\
 	defs.cc					\
 	dh.cc					\
 	dist.cc					\
+	fulltree.cc				\
 	gdata.cc				\
 	gf.cc					\
 	ggf.cc					\
diff --git a/Carpet/CarpetLib/src/region.cc b/Carpet/CarpetLib/src/region.cc
index caf535598..0230d373d 100644
--- a/Carpet/CarpetLib/src/region.cc
+++ b/Carpet/CarpetLib/src/region.cc
@@ -1,5 +1,8 @@
+#include <cassert>
+#include <cstdlib>
 #include <iostream>
 
+#include "bboxset.hh"
 #include "defs.hh"
 #include "region.hh"
 
@@ -7,6 +10,69 @@ using namespace std;
 
 
 
+region_t::region_t ()
+  : processor (-1), processors (NULL)
+{
+  assert (invariant());
+}
+
+region_t::region_t (region_t const & a)
+{
+  assert (a.invariant());
+  extent = a.extent;
+  outer_boundaries = a.outer_boundaries;
+  map = a.map;
+  processor = a.processor;
+  if (a.processors == NULL) {
+    processors = NULL;
+  } else {
+    processors = new ipfulltree (*a.processors);
+  }
+  assert (invariant());
+}
+
+region_t &
+region_t::operator= (region_t const & a)
+{
+  assert (invariant());
+  if (processors != NULL) {
+    delete processors;
+  }
+  assert (a.invariant());
+  extent = a.extent;
+  outer_boundaries = a.outer_boundaries;
+  map = a.map;
+  processor = a.processor;
+  if (a.processors == NULL) {
+    processors = NULL;
+  } else {
+    processors = new ipfulltree (*a.processors);
+  }
+  assert (invariant());
+  return *this;
+}
+
+
+region_t::~region_t ()
+{
+  assert (invariant());
+  if (processors != NULL) {
+    delete processors;
+  }
+}
+
+
+
+bool
+region_t::invariant () const
+{
+  if (processor >= 0 and processors != NULL) return false;
+  return true;
+}
+
+
+
+// Compare two regions for equality.
 bool
 operator== (region_t const & a, region_t const & b)
 {
@@ -14,7 +80,113 @@ operator== (region_t const & a, region_t const & b)
     a.extent == b.extent and
     all (all (a.outer_boundaries == b.outer_boundaries)) and
     a.map == b.map and
-    a.processor == b.processor;
+    a.processor == b.processor and
+    ((a.processors == NULL and b.processors == NULL) or
+     (a.processors != NULL and b.processors != NULL and
+      *a.processors == *b.processors));
+}
+
+
+
+// Combine a collection of regions.  Regions can be combined if they
+// abutt on boundaries which are not outer boundaries, ignoring the
+// processor distribution.  This should lead to a canonical
+// representations of collections of regions.
+//
+// We use vectors to represent the collection, but we could also use
+// other containers.  oldregs is read, newregs is added-to.  newregs
+// is not cleared.
+void
+combine_regions (vector<region_t> const & oldregs,
+                 vector<region_t> & newregs)
+{
+  // Find the union of all regions' bounding boxes, and the union of
+  // all regions' outer boundaries.  Represent the boundaries as the
+  // outermost layer of grid points of the corresponding bounding
+  // boxes.
+  int const m = oldregs.empty() ? -1 : oldregs.AT(0).map;
+  ibset comps;
+  ibset cobnds[2][dim];
+  for (vector<region_t>::const_iterator
+         ri = oldregs.begin(); ri != oldregs.end(); ++ ri)
+  {
+    region_t const & reg = * ri;
+    assert (reg.map == m);
+    comps += reg.extent;
+    for (int f = 0; f < 2; ++ f) {
+      for (int d = 0; d < dim; ++ d) {
+        if (reg.outer_boundaries[f][d]) {
+          ibbox bnd = reg.extent;
+          ivect lo = bnd.lower();
+          ivect up = bnd.upper();
+          if (f==0) {
+            up[d] = lo[d];
+          } else {
+            lo[d] = up[d];
+          }
+          bnd = ibbox (lo, up, bnd.stride());
+          cobnds[f][d] += bnd;
+        }
+      }
+    }
+  }
+  comps.normalize();
+  for (int f = 0; f < 2; ++ f) {
+    for (int d = 0; d < dim; ++ d) {
+      cobnds[f][d].normalize();
+    }
+  }
+  
+  // Reserve (generous) memory for the result
+  size_t const needsize = newregs.size() + comps.setsize();
+  if (newregs.capacity() < needsize) {
+    newregs.reserve (1000 + 2 * needsize);
+  }
+  
+  // Insert the regions
+  for (ibset::const_iterator ci = comps.begin(); ci != comps.end(); ++ ci) {
+    ibbox const & c = * ci;
+    b2vect obnds;
+    for (int f = 0; f < 2; ++ f) {
+      for (int d = 0; d < dim; ++ d) {
+        obnds[f][d] = cobnds[f][d].intersects (c);
+        if (obnds[f][d]) {
+          ivect lo = c.lower();
+          ivect up = c.upper();
+          if (f) lo[d]=up[d]; else up[d]=lo[d];
+          ibbox const cbnds (lo, up, c.stride());
+          if (not ((cobnds[f][d] & ibset(c)) == ibset(cbnds))) {
+            cout << "cobnds[f][d] = " << cobnds[f][d] << endl
+                 << "ibset(c) = " << ibset(c) << endl
+                 << "(cobnds[f][d] & ibset(c)) = " << (cobnds[f][d] & ibset(c)) << endl
+                 << "ibset(cbnds) = " << ibset(cbnds) << endl;
+          }
+          assert ((cobnds[f][d] & ibset(c)) == ibset(cbnds));
+        }
+      }
+    }
+    
+    region_t reg;
+    reg.extent           = c;
+    reg.outer_boundaries = obnds;
+    reg.map              = m;
+    reg.processor        = -1;
+    reg.processors       = NULL;
+    newregs.push_back (reg);
+  }
+}
+
+
+
+size_t memoryof (region_t const & reg)
+{
+  return
+    memoryof (reg.extent) +
+    memoryof (reg.outer_boundaries) +
+    memoryof (reg.map) +
+    memoryof (reg.processor) +
+    memoryof (reg.processors) +
+    (reg.processors != NULL ? memoryof (*reg.processors) : 0);
 }
 
 
@@ -59,6 +231,8 @@ operator>> (istream & is, region_t & reg)
   skipws (is);
   consume (is, ')');
   
+  reg.processors = NULL;
+  
   return is;
 }
 
@@ -77,9 +251,20 @@ operator<< (ostream & os, region_t const & reg)
 
 
 
+// Compare two pseudoregions for equality.
+bool
+operator== (pseudoregion_t const & a, pseudoregion_t const & b)
+{
+  return
+    a.extent == b.extent and
+    a.component == b.component;
+}
+
+
+
 ostream & operator<< (ostream & os, pseudoregion_t const & p)
 {
-  return os << p.extent << "/p:" << p.processor;
+  return os << p.extent << "/c:" << p.component;
 }
 
 ostream & operator<< (ostream & os, sendrecv_pseudoregion_t const & srp)
diff --git a/Carpet/CarpetLib/src/region.hh b/Carpet/CarpetLib/src/region.hh
index 033e73a05..66037bdc7 100644
--- a/Carpet/CarpetLib/src/region.hh
+++ b/Carpet/CarpetLib/src/region.hh
@@ -2,20 +2,29 @@
 #define REGION_HH
 
 #include <iostream>
+#include <vector>
 
 #include "defs.hh"
 #include "bbox.hh"
+#include "fulltree.hh"
 #include "vect.hh"
 
 
 
 // Region description
 struct region_t {
-  ibbox  extent;                            // extent
-  b2vect outer_boundaries;                  // outer boundaries
-  int    map;                               // map to which this
-                                            // region belongs
-  int    processor;                         // processor number
+  ibbox        extent;           // extent
+  b2vect       outer_boundaries; // outer boundaries
+  int          map;              // map to which this region belongs
+  int          processor;        // processor number
+  ipfulltree * processors;       // processor decomposition
+  
+  region_t ();
+  region_t (region_t const & a);
+  region_t & operator= (region_t const & a);
+  ~region_t ();
+  
+  bool invariant () const;
 };
 
 
@@ -29,54 +38,62 @@ bool operator!= (region_t const & a, region_t const & b)
 
 
 
-inline size_t memoryof (region_t const & reg)
-{
-  return
-    memoryof (reg.extent) +
-    memoryof (reg.outer_boundaries) +
-    memoryof (reg.map) +
-    memoryof (reg.processor);
-}
+void
+combine_regions (vector<region_t> const & oldregs,
+                 vector<region_t> & newregs);
+
+
+
+size_t memoryof (region_t const & reg);
 
 istream & operator>> (istream & is, region_t       & reg);
 ostream & operator<< (ostream & os, region_t const & reg);
 
 
 
-// A pseudoregion is almost a region; it is a bbox that lives on a
-// certain processor.  Pseudoregions are a compact way to store
-// information about what processors needs to send data to what other
-// processors during synchronisation or regridding.
+// A pseudoregion is almost a region; it is a bbox that belongs to a
+// certain component.  Pseudoregions are a compact way to store
+// information about what components needs to send data to what other
+// components during synchronisation or regridding.
 struct pseudoregion_t {
   ibbox extent;
-  int processor;
+  int component;
   pseudoregion_t ()
   {
   }
-  pseudoregion_t (ibbox const & extent_, int const processor_)
-    : extent (extent_), processor (processor_)
+  pseudoregion_t (ibbox const & extent_, int const component_)
+    : extent (extent_), component (component_)
   {
   }
 };
 
+bool operator== (pseudoregion_t const & a, pseudoregion_t const & b);
+inline
+bool operator!= (pseudoregion_t const & a, pseudoregion_t const & b)
+{
+  return not (a == b);
+}
+
 inline size_t memoryof (pseudoregion_t const & p)
 {
   return
     memoryof (p.extent) +
-    memoryof (p.processor);
+    memoryof (p.component);
 }
 
 ostream & operator<< (ostream & os, pseudoregion_t const & p);
 
+
+
 struct sendrecv_pseudoregion_t {
   pseudoregion_t send, recv;
   sendrecv_pseudoregion_t ()
   {
   }
-  sendrecv_pseudoregion_t (ibbox const & send_extent, int const send_processor,
-                           ibbox const & recv_extent, int const recv_processor)
-    : send (pseudoregion_t (send_extent, send_processor)),
-      recv (pseudoregion_t (recv_extent, recv_processor))
+  sendrecv_pseudoregion_t (ibbox const & send_extent, int const send_component,
+                           ibbox const & recv_extent,  int const recv_component)
+    : send (pseudoregion_t (send_extent, send_component)),
+      recv (pseudoregion_t (recv_extent, recv_component))
   {
   }
 };
diff --git a/Carpet/CarpetRegrid/interface.ccl b/Carpet/CarpetRegrid/interface.ccl
index 1cece9380..ac9d5b112 100644
--- a/Carpet/CarpetRegrid/interface.ccl
+++ b/Carpet/CarpetRegrid/interface.ccl
@@ -52,9 +52,11 @@ USES FUNCTION ConvertFromPhysicalBoundary
 
 # The true prototype of the routine below:
 #   int Carpet_Regrid (cGH const * cctkGH,
+#                      gh::rregs * superregss,
 #                      gh::mregs * regsss,
 #                      int         force);
 CCTK_INT FUNCTION Carpet_Regrid (CCTK_POINTER_TO_CONST IN cctkGH, \
+                                 CCTK_POINTER IN superregss,      \
                                  CCTK_POINTER IN regsss,          \
 				 CCTK_INT IN force)
 PROVIDES FUNCTION Carpet_Regrid WITH CarpetRegrid_Regrid LANGUAGE C
diff --git a/Carpet/CarpetRegrid/src/automatic.cc b/Carpet/CarpetRegrid/src/automatic.cc
index 16c12063c..442224f2c 100644
--- a/Carpet/CarpetRegrid/src/automatic.cc
+++ b/Carpet/CarpetRegrid/src/automatic.cc
@@ -26,15 +26,12 @@ namespace CarpetRegrid {
   
   int Automatic (cGH const * const cctkGH,
                  gh const & hh,
-                 gh::mregs & regsss)
+                 gh::rregs & regss)
   {
     DECLARE_CCTK_PARAMETERS;
     
     assert (refinement_levels >= 1);
     
-    assert (regsss.size() >= 1);
-    vector<vector<region_t> > regss = regsss.at(0);
-    
     
     
     const int vi = CCTK_VarIndex (errorvar);
@@ -60,9 +57,6 @@ namespace CarpetRegrid {
       (cctkGH, hh, reflevel, min(reflevels+1, (int)refinement_levels),
        errorgf, regs);
     
-    // make multiprocessor aware
-    SplitRegions (cctkGH, regs);
-    
     if (regs.size() == 0) {
       // remove all finer levels
       regss.resize(reflevel+1);
@@ -77,9 +71,6 @@ namespace CarpetRegrid {
       }
     }
     
-    // make multigrid aware
-    MakeMultigridBoxes (cctkGH, Carpet::map, regss, regsss);
-    
     return 1;
   }
   
diff --git a/Carpet/CarpetRegrid/src/baselevel.cc b/Carpet/CarpetRegrid/src/baselevel.cc
index 105b9f72c..a7229fdad 100644
--- a/Carpet/CarpetRegrid/src/baselevel.cc
+++ b/Carpet/CarpetRegrid/src/baselevel.cc
@@ -19,14 +19,12 @@ namespace CarpetRegrid {
   
   int BaseLevel (cGH const * const cctkGH,
                  gh const & hh,
-                 gh::mregs & regsss)
+                 gh::rregs & regss)
   {
     DECLARE_CCTK_PARAMETERS;
     
     assert (refinement_levels == 1);
     
-    assert (regsss.size() == 1);
-    
     return 0;
   }
   
diff --git a/Carpet/CarpetRegrid/src/centre.cc b/Carpet/CarpetRegrid/src/centre.cc
index 5bda85c85..078dd0181 100644
--- a/Carpet/CarpetRegrid/src/centre.cc
+++ b/Carpet/CarpetRegrid/src/centre.cc
@@ -19,7 +19,7 @@ namespace CarpetRegrid {
   
   int Centre (cGH const * const cctkGH,
               gh const & hh,
-              gh::mregs & regsss)
+              gh::rregs & regss)
   {
     DECLARE_CCTK_PARAMETERS;
     
@@ -28,9 +28,6 @@ namespace CarpetRegrid {
     // do nothing if the levels already exist
     if (reflevel == refinement_levels) return 0;
     
-    assert (regsss.size() >= 1);
-    vector<vector<region_t> > regss = regsss.at(0);
-    
     regss.resize (refinement_levels);
     
     bvect const symmetric (symmetry_x, symmetry_y, symmetry_z);
@@ -40,7 +37,7 @@ namespace CarpetRegrid {
     ivect rlb  = hh.baseextents.at(0).at(0).lower();
     ivect rub  = hh.baseextents.at(0).at(0).upper();
     
-    assert (! smart_outer_boundaries);
+    assert (not smart_outer_boundaries);
     
     for (size_t rl=1; rl<regss.size(); ++rl) {
       
@@ -69,16 +66,10 @@ namespace CarpetRegrid {
       
       regs.at(0).map = Carpet::map;
       
-      // make multiprocessor aware
-      SplitRegions (cctkGH, regs);
-      
       regss.at(rl) = regs;
       
     } // for rl
     
-    // make multigrid aware
-    MakeMultigridBoxes (cctkGH, Carpet::map, regss, regsss);
-    
     return 1;
   }
   
diff --git a/Carpet/CarpetRegrid/src/manualcoordinatelist.cc b/Carpet/CarpetRegrid/src/manualcoordinatelist.cc
index 7becee434..f32ec597c 100644
--- a/Carpet/CarpetRegrid/src/manualcoordinatelist.cc
+++ b/Carpet/CarpetRegrid/src/manualcoordinatelist.cc
@@ -23,7 +23,7 @@ namespace CarpetRegrid {
   
   int ManualCoordinateList (cGH const * const cctkGH,
                             gh const & hh,
-                            gh::mregs & regsss)
+                            gh::rregs & regss)
   {
     DECLARE_CCTK_PARAMETERS;
     int ierr;
@@ -33,9 +33,6 @@ namespace CarpetRegrid {
     // do nothing if the levels already exist
     if (reflevel == refinement_levels && !tracking) return 0;
     
-    assert (regsss.size() >= 1);
-    vector<vector<region_t> > regss = regsss.at(0);
-    
     jjvect nboundaryzones, is_internal, is_staggered, shiftout;
     ierr = GetBoundarySpecification
       (2*dim, &nboundaryzones[0][0], &is_internal[0][0],
@@ -116,7 +113,7 @@ namespace CarpetRegrid {
         } // for c
       } // for rl
       
-    } else {                    // if ! smart_outer_boundaries
+    } else {                    // if not smart_outer_boundaries
       
       vector<vector<bbvect> > newobss1;
       if (strcmp(outerbounds, "") != 0) {
@@ -159,7 +156,7 @@ namespace CarpetRegrid {
         }
       }
 
-    } // if ! smart_outer_boundaries
+    } // if not smart_outer_boundaries
     
     for (int rl=1; rl<refinement_levels; ++rl) {
       
@@ -224,17 +221,11 @@ namespace CarpetRegrid {
         }
         
       } // if merge_overlapping_components
-
-      // make multiprocessor aware
-      SplitRegions (cctkGH, regs);
       
       regss.at(rl) = regs;
       
     } // for rl
     
-    // make multigrid aware
-    MakeMultigridBoxes (cctkGH, Carpet::map, regss, regsss);
-    
     return 1;
   }
   
diff --git a/Carpet/CarpetRegrid/src/manualcoordinates.cc b/Carpet/CarpetRegrid/src/manualcoordinates.cc
index e6d2d2ae0..1c4cea5b0 100644
--- a/Carpet/CarpetRegrid/src/manualcoordinates.cc
+++ b/Carpet/CarpetRegrid/src/manualcoordinates.cc
@@ -20,7 +20,7 @@ namespace CarpetRegrid {
   
   int ManualCoordinates (cGH const * const cctkGH,
                          gh const & hh,
-                         gh::mregs & regsss)
+                         gh::rregs & regss)
   {
     DECLARE_CCTK_PARAMETERS;
     
@@ -32,9 +32,6 @@ namespace CarpetRegrid {
     // do nothing if the levels already exist
     if (reflevel == refinement_levels) return 0;
     
-    assert (regsss.size() >= 1);
-    vector<vector<region_t> > regss = regsss.at(0);
-    
     regss.resize (refinement_levels);
     
     vector<rvect> lower(3), upper(3);
@@ -59,16 +56,10 @@ namespace CarpetRegrid {
         (cctkGH, hh, rl, refinement_levels,
          lower.at(rl-1), upper.at(rl-1), reg, regs);
       
-      // make multiprocessor aware
-      SplitRegions (cctkGH, regs);
-      
       regss.at(rl) = regs;
       
     } // for rl
     
-    // make multigrid aware
-    MakeMultigridBoxes (cctkGH, Carpet::map, regss, regsss);
-    
     return 1;
   }
   
diff --git a/Carpet/CarpetRegrid/src/manualgridpointlist.cc b/Carpet/CarpetRegrid/src/manualgridpointlist.cc
index 60968e59d..fe5bdb4ab 100644
--- a/Carpet/CarpetRegrid/src/manualgridpointlist.cc
+++ b/Carpet/CarpetRegrid/src/manualgridpointlist.cc
@@ -23,7 +23,7 @@ namespace CarpetRegrid {
   
   int ManualGridpointList (cGH const * const cctkGH,
                            gh const & hh,
-                           gh::mregs & regsss)
+                           gh::rregs & regss)
   {
     DECLARE_CCTK_PARAMETERS;
     
@@ -32,9 +32,6 @@ namespace CarpetRegrid {
     // do nothing if the levels already exist
     if (reflevel == refinement_levels) return 0;
     
-    assert (regsss.size() >= 1);
-    vector<vector<region_t> > regss = regsss.at(0);
-    
     regss.resize (refinement_levels);
     
     vector<vector<ibbox> > newbbss;
@@ -116,16 +113,10 @@ namespace CarpetRegrid {
            ext.lower(), ext.upper(), reg, regs);
       }
       
-      // make multiprocessor aware
-      SplitRegions (cctkGH, regs);
-      
       regss.at(rl) = regs;
       
     } // for rl
     
-    // make multigrid aware
-    MakeMultigridBoxes (cctkGH, Carpet::map, regss, regsss);
-    
     return 1;
   }
   
diff --git a/Carpet/CarpetRegrid/src/manualgridpoints.cc b/Carpet/CarpetRegrid/src/manualgridpoints.cc
index 6481247c9..3df1f2119 100644
--- a/Carpet/CarpetRegrid/src/manualgridpoints.cc
+++ b/Carpet/CarpetRegrid/src/manualgridpoints.cc
@@ -22,7 +22,7 @@ namespace CarpetRegrid {
   
   int ManualGridpoints (cGH const * const cctkGH,
                         gh const & hh,
-                        gh::mregs & regsss)
+                        gh::rregs & regss)
   {
     DECLARE_CCTK_PARAMETERS;
     
@@ -34,9 +34,6 @@ namespace CarpetRegrid {
     // do nothing if the levels already exist
     if (reflevel == refinement_levels) return 0;
     
-    assert (regsss.size() >= 1);
-    vector<vector<region_t> > regss = regsss.at(0);
-    
     regss.resize (refinement_levels);
     
     vector<ivect> ilower(3), iupper(3);
@@ -47,7 +44,7 @@ namespace CarpetRegrid {
     ilower.at(2) = ivect (l3ixmin, l3iymin, l3izmin);
     iupper.at(2) = ivect (l3ixmax, l3iymax, l3izmax);
     
-    assert (! smart_outer_boundaries);
+    assert (not smart_outer_boundaries);
     
     for (size_t rl=1; rl<regss.size(); ++rl) {
       
@@ -61,16 +58,10 @@ namespace CarpetRegrid {
         (cctkGH, hh, rl,refinement_levels,
          ilower.at(rl-1), iupper.at(rl-1), reg, regs);
       
-      // make multiprocessor aware
-      SplitRegions (cctkGH, regs);
-      
       regss.at(rl) = regs;
       
     } // for rl
     
-    // make multigrid aware
-    MakeMultigridBoxes (cctkGH, Carpet::map, regss, regsss);
-    
     return 1;
   }
   
diff --git a/Carpet/CarpetRegrid/src/moving.cc b/Carpet/CarpetRegrid/src/moving.cc
index 41980ea45..075fdbb93 100644
--- a/Carpet/CarpetRegrid/src/moving.cc
+++ b/Carpet/CarpetRegrid/src/moving.cc
@@ -20,16 +20,13 @@ namespace CarpetRegrid {
   
   int Moving (cGH const * const cctkGH,
               gh const & hh,
-              gh::mregs & regsss)
+              gh::rregs & regss)
   {
     DECLARE_CCTK_ARGUMENTS;
     DECLARE_CCTK_PARAMETERS;
     
     assert (refinement_levels >= 1);
     
-    assert (regsss.size() >= 1);
-    vector<vector<region_t> > regss = regsss.at(0);
-    
     regss.resize (refinement_levels);
     
     bvect const symmetric (symmetry_x, symmetry_y, symmetry_z);
@@ -60,16 +57,10 @@ namespace CarpetRegrid {
       ManualCoordinates_OneLevel
         (cctkGH, hh, rl, refinement_levels, rlb, rub, reg, regs);
       
-      // make multiprocessor aware
-      SplitRegions (cctkGH, regs);
-      
       regss.at(rl) = regs;
       
     } // for rl
     
-    // make multigrid aware
-    MakeMultigridBoxes (cctkGH, Carpet::map, regss, regsss);
-    
     return 1;
   }
   
diff --git a/Carpet/CarpetRegrid/src/regrid.cc b/Carpet/CarpetRegrid/src/regrid.cc
index 2dd2a33a9..0a5567e2b 100644
--- a/Carpet/CarpetRegrid/src/regrid.cc
+++ b/Carpet/CarpetRegrid/src/regrid.cc
@@ -19,6 +19,7 @@ namespace CarpetRegrid {
   using namespace Carpet;
   
   CCTK_INT CarpetRegrid_Regrid (CCTK_POINTER_TO_CONST const cctkGH_,
+                                CCTK_POINTER const superregss_,
                                 CCTK_POINTER const regsss_,
 				CCTK_INT force)
   {
@@ -26,6 +27,7 @@ namespace CarpetRegrid {
     
     const cGH * const cctkGH = (const cGH *) cctkGH_;
     
+    gh::rregs & superregss = * (gh::rregs *) superregss_;
     gh::mregs & regsss = * (gh::mregs *) regsss_;
     
     gh const & hh = *vhh.at(Carpet::map);
@@ -151,44 +153,54 @@ namespace CarpetRegrid {
     
     if (CCTK_EQUALS(refined_regions, "none")) {
       
-      do_recompose = BaseLevel (cctkGH, hh, regsss);
+      do_recompose = BaseLevel (cctkGH, hh, superregss);
                  
     } else if (CCTK_EQUALS(refined_regions, "centre")) {
       
-      do_recompose = Centre (cctkGH, hh, regsss);
+      do_recompose = Centre (cctkGH, hh, superregss);
                  
     } else if (CCTK_EQUALS(refined_regions, "manual-gridpoints")) {
       
       do_recompose
-        = ManualGridpoints (cctkGH, hh, regsss);
+        = ManualGridpoints (cctkGH, hh, superregss);
                  
     } else if (CCTK_EQUALS(refined_regions, "manual-coordinates")) {
       
       do_recompose
-        = ManualCoordinates (cctkGH, hh, regsss);
+        = ManualCoordinates (cctkGH, hh, superregss);
                  
     } else if (CCTK_EQUALS(refined_regions, "manual-gridpoint-list")) {
       
       do_recompose
-        = ManualGridpointList (cctkGH, hh, regsss);
+        = ManualGridpointList (cctkGH, hh, superregss);
                  
     } else if (CCTK_EQUALS(refined_regions, "manual-coordinate-list")) {
       
       do_recompose
-        = ManualCoordinateList (cctkGH, hh, regsss);
+        = ManualCoordinateList (cctkGH, hh, superregss);
                  
     } else if (CCTK_EQUALS(refined_regions, "moving")) {
       
-      do_recompose = Moving (cctkGH, hh, regsss);
+      do_recompose = Moving (cctkGH, hh, superregss);
                  
     } else if (CCTK_EQUALS(refined_regions, "automatic")) {
       
-      do_recompose = Automatic (cctkGH, hh, regsss);
+      do_recompose = Automatic (cctkGH, hh, superregss);
                  
     } else {
       assert (0);
     }
     
+    // make multiprocessor aware
+    vector<vector<region_t> > regss(superregss.size());
+    for (size_t rl=0; rl<superregss.size(); ++rl) {
+#warning "TODO: delete .processors"
+      SplitRegions (cctkGH, superregss.at(rl), regss.at(rl));
+    }
+    
+    // make multigrid aware
+    MakeMultigridBoxes (cctkGH, Carpet::map, regss, regsss);
+    
     assert (regsss.size() > 0);
     
     return do_recompose;
diff --git a/Carpet/CarpetRegrid/src/regrid.hh b/Carpet/CarpetRegrid/src/regrid.hh
index b7d727317..f37c9cb1d 100644
--- a/Carpet/CarpetRegrid/src/regrid.hh
+++ b/Carpet/CarpetRegrid/src/regrid.hh
@@ -31,8 +31,10 @@ namespace CarpetRegrid {
     
     // Aliased functions
 //     CCTK_INT CarpetRegrid_Regrid (const cGH * const cctkGH,
+//                                   gh<dim>::rregs * superregss,
 //                                   gh<dim>::mregs * regsss);
     CCTK_INT CarpetRegrid_Regrid (CCTK_POINTER_TO_CONST const cctkGH_,
+                                  CCTK_POINTER const superregss_,
                                   CCTK_POINTER const regsss_,
 				  CCTK_INT force);
   }
@@ -41,15 +43,15 @@ namespace CarpetRegrid {
   
   int BaseLevel (cGH const * const cctkGH,
                  gh const & hh,
-                 gh::mregs & regsss);
+                 gh::rregs & regss);
 
   int Centre (cGH const * const cctkGH,
               gh const & hh,
-              gh::mregs & regsss);
+              gh::rregs & regss);
 
   int ManualGridpoints (cGH const * const cctkGH,
                         gh const & hh,
-                        gh::mregs & regsss);
+                        gh::rregs & regss);
   
   void ManualGridpoints_OneLevel (const cGH * const cctkGH,
                                   const gh & hh,
@@ -62,7 +64,7 @@ namespace CarpetRegrid {
 
   int ManualCoordinates (cGH const * const cctkGH,
                          gh const & hh,
-                         gh::mregs & regsss);
+                         gh::rregs & regss);
   
   void ManualCoordinates_OneLevel (const cGH * const cctkGH,
                                    const gh & hh,
@@ -84,19 +86,19 @@ namespace CarpetRegrid {
 
   int ManualGridpointList (cGH const * const cctkGH,
                            gh const & hh,
-                           gh::mregs & regsss);
+                           gh::rregs & regss);
   
   int ManualCoordinateList (cGH const * const cctkGH,
                             gh const & hh,
-                            gh::mregs & regsss);
+                            gh::rregs & regss);
 
   int Moving (cGH const * const cctkGH,
               gh const & hh,
-              gh::mregs & regsss);
+              gh::rregs & regss);
   
   int Automatic (cGH const * const cctkGH,
                  gh const & hh,
-                 gh::mregs & regsss);
+                 gh::rregs & regss);
   
   void Automatic_OneLevel (const cGH * const cctkGH,
                            const gh & hh,
diff --git a/Carpet/CarpetRegrid2/interface.ccl b/Carpet/CarpetRegrid2/interface.ccl
index aa6f799da..6eb913280 100644
--- a/Carpet/CarpetRegrid2/interface.ccl
+++ b/Carpet/CarpetRegrid2/interface.ccl
@@ -87,21 +87,25 @@ USES FUNCTION MultiPatch_ConvertFromPhysicalBoundary
 
 # The true prototype of the routine below:
 #   int Carpet_Regrid (cGH const * cctkGH,
+#                      gh::rregs * superregss,
 #                      gh::mregs * regsss,
 #                      int         force);
-CCTK_INT FUNCTION Carpet_Regrid       \
-    (CCTK_POINTER_TO_CONST IN cctkGH, \
-     CCTK_POINTER          IN regsss, \
+CCTK_INT FUNCTION Carpet_Regrid           \
+    (CCTK_POINTER_TO_CONST IN cctkGH,     \
+     CCTK_POINTER          IN superregss, \
+     CCTK_POINTER          IN regsss,     \
      CCTK_INT              IN force)
 PROVIDES FUNCTION Carpet_Regrid WITH CarpetRegrid2_Regrid LANGUAGE C
 
 # The true prototype of the routine below:
 #   int Carpet_Regrid (cGH const          * cctkGH,
+#                      vector <gh::rregs> * superregsss,
 #                      vector <gh::mregs> * regssss,
 #                      int                  force)
-CCTK_INT FUNCTION Carpet_RegridMaps    \
-    (CCTK_POINTER_TO_CONST IN cctkGH,  \
-     CCTK_POINTER          IN regssss, \
+CCTK_INT FUNCTION Carpet_RegridMaps        \
+    (CCTK_POINTER_TO_CONST IN cctkGH,      \
+     CCTK_POINTER          IN superregsss, \
+     CCTK_POINTER          IN regssss,     \
      CCTK_INT              IN force)
 PROVIDES FUNCTION Carpet_RegridMaps WITH CarpetRegrid2_RegridMaps LANGUAGE C
 
diff --git a/Carpet/CarpetRegrid2/src/regrid.cc b/Carpet/CarpetRegrid2/src/regrid.cc
index c08705d29..eaf993bc0 100644
--- a/Carpet/CarpetRegrid2/src/regrid.cc
+++ b/Carpet/CarpetRegrid2/src/regrid.cc
@@ -239,11 +239,13 @@ namespace CarpetRegrid2 {
   extern "C" {
     CCTK_INT
     CarpetRegrid2_Regrid (CCTK_POINTER_TO_CONST const cctkGH_,
+                          CCTK_POINTER          const superregss_,
                           CCTK_POINTER          const regsss_,
                           CCTK_INT              const force);
     
     CCTK_INT
     CarpetRegrid2_RegridMaps (CCTK_POINTER_TO_CONST const cctkGH_,
+                              CCTK_POINTER          const superregsss_,
                               CCTK_POINTER          const regssss_,
                               CCTK_INT              const force);
   }
@@ -759,6 +761,7 @@ namespace CarpetRegrid2 {
   
   CCTK_INT
   CarpetRegrid2_Regrid (CCTK_POINTER_TO_CONST const cctkGH_,
+                        CCTK_POINTER          const superregss_,
                         CCTK_POINTER          const regsss_,
                         CCTK_INT              const force)
   {
@@ -882,16 +885,15 @@ namespace CarpetRegrid2 {
     
     if (do_recompose) {
       
-      gh::mregs & regsss = * static_cast <gh::mregs *> (regsss_);
+      gh::rregs & superregss = * static_cast <gh::rregs *> (superregss_);
+      gh::mregs & regsss     = * static_cast <gh::mregs *> (regsss_);
       
-      // Make multigrid unaware
-      vector <vector <region_t> > regss = regsss.at(0);
-      
-      Regrid (cctkGH, regss);
+      Regrid (cctkGH, superregss);
       
       // Make multiprocessor aware
+      vector <vector <region_t> > regss;
       for (size_t rl = 0; rl < regss.size(); ++ rl) {
-        Carpet::SplitRegions (cctkGH, regss.at(rl));
+        Carpet::SplitRegions (cctkGH, superregss.at(rl), regss.at(rl));
       } // for rl
       
       // Make multigrid aware
@@ -919,6 +921,7 @@ namespace CarpetRegrid2 {
   
   CCTK_INT
   CarpetRegrid2_RegridMaps (CCTK_POINTER_TO_CONST const cctkGH_,
+                            CCTK_POINTER          const superregsss_,
                             CCTK_POINTER          const regssss_,
                             CCTK_INT              const force)
   {
@@ -1040,23 +1043,21 @@ namespace CarpetRegrid2 {
     
     if (do_recompose) {
       
+      vector <gh::rregs> & superregsss =
+        * static_cast <vector <gh::rregs>  *> (superregsss_);
       vector <gh::mregs> & regssss =
         * static_cast <vector <gh::mregs>  *> (regssss_);
       
-      // Make multigrid unaware
-      vector< vector <vector <region_t> > > regsss (Carpet::maps);
-      for (int m = 0; m < maps; ++ m) {
-        regsss.at(m) = regssss.at(m).at(0);
-      }
-      
       BEGIN_MAP_LOOP (cctkGH, CCTK_GF) {
-        Regrid (cctkGH, regsss.at(Carpet::map));
+        Regrid (cctkGH, superregsss.at(Carpet::map));
       } END_MAP_LOOP;
       
+      vector< vector <vector <region_t> > > regsss (Carpet::maps);
+      
       // Count levels
       vector <int> rls (maps);
       for (int m = 0; m < maps; ++ m) {
-        rls.at(m) = regsss.at(m).size();
+        rls.at(m) = superregsss.at(m).size();
       }
       int maxrl = 0;
       for (int m = 0; m < maps; ++ m) {
@@ -1064,17 +1065,21 @@ namespace CarpetRegrid2 {
       }
       // All maps must have the same number of levels
       for (int m = 0; m < maps; ++ m) {
+        superregsss.at(m).resize (maxrl);
         regsss.at(m).resize (maxrl);
       }
       
       // Make multiprocessor aware
       for (int rl = 0; rl < maxrl; ++ rl) {
-        vector <vector <region_t> > regss (maps);
+        vector <vector <region_t> > superregss (maps);
         for (int m = 0; m < maps; ++ m) {
-          regss.at(m) = regsss.at(m).at(rl);
+          superregss.at(m) = superregsss.at(m).at(rl);
         }
-        Carpet::SplitRegionsMaps (cctkGH, regss);
+        vector <vector <region_t> > regss (maps);
+        Carpet::SplitRegionsMaps (cctkGH, superregss, regss);
+        
         for (int m = 0; m < maps; ++ m) {
+          superregsss.at(m).at(rl) = superregss.at(m);
           regsss.at(m).at(rl) = regss.at(m);
         }
       } // for rl
diff --git a/Carpet/CarpetWeb/get-carpet.html b/Carpet/CarpetWeb/get-carpet.html
index 00326b11b..54d89af72 100644
--- a/Carpet/CarpetWeb/get-carpet.html
+++ b/Carpet/CarpetWeb/get-carpet.html
@@ -253,7 +253,7 @@ darcs pull</pre>
     <p>The <a href="http://git.or.cz/">git web site</a> contains
     introductions and documentation for git.  The Linux kernel
     developers also maintain
-    a <a href="http://www.kernel.org/pub/software/scm/git/docs/tutorial.html">tutorial</a> for
+    a <a href="http://www.kernel.org/pub/software/scm/git/docs/gittutorial.html">tutorial</a> for
     git.  Git should be available for all modern operating systems.
     It is also not difficult to install manually.</p>
 
diff --git a/Carpet/CarpetWeb/publications.html b/Carpet/CarpetWeb/publications.html
index aacf1bb04..c67ec709d 100644
--- a/Carpet/CarpetWeb/publications.html
+++ b/Carpet/CarpetWeb/publications.html
@@ -44,7 +44,8 @@ Springer, Berlin 2003
 <p>Erik Schnetter, Scott H. Hawley, Ian Hawke,<br />
 <i>Evolutions in 3D numerical relativity using fixed mesh
   refinement,</i><br />
-<a href="http://www.iop.org/EJ/abstract/0264-9381/21/6/014">Class. Quantum Grav. <b>21</b>, 1465-1488 (2004)</a>,<br />
+<a href="http://dx.doi.org/10.1088/0264-9381/21/6/014">Class. Quantum
+  Grav. <b>21</b>, 1465-1488 (2004)</a>,<br />
 <a href="http://arxiv.org/abs/gr-qc/0310042">arXiv:gr-qc/0310042</a>.</p>
 </li>
 
@@ -83,7 +84,8 @@ Springer, Berlin 2003
 <p>Erik Schnetter, Scott H. Hawley, Ian Hawke,<br />
 <i>Evolutions in 3D numerical relativity using fixed mesh
   refinement,</i><br />
-<a href="http://www.iop.org/EJ/abstract/0264-9381/21/6/014">Class. Quantum Grav. <b>21</b>, 1465-1488 (2004)</a>,<br />
+<a href="http://dx.doi.org/10.1088/0264-9381/21/6/014">Class. Quantum
+  Grav. <b>21</b>, 1465-1488 (2004)</a>,<br />
 <a href="http://arxiv.org/abs/gr-qc/0310042">arXiv:gr-qc/0310042</a>.</p>
 </li>
 
@@ -137,7 +139,7 @@ Seidel, Ryoji Takahashi, Jonathan Thornburg, Jason Ventrella,<br />
   satisfying summation by parts, and applications in three-dimensional
   multi-block evolutions,</i><br />
 <a
-   href="http://www.springerlink.com/content/l724hr0846n2/">J. Sci. Comp. <b>32</b>,
+   href="http://www.springerlink.com/content/l724hr0846n2/">J. Sci. Comput. <b>32</b>,
    109-145 (2007)</a>,<br />
 <a href="http://arxiv.org/abs/gr-qc/0512001">arXiv:gr-qc/0512001</a>.</p>
 </li>
@@ -148,7 +150,7 @@ Seidel, Ryoji Takahashi, Jonathan Thornburg, Jason Ventrella,<br />
 <p>Erik Schnetter, Peter Diener, Ernst Nils Dorband, Manuel Tiglio,<br />
 <i>A multi-block infrastructure for three-dimensional time-dependent
   numerical relativity</i>,<br />
-<a href="http://www.iop.org/EJ/abstract/0264-9381/23/16/S14/">Class. Quantum
+<a href="http://dx.doi.org/10.1088/0264-9381/23/16/S14">Class. Quantum
   Grav. <b>23</b>, S553-S578 (2006)</a>,<br />
 <a href="http://arxiv.org/abs/gr-qc/0602104">arXiv:gr-qc/0602104</a>.</p>
 </li>
@@ -156,10 +158,10 @@ Seidel, Ryoji Takahashi, Jonathan Thornburg, Jason Ventrella,<br />
 <li>
 <!-- Sopuerta:2006bw -->
 <!-- received 2006-03-20 -->
-<p>Carlos F. Sopuerta and Ulrich Sperhake and Pablo Laguna,<br />
+<p>Carlos F. Sopuerta, Ulrich Sperhake, Pablo Laguna,<br />
 <i>Hydro-without-Hydro Framework for Simulations of Black Hole-Neutron
   Star Binaries</i>,<br />
-<a href="http://www.iop.org/EJ/abstract/0264-9381/23/16/S15/">Class. Quantum
+<a href="http://dx.doi.org/10.1088/0264-9381/23/16/S15">Class. Quantum
   Grav. <b>23</b>, S579-S598 (2006)</a>,<br />
 <a href="http://arxiv.org/abs/gr-qc/0605018">arXiv:gr-qc/0605018</a>.</p>
 </li>
@@ -229,8 +231,8 @@ Janka, Ian Hawke, Burkhard Zink, Erik Schnetter,<br />
   Thornburg, Béla Szilágyi,<br />
 <i>Characteristic evolutions in numerical relativity using six angular
   patches,</i><br />
-<a
-   href="http://www.iop.org/EJ/abstract/0264-9381/24/12/S21/">Class. Quantum Grav. <b>24</b>, S327-S339 (2007)</a>,<br />
+<a href="http://dx.doi.org/10.1088/0264-9381/24/12/S21">Class. Quantum
+  Grav. <b>24</b>, S327-S339 (2007)</a>,<br />
 <a href="http://arxiv.org/abs/gr-qc/0610019">arXiv:gr-qc/0610019</a>.</p>
 </li>
 
@@ -260,8 +262,8 @@ Erik Schnetter, Ewald Müller,<br />
 <p>Christian D. Ott, Harald Dimmelmeier, Andreas Marek, Hans-Thomas
 Janka, Burkhard Zink, Ian Hawke, Erik Schnetter,<br />
 <i>Rotating collapse of stellar iron cores in general relativity,</i><br />
-<a
-   href="http://www.iop.org/EJ/abstract/0264-9381/24/12/S10/">Class. Quantum Grav. <b>24</b> S139-S154 (2007)</a>,<br />
+<a href="http://dx.doi.org/10.1088/0264-9381/24/12/S10">Class. Quantum
+  Grav. <b>24</b> S139-S154 (2007)</a>,<br />
 <a href="http://arxiv.org/abs/astro-ph/0612638">arXiv:astro-ph/0612638</a>.</p>
 </li>
 
@@ -270,8 +272,8 @@ Janka, Burkhard Zink, Ian Hawke, Erik Schnetter,<br />
 <p>Luca Baiotti, Ian Hawke, Luciano Rezzolla,<br />
 <i>On the gravitational radiation from the collapse of neutron stars
   to rotating black holes,</i><br />
-<a
-   href="http://www.iop.org/EJ/abstract/0264-9381/24/12/S13/">Class. Quantum Grav. <b>24</b> S187-S206 (2007)</a>,<br />
+<a href="http://dx.doi.org/10.1088/0264-9381/24/12/S13">Class. Quantum
+  Grav. <b>24</b> S187-S206 (2007)</a>,<br />
 <a href="http://arxiv.org/abs/gr-qc/0701043">arXiv:gr-qc/0701043</a>.</p>
 </li>
 
@@ -281,8 +283,8 @@ Janka, Burkhard Zink, Ian Hawke, Erik Schnetter,<br />
 Palenzuela, Erik Schnetter, Manuel Tiglio,<br />
 <i>How far away is far enough for extracting numerical waveforms, and
   how much do they depend on the extraction method?</i><br />
-<a
-   href="http://www.iop.org/EJ/abstract/0264-9381/24/12/S22/">Class. Quantum Grav. <b>24</b>, S341-S368 (2007)</a>,<br />
+<a href="http://dx.doi.org/10.1088/0264-9381/24/12/S22">Class. Quantum
+  Grav. <b>24</b>, S341-S368 (2007)</a>,<br />
 <a href="http://arxiv.org/abs/gr-qc/0612149">arXiv:gr-qc/0612149</a>.</p>
 </li>
 
@@ -292,8 +294,8 @@ Palenzuela, Erik Schnetter, Manuel Tiglio,<br />
   Jeffrey Winicour,<br />
 <i>An explicit harmonic code for black-hole evolution using
   excision,</i><br />
-<a
-   href="http://www.iop.org/EJ/abstract/0264-9381/24/12/S18/">Class. Quantum Grav. <b>24</b>, S275-S293 (2007)</a>,<br />
+<a href="http://dx.doi.org/10.1088/0264-9381/24/12/S18">Class. Quantum
+  Grav. <b>24</b>, S275-S293 (2007)</a>,<br />
 <a href="http://arxiv.org/abs/gr-qc/0612150">arXiv:gr-qc/0612150</a>.</p>
 </li>
 
@@ -303,7 +305,7 @@ Palenzuela, Erik Schnetter, Manuel Tiglio,<br />
 <p>Frank Herrmann, Ian Hinder, Deirdre M. Shoemaker, Pablo Laguna,<br />
 <i>Unequal mass binary black hole plunges and gravitational
   recoil,</i><br />
-<a href="http://www.iop.org/EJ/abstract/0264-9381/24/12/S04/">Class. Quantum
+<a href="http://dx.doi.org/10.1088/0264-9381/24/12/S04">Class. Quantum
   Grav. <b>24</b>, S33-S42 (2007)</a>,<br />
 <a href="http://arxiv.org/abs/gr-qc/0601026">arXiv:gr-qc/0601026</a>.</p>
 </li>
@@ -313,7 +315,7 @@ Palenzuela, Erik Schnetter, Manuel Tiglio,<br />
 <p>Bruno Giacomazzo, Luciano Rezzolla,<br />
 <i>WhiskyMHD: a new numerical code for general relativistic
   magnetohydrodynamics,</i><br />
-<a href="http://www.iop.org/EJ/abstract/0264-9381/24/12/S16/">Class. Quantum
+<a href="http://dx.doi.org/10.1088/0264-9381/24/12/S16">Class. Quantum
   Grav. <b>24</b>, S235-S258 (2007)</a>,<br />
 <a href="http://arxiv.org/abs/gr-qc/0701109">arXiv:gr-qc/0701109</a>.</p>
 </li>
@@ -323,7 +325,7 @@ Palenzuela, Erik Schnetter, Manuel Tiglio,<br />
 <p>Pedro Marronetti, Wolfgang Tichy, Bernd Brügmann, José González,
   Mark Hannam, Sascha Husa, Ulrich Sperhake,<br />
 <i>Binary black holes on a budget: simulations using workstations,</i><br />
-<a href="http://www.iop.org/EJ/abstract/0264-9381/24/12/S05/">Class. Quantum
+<a href="http://dx.doi.org/10.1088/0264-9381/24/12/S05">Class. Quantum
   Grav. <b>24</b>, S45-S58 (2007)</a>,<br />
 <a href="http://arxiv.org/abs/gr-qc/0701123">arXiv:gr-qc/0701123</a>.</p>
 </li>
@@ -355,7 +357,7 @@ Palenzuela, Erik Schnetter, Manuel Tiglio,<br />
 <p>Manuela Campanelli, Carlos O. Lousto, Yosef Zlochower, David
   Merritt,<br />
 <i>Large Merger Recoils and Spin Flips From Generic Black-Hole Binaries,</i><br />
-<a href="http://www.journals.uchicago.edu/doi/abs/10.1086/516712">Astrophys. J. <b>659</b>,
+<a href="http://www.journals.uchicago.edu/doi/abs/10.1086/516712">Astrophys. J. Lett.<b>659</b>,
   L5-L8 (2007)</a>,<br />
 <a href="http://arxiv.org/abs/gr-qc/0701164">arXiv:gr-qc/0701164</a>.</p>
 </li>
@@ -388,8 +390,8 @@ Palenzuela, Erik Schnetter, Manuel Tiglio,<br />
 <p>Jonathan Thornburg, Peter Diener, Denis Pollney, Luciano Rezzolla,
   Erik Schnetter, Ed Seidel, Ryoji Takahashi,<br />
 <i>Are moving punctures equivalent to moving black holes?,</i><br />
-<a
-   href="http://www.iop.org/EJ/abstract/0264-9381/24/15/009/">Class. Quantum Grav. <b>24</b>, 3911-3918 (2007)</a>,<br />
+<a href="http://dx.doi.org/10.1088/0264-9381/24/15/009">Class. Quantum
+  Grav. <b>24</b>, 3911-3918 (2007)</a>,<br />
 <a href="http://arxiv.org/abs/gr-qc/0701038">arXiv:gr-qc/0701038</a>.</p>
 </li>
 
@@ -402,8 +404,7 @@ Palenzuela, Erik Schnetter, Manuel Tiglio,<br />
  Sperhake, Jonathan Thornburg,<br />
 <i>Phenomenological template family for black-hole coalescence
   waveforms,</i><br />
-<a
-   href="http://www.iop.org/EJ/abstract/0264-9381/24/19/S31">Class. Quantum
+<a href="http://dx.doi.org/10.1088/0264-9381/24/19/S31">Class. Quantum
    Grav. <b>24</b> S689-S699 (2007)</a>,<br />
 <a href="http://arxiv.org/abs/0704.3764">arXiv:0704.3764
   [gr-qc]</a>.</p>
@@ -422,6 +423,16 @@ Palenzuela, Erik Schnetter, Manuel Tiglio,<br />
 </li>
 
 <li>
+<!-- received 2007-06-18 -->
+<p>Frank Herrmann, Ian Hinder, Deirdre M. Shoemaker, Pablo Laguna,
+  Richard A. Matzner,<br />
+<i>Binary Black Holes: Spin Dynamics and Gravitational Recoil,</i><br />
+<a href="http://link.aps.org/abstract/PRD/v76/e084032">Phys. Rev. D <b>76</b>, 084032 (2007)</a>,<br />
+<a href="http://arxiv.org/abs/0706.2541">arXiv:0706.2541
+  [gr-qc]</a>.</p>
+</li>
+
+<li>
 <!-- received 2007-07-06 -->
 <p>Badri Krishnan, Carlos O. Lousto, Yosef Zlochower,<br />
 <i>Quasi-Local Linear Momentum in Black-Hole Binaries,</i><br />
@@ -445,6 +456,19 @@ Palenzuela, Erik Schnetter, Manuel Tiglio,<br />
 </li>
 
 <li>
+<!-- Rezzolla2007a -->
+<!-- received 2007-08-29 -->
+<p>Luciano Rezzolla, Ernst Nils Dorband, Christian Reisswig, Peter
+  Diener, Denis Pollney, Erik Schnetter, Béla Szilágyi,<br />
+<i>Spin Diagrams for Equal-Mass Black-Hole Binaries with Aligned
+  Spins,</i><br />
+<a href="http://www.journals.uchicago.edu/doi/abs/10.1086/587679">Astrophys. J. <b>679</b>,
+  1422-1426 (2008)</a>,<br />
+<a href="http://arxiv.org/abs/0708.3999">arXiv:0708.3999
+  [gr-qc]</a>.</p>
+</li>
+
+<li>
 <!-- received 2007-08-30 -->
 <p>Carlos O. Lousto, Yosef Zlochower,<br />
 <i>Further insight into gravitational recoil,</i><br />
@@ -456,50 +480,20 @@ Palenzuela, Erik Schnetter, Manuel Tiglio,<br />
 </li>
 
 <li>
-<p>Frank Herrmann, Ian Hinder, Deirdre M. Shoemaker, Pablo Laguna,
-  Richard A. Matzner,<br />
-<i>Binary Black Holes: Spin Dynamics and Gravitational Recoil,</i><br />
-<a href="http://prd.aps.org/">Phys. Rev. D (accepted for
-  publication)</a>,<br />
-<a href="http://arxiv.org/abs/0706.2541">arXiv:0706.2541
-  [gr-qc]</a>.</p>
-</li>
-
-<li>
+<!-- received 2007-09-06 -->
 <p>Denis Pollney, Christian Reisswig, Luciano Rezzolla, Bela Szilagyi,
   Marcus Ansorg, Barrett Deris, Peter Diener, Ernst Nils Dorband,
   Michael Koppitz, Alessandro Nagar, Erik Schnetter,<br />
 <i>Recoil velocities from equal-mass binary black-hole mergers: a
   systematic investigation of spin-orbit aligned configurations,</i><br />
-<a href="http://prd.aps.org/">Phys. Rev. D (accepted for
-  publication)</a>,<br />
+<a href="http://link.aps.org/abstract/PRD/v76/e124002">Phys. Rev. D <b>76</b>,
+  124002 (2007)</a>,<br />
 <a href="http://arxiv.org/abs/0707.2559">arXiv:0707.2559
   [gr-qc]</a>.</p>
 </li>
 
 <li>
-<!-- Rezzolla2007a -->
-<p>Luciano Rezzolla, Ernst Nils Dorband, Christian Reisswig, Peter
-  Diener, Denis Pollney, Erik Schnetter, Béla Szilágyi,<br />
-<i>Spin Diagrams for Equal-Mass Black-Hole Binaries with Aligned
-  Spins,</i><br />
-<a href="http://www.journals.uchicago.edu/toc/apj/current">Astrophys. J. (accepted
-  for publication)</a>,<br />
-<a href="http://arxiv.org/abs/0708.3999">arXiv:0708.3999
-  [gr-qc]</a>.</p>
-</li>
-
-<li>
-<!-- uses Carpet only indirectly -->
-<p>Latham Boyle, Michael Kesden, Samaya Nissanke,<br />
-<i>Binary black hole merger: symmetry and the spin expansion,</i><br />
-<a href="http://prl.aps.org/">Phys. Rev. Lett. (accepted for
-  publication)</a>,<br />
-<a href="http://arxiv.org/abs/0709.0299">arXiv:0709.0299
-  [gr-qc]</a>.</p>
-</li>
-
-<li>
+<!-- received 2007-10-15 -->
 <p>Parameswaran Ajith, Stanislav Babak, Yanbei Chen, Martin Hewitson,
   Badri Krishnan, Alicia M. Sintes, John T. Whelan, Bernd Brügmann,
   Peter Diener, Ernst Nils Dorband, José González, Mark Hannam, Sascha Husa,
@@ -507,33 +501,68 @@ Palenzuela, Erik Schnetter, Manuel Tiglio,<br />
   Jonathan Thornburg,<br />
 <i>A template bank for gravitational waveforms from coalescing binary
   black holes: I. non-spinning binaries,</i><br />
-<a href="http://prd.aps.org/">Phys. Rev. D (accepted for
-  publication)</a>,<br />
+<a href="http://link.aps.org/abstract/PRD/v77/e104017">Phys. Rev. D <b>77</b>,
+  104017 (2008)</a>,<br />
 <a href="http://arxiv.org/abs/0710.2335">arXiv:0710.2335
   [gr-qc]</a>.</p>
 </li>
 
 <li>
+<!-- received 2007-10-18 -->
 <p>Luciano Rezzolla, Peter Diener, Ernst Nils Dorband, Denis Pollney,
   Christian Reisswig, Erik Schnetter, Jennifer Seiler,<br />
 <i>The final spin from the coalescence of aligned-spin black-hole
   binaries,</i><br />
-<a href="http://www.journals.uchicago.edu/toc/apjl/current">ApJ
-  letters (accepted for publication)</a>,<br />
+<a href="http://www.journals.uchicago.edu/doi/abs/10.1086/528935">Astrophys. J. Lett.<b>674</b>,
+  L29-L32 (2008)</a>,<br />
 <a href="http://arxiv.org/abs/0710.3345">arXiv:0710.3345
   [gr-qc]</a>.</p>
 </li>
 
 <li>
+<!-- received 2007-10-22 -->
+<!-- uses Carpet only indirectly -->
+<p>Latham Boyle, Michael Kesden, Samaya Nissanke,<br />
+<i>Binary black hole merger: symmetry and the spin expansion,</i><br />
+<a href="http://link.aps.org/abstract/PRL/v100/e151101">Phys. Rev. Lett. <b>100</b>, 151101 (2008)</a>,<br />
+<a href="http://arxiv.org/abs/0709.0299">arXiv:0709.0299
+  [gr-qc]</a>.</p>
+</li>
+
+<li>
+<!-- received 2007-11-07 -->
+<p>Emanuele Berti, Vitor Cardoso, José A. González, Ulrich Sperhake,
+  Bernd Brügmann,<br />
+<i>Multipolar analysis of spinning binaries,</i><br />
+<a href="http://dx.doi.org/10.1088/0264-9381/25/11/114035">Class. Quantum
+  Grav. <b>25</b>, 114035 (2008)</a>,<br />
+<a href="http://arxiv.org/abs/0711.1097">arXiv:0711.1097
+  [gr-qc]</a>.</p>
+</li>
+
+<li>
+<!-- received 2007-12-04 -->
 <p>Burkhard Zink, Erik Schnetter, Manuel Tiglio,<br />
 <i>Multi-patch methods in general relativistic astrophysics -
   I. Hydrodynamical flows on fixed backgrounds,</i><br />
-<a href="http://prd.aps.org/">Phys. Rev. D (accepted for
-  publication)</a>,<br />
+<a href="http://link.aps.org/abstract/PRD/v77/e103015">Phys. Rev. D <b>77</b>,
+  103015 (2008)</a>,<br />
 <a href="http://arxiv.org/abs/0712.0353">arXiv:0712.0353
   [astro-ph]</a>.</p>
 </li>
 
+<li>
+<!-- received 2007-12-14 -->
+<p>Deirdre M. Shoemaker, Birjoo Vaishnav, Ian Hinder, Frank
+  Herrmann,<br />
+<i>Numerical relativity meets data analysis: spinning binary black
+  hole case,</i><br />
+<a href="http://dx.doi.org/10.1088/0264-9381/25/11/114047">Class. Quantum
+  Grav. <b>25</b>, 114047 (2008)</a>,<br />
+<a href="http://arxiv.org/abs/0802.4427">arXiv:0802.4427
+  [gr-qc]</a>.</p>
+</li>
+
 </ol>
 
 
@@ -546,7 +575,7 @@ Palenzuela, Erik Schnetter, Manuel Tiglio,<br />
 
 <li>
 <p>Burkhard Zink, Nikolaos Stergioulas, Ian Hawke, Christian D. Ott,
-  Erik Schnetter, and Ewald Müller,<br />
+  Erik Schnetter, Ewald Müller,<br />
 <i>Rotational instabilities in supermassive stars: a new way to form
   supermassive black holes,</i>,<br /> in N. K. Spyrou,
   N. Stergioulas, and C. Tsagas, editors, <i>International Scientific
@@ -588,16 +617,14 @@ Palenzuela, Erik Schnetter, Manuel Tiglio,<br />
 <p>Luca Baiotti, Ian Hawke, Luciano Rezzolla, Erik Schnetter,
 <i>Details on the gravitational-wave emission from rotating
   gravitational collapse in 3D,</i><br />
-<a
-   href="http://www.iop.org/EJ/abstract/1742-6596/66/1/012045/">J. Phys.:
-   Conf. Ser. <b>66</b>, 012045 (2007)</a>.</p>
+<a href="http://dx.doi.org/10.1088/1742-6596/66/1/012045">J. Phys.:
+  Conf. Ser. <b>66</b>, 012045 (2007)</a>.</p>
 </li>
 
 <li>
 <p>Ulrich Sperhake,<br />
 <i>Black-hole binary evolutions with the LEAN code,</i><br />
-<a
-   href="http://www.iop.org/EJ/abstract/1742-6596/66/1/012049/">J. Phys.:
+<a href="http://dx.doi.org/10.1088/1742-6596/66/1/012049">J. Phys.:
    Conf. Ser. <b>66</b> 012049 (2007)</a>.</p>
 </li>
 
@@ -605,8 +632,7 @@ Palenzuela, Erik Schnetter, Manuel Tiglio,<br />
 <!-- uses Carpet only indirectly -->
 <p>José A. Font,<br />
 <i>Current status of relativistic core collapse simulations,</i><br />
-<a
-   href="http://www.iop.org/EJ/abstract/1742-6596/66/1/012063/">J. Phys.:
+<a href="http://dx.doi.org/10.1088/1742-6596/66/1/012063">J. Phys.:
    Conf. Ser. <b>66</b>, 012063 (2007)</a>.</p>
 </li>
 
@@ -615,8 +641,7 @@ Palenzuela, Erik Schnetter, Manuel Tiglio,<br />
   Erik Schnetter, Ewald Müller,<br />
 <i>Supermassive Black Hole Formation through Rotational
   Instabilities,</i><br />
-<a
-   href="http://www.iop.org/EJ/abstract/1742-6596/68/1/012050/">J. Phys.:
+<a href="http://dx.doi.org/10.1088/1742-6596/68/1/012050">J. Phys.:
    Conf. Ser. <b>68</b>, 012050 (2007)</a>.</p>
 </li>
 
@@ -740,15 +765,6 @@ triangulations using finite element methods,</i><br />
 </li>
 
 <li>
-<p>Deirdre M. Shoemaker, Birjoo Vaishnav, Ian Hinder, Frank
-  Herrmann,<br />
-<i>Numerical Relativity meets Data Analysis: Spinning Binary Black
-  Hole Case,</i><br />
-<a href="http://arxiv.org/abs/0802.4427">arXiv:0802.4427
-  [gr-qc]</a>.</p>
-</li>
-
-<li>
 <p>Sergio Dain, Carlos O. Lousto, Yosef Zlochower,<br />
 <i>Extra-Large Remnant Recoil Velocities and Spins from
   Near-Extremal-Bowen-York-Spin Black-Hole Binaries,</i><br />
@@ -756,6 +772,47 @@ triangulations using finite element methods,</i><br />
   [gr-qc]</a>.</p>
 </li>
 
+<li>
+<p>Accurate evolutions of inspiralling neutron-star binaries: prompt
+  and delayed collapse to black hole,<br />
+<i>Luca Baiotti, Bruno Giacomazzo, Luciano Rezzolla,</i><br />
+<a href="http://arxiv.org/abs/0804.0594">arXiv:0804.0594
+  [gr-qc]</a>.</p>
+</li>
+
+<li>
+  <p>Modeling gravitational recoil from precessing highly-spinning
+  unequal-mass black-hole binaries,<br />
+<i>Carlos O. Lousto, Yosef Zlochower,</i><br />
+<a href="http://arxiv.org/abs/0805.0159">arXiv:0805.0159
+  [gr-qc]</a>.</p>
+</li>
+
+<li>
+  <p>Transformation of the multipolar components of gravitational
+    radiation under rotations and boosts,<br />
+    <i>Leonardo Gualtieri, Emanuele Berti, Vitor Cardoso, Ulrich
+      Sperhake,</i><br />
+    <a href="http://arxiv.org/abs/0805.1017">arXiv:0805.1017
+      [gr-qc]</a>.</p>
+</li>
+
+<li>
+  <p>Comparisons of eccentric binary black hole simulations with
+    post-Newtonian models,<br />
+    <i>Ian Hinder, Frank Herrmann, Pablo Lagona, Deirdre Shoemaker,</i><br />
+    <a href="http://arxiv.org/abs/0806.1037">arXiv:0806.1037
+      [gr-qc]</a>.</p>
+</li>
+
+<li>
+  <p>The high-energy collision of two black holes,<br />
+    <i>Ulrich Sperhake, Vitor Cardoso, Frans Pretorius, Emanuele
+    Berti, José A. González,</i><br />
+    <a href="http://arxiv.org/abs/0806.1738">arXiv:0806.1738
+      [gr-qc]</a>.</p>
+</li>
+
 </ol>
 
 
@@ -834,6 +891,13 @@ PhD thesis, Max-Planck Institute for Gravitational Physics (AEI) and
   [gr-qc]</a>.</p>
 </li>
 
+<li>
+<p>Michael Jasiulek,<br />
+Spin Measures on Isolated and Dynamical Horizons in Numerical
+  Relativity,<br />
+Diplomarbeit, Humboldt-Universität zu Berlin, 2008.</p>
+</li>
+
 </ol>
 
 
@@ -895,7 +959,7 @@ PhD thesis, SISSA, 2006.
 <p>
 <!-- Created: Sun Feb 26 2006 -->
 <!-- hhmts start -->
-Last modified: Mon Mar 31 2008
+Last modified: Fri Jun 6 2008
 <!-- hhmts end -->
 </p>
 
diff --git a/README b/README
new file mode 100644
index 000000000..e69de29bb
--- /dev/null
+++ b/README