Import the recently announced changes:

Import the recently announced changes:
1. Carpet has now an infrastructure for multiple maps (aka "grid
patches").  Instead of a single grid hierarchy there can now be
several.  This is largely untested, because the remainder of Cactus
cannot handle multiple coordinate systems.

2. The order in which the schedule bins are called has changed.  As Ian
Hawke pointed out, the previous order during time evolution was
inconsistent.  The initial data ordering did not allow for recovering
and was not usable for progressively solving elliptic equations for
initial data.

3. Carpet now supports convergence levels.  The convergence level
specifies by how many factors of two the resolution in the parameter
file should be coarsened (or refined, if negative).  This should make
convergence tests and test runs much easier.  It is, in principle, also
possible to run several convergence levels at once.  This has not been
tested because the remainder of Cactus cannot handle multiple
resolutions.  This will be necessary for a multigrid solver, and also
for having a shadow hierarchy to determine where to refine adaptively.

4. Carpet works together with the new CoordBase domain specification
parameters.  Without these, using convergence levels will lead to very
strange results.

5. The "modes" have changed.  There are now:
        meta mode:      the whole simulation
        global mode:    one convergence level
        level mode:     one refinement level
        singlemap mode: one map on one refinement level
        local mode:     as previously
The whole mode handling has been cleaned up.

6. The regridding thorn has been cleaned up.

7. The kind of prolongation stencil is now determined in Carpet, i.e. at
a fairly hight level, instead of in CarpetLib.

8. The low-order prolongation operators have been made much more
efficient (as have previously the higher-order ones).

9. Assorted smaller changes.

For Carpet users, there should be no major incompatibilities.  The major
improvements are 3 and 4 combined.  Here is an example:

CoordBase::domainsize    = extent
CoordBase::spacing       = gridspacing
CoordBase::zero_origin_x = yes
CoordBase::zero_origin_y = yes
CoordBase::zero_origin_z = yes
CoordBase::xextent       = 20.0
CoordBase::yextent       = 20.0
CoordBase::zextent       = 20.0
CoordBase::dx            = 1.0
CoordBase::dy            = 1.0
CoordBase::dz            = 1.0

CoordBase::boundary_shiftout_x_lower = 1
CoordBase::boundary_shiftout_y_lower = 1
CoordBase::boundary_shiftout_z_lower = 1

Carpet::domain_from_coordbase   = yes
Carpet::convergence_level       = 0

grid::type         = coordbase
grid::domain       = octant
grid::avoid_origin = no

This gives you a grid that extends from the origin ("zero_origin") up to
20.0 with a grid spacing of 1.0.  Symmetry zones and boundary zones are
added automatically.  The "shiftout" says that there is no boundary
point on the origin.  The staggering parameters (not shown) default to
"no".  In order to change the resolution, only the convergence level
has to be adjusted.  Note that the old way of specifying the domain
extent still works.

For Carpet developers, one major change is the new mode handling.  As
described in 5, the looping macros (that loop over all refinement
levels, or all components) have changed.

darcs-hash:20040125135727-07bb3-51c9647c1b5080e7e180b52a1b81fa155cfd19e9.gz

12 files changed, 512 insertions, 1075 deletions

diff --git a/Carpet/CarpetRegrid/interface.ccl b/Carpet/CarpetRegrid/interface.ccl
index 3c6352fda..efc701a73 100644
--- a/Carpet/CarpetRegrid/interface.ccl
+++ b/Carpet/CarpetRegrid/interface.ccl
@@ -1,5 +1,5 @@
 # Interface definition for thorn CarpetRegrid
-# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/interface.ccl,v 1.5 2004/01/15 09:45:58 cott Exp $
+# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/interface.ccl,v 1.6 2004/01/25 14:57:30 schnetter Exp $
 
 implements: CarpetRegrid
 
@@ -14,8 +14,31 @@ uses include header: vect.hh
 uses include header: gf.hh
 uses include header: gh.hh
 
-CCTK_INT FUNCTION RegridLevel (CCTK_POINTER_TO_CONST IN cctkGH, \
-                               CCTK_INT IN current_reflevel, \
+
+
+# The true prototype of the routine below:
+#   int Carpet_Regrid (const cGH * cctkGH,
+#                      gh<dim>::rexts  * bbsss,
+#                      gh<dim>::rbnds  * obss,
+#                      gh<dim>::rprocs * pss);
+CCTK_INT FUNCTION Carpet_Regrid (CCTK_POINTER_TO_CONST IN cctkGH,  \
+                                 CCTK_INT IN reflevel,             \
+                                 CCTK_INT IN map,                  \
+                                 CCTK_INT IN size,                 \
+                                 CCTK_INT IN ARRAY nboundaryzones, \
+                                 CCTK_INT IN ARRAY is_internal,    \
+                                 CCTK_INT IN ARRAY is_staggered,   \
+                                 CCTK_INT IN ARRAY shiftout,       \
+                                 CCTK_POINTER IN bsss,             \
+                                 CCTK_POINTER IN obss,             \
+                                 CCTK_POINTER IN pss)
+PROVIDES FUNCTION Carpet_Regrid WITH CarpetRegrid_Regrid LANGUAGE C
+
+
+
+
+CCTK_INT FUNCTION RegridLevel (CCTK_POINTER_TO_CONST IN cctkGH,  \
+                               CCTK_INT IN current_reflevel,     \
 			       CCTK_INT IN current_max_reflevel, \
 			       CCTK_INT IN max_reflevels)
 USES FUNCTION RegridLevel
diff --git a/Carpet/CarpetRegrid/param.ccl b/Carpet/CarpetRegrid/param.ccl
index 9f2dc9bf9..1de876c73 100644
--- a/Carpet/CarpetRegrid/param.ccl
+++ b/Carpet/CarpetRegrid/param.ccl
@@ -1,5 +1,5 @@
 # Parameter definitions for thorn CarpetRegrid
-# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/param.ccl,v 1.13 2004/01/13 13:50:05 hawke Exp $
+# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/param.ccl,v 1.14 2004/01/25 14:57:30 schnetter Exp $
 
 
 
@@ -17,16 +17,19 @@ CCTK_INT regrid_every "Regrid every n time steps" STEERABLE=always
    1:* :: "regrid every n time steps"
 } 0
 
-CCTK_INT activate_newlevels_on_regrid "When regridding, activate this many new levels (if possible)"
-# Note that this feature will STEER the parameter refinement_levels!!!
-#   It does this by adding activate_newlevels_on_regrid to refinement levels
+
+
+CCTK_INT activate_newlevels_on_regrid "When regridding, activate this many new levels (if possible).  Note that this will steer the parameter refinement_levels." STEERABLE=always
 {
-  : :: "Number of new levels to activate.  Negative numbers mean to de-activate"
+  : :: "Number of new levels to activate (negative numbers deactivate)"
 } 0
 
+CCTK_INT activate_next "The next iteration at which new levels should be activated" STEERABLE=always
+{
+  0: :: "Note that this parameter is steered when new levels are activated"
+} 1
 
-
-BOOLEAN outside_boundary_points[3] "On finer grids, where the upper grid boundary is adjacent to the outer boundary, put points outside the outer boundary (needed e.g. for periodicity)"
+BOOLEAN regrid_from_function "Regridding criteria from the aliased function?"
 {
 } "no"
 
@@ -45,6 +48,20 @@ KEYWORD refined_regions "Regions where the grid is refined" STEERABLE=always
 
 
 
+# Region specifications for centre refinement
+
+BOOLEAN symmetry_x "Refine the lower half in x-direction"
+{
+} "no"
+BOOLEAN symmetry_y "Refine the lower half in y-direction"
+{
+} "no"
+BOOLEAN symmetry_z "Refine the lower half in z-direction"
+{
+} "no"
+
+
+
 # Region specifications for manual gridpoint refinement
 
 CCTK_INT l1ixmin "Lower boundary of level 1 box in x-direction" STEERABLE=always
@@ -268,10 +285,3 @@ CCTK_STRING errorvar "Name of grid function that contains the error" STEERABLE=a
 {
   ".*" :: "must be the name of a grid function"
 } ""
-
-
-# Should we set the regridding criteria from an aliased function?
-
-BOOLEAN regrid_from_function "Regridding criteria from the aliased function?"
-{
-} "no"
diff --git a/Carpet/CarpetRegrid/src/automatic.cc b/Carpet/CarpetRegrid/src/automatic.cc
index a3b4da055..6add8c770 100644
--- a/Carpet/CarpetRegrid/src/automatic.cc
+++ b/Carpet/CarpetRegrid/src/automatic.cc
@@ -16,7 +16,7 @@
 #include "regrid.hh"
 
 extern "C" {
-  static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/src/automatic.cc,v 1.5 2004/08/04 16:25:58 schnetter Exp $";
+  static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/src/automatic.cc,v 1.1 2004/01/25 14:57:30 schnetter Exp $";
   CCTK_FILEVERSION(Carpet_CarpetRegrid_automatic_cc);
 }
 
@@ -31,12 +31,22 @@ namespace CarpetRegrid {
   
   int Automatic (cGH const * const cctkGH,
                  gh<dim> const & hh,
+                 int const reflevel,
+                 int const map,
+                 int const size,
+                 jjvect const & nboundaryzones,
+                 jjvect const & is_internal,
+                 jjvect const & is_staggered,
+                 jjvect const & shiftout,
                  gh<dim>::rexts  & bbsss,
                  gh<dim>::rbnds  & obss,
                  gh<dim>::rprocs & pss)
   {
     DECLARE_CCTK_PARAMETERS;
     
+    assert (reflevel>=0 && reflevel<maxreflevels);
+    assert (map>=0 && map<maps);
+    
     assert (refinement_levels >= 1);
     
     assert (bbsss.size() >= 1);
@@ -54,19 +64,15 @@ namespace CarpetRegrid {
     assert (CCTK_VarTypeI(vi) == CCTK_VARIABLE_REAL);
     assert (CCTK_GroupDimI(gi) == dim);
     
-    assert (arrdata.at(gi).at(Carpet::map).data.at(vi-v1));
-    const gf<CCTK_REAL,dim>& errorgf
+    assert (arrdata.at(gi).at(map).data.at(vi-v1));
+    const gf<CCTK_REAL,dim>& errorvar
       = (*dynamic_cast<const gf<CCTK_REAL,dim>*>
-         (arrdata.at(gi).at(Carpet::map).data.at(vi-v1)));
-    
-    assert (! smart_outer_boundaries);
+         (arrdata.at(gi).at(map).data.at(vi-v1)));
     
     vector<ibbox> bbs;
     gh<dim>::cbnds obs;
     Automatic_OneLevel
-      (cctkGH, hh,
-       reflevel, min(reflevels+1, maxreflevels),
-       minwidth, minfraction, maxerror, errorgf,
+      (cctkGH, hh, reflevel, minwidth, minfraction, maxerror, errorvar,
        bbs, obs);
     
     // make multiprocessor aware
@@ -75,7 +81,10 @@ namespace CarpetRegrid {
     
     // make multigrid aware
     vector<vector<ibbox> > bbss;
-    MakeMultigridBoxes (cctkGH, bbs, obs, bbss);
+    MakeMultigridBoxes
+      (cctkGH,
+       size, nboundaryzones, is_internal, is_staggered, shiftout,
+       bbs, obs, bbss);
     
     
     
@@ -106,32 +115,31 @@ namespace CarpetRegrid {
   
   void Automatic_OneLevel (const cGH * const cctkGH,
                            const gh<dim> & hh,
-                           const int rl,
-                           const int numrl,
+                           const int reflevel,
                            const int minwidth,
                            const CCTK_REAL minfraction,
                            const CCTK_REAL maxerror,
-                           const gf<CCTK_REAL,dim> & errorgf,
+                           const gf<CCTK_REAL,dim> & errorvar,
                            vector<ibbox> & bbs,
                            vector<bbvect> & obs)
   {
-    if (rl+1 >= numrl) return;
+    if (reflevel+1 >= maxreflevels) return;
     
     // Arbitrary
     const int tl = 0;
     const int ml = 0;
     
-//     cout << endl << "MRA: Choosing regions to refine in " << hh.components(rl) << " components" << endl;
+//     cout << endl << "MRA: Choosing regions to refine in " << hh.components(reflevel) << " components" << endl;
     
     list<ibbox> bbl;
-    for (int c=0; c<hh.components(rl); ++c) {
-      const ibbox region = hh.extents.at(rl).at(c).at(ml);
+    for (int c=0; c<hh.components(reflevel); ++c) {
+      const ibbox region = hh.extents.at(reflevel).at(c).at(ml);
       assert (! region.empty());
       
-      const data<CCTK_REAL,dim>& errordata = *errorgf(tl,rl,c,ml);
+      const data<CCTK_REAL,dim>& errdata = *errorvar(tl,reflevel,c,ml);
       
       Automatic_Recursive (cctkGH, hh, minwidth, minfraction, maxerror,
-                           errordata, bbl, region);
+                           errdata, bbl, region);
     }
     
 //     int numpoints = 0;
@@ -177,7 +185,7 @@ namespace CarpetRegrid {
                             const int minwidth,
                             const CCTK_REAL minfraction,
                             const CCTK_REAL maxerror,
-                            const data<CCTK_REAL,dim> & errordata,
+                            const data<CCTK_REAL,dim> & errorvar,
                             list<ibbox> & bbl,
                             const ibbox & region)
   {
@@ -188,12 +196,8 @@ namespace CarpetRegrid {
     // (this doesn't work yet on multiple processors)
     assert (CCTK_nProcs(cctkGH)==1);
     int cnt = 0;
-    {
-      ibbox::iterator it=region.begin();
-      do {
-        if (errordata[*it] > maxerror) ++cnt;
-        ++it;
-      } while (it!=region.end());
+    for (ibbox::iterator it=region.begin(); it!=region.end(); ++it) {
+      if (errorvar[*it] > maxerror) ++cnt;
     }
     const CCTK_REAL fraction = (CCTK_REAL)cnt / region.size();
     const int width = maxval(region.shape() / region.stride());
@@ -227,9 +231,9 @@ namespace CarpetRegrid {
       assert (region1 + region2 == region);
       list<ibbox> bbl1, bbl2;
       Automatic_Recursive (cctkGH, hh, minwidth, minfraction, maxerror,
-                           errordata, bbl1, region1);
+                           errorvar, bbl1, region1);
       Automatic_Recursive (cctkGH, hh, minwidth, minfraction, maxerror,
-                           errordata, bbl2, region2);
+                           errorvar, bbl2, region2);
       // Combine regions if possible
       up2 += str-str/reffact;
       up2[dir] = lo2[dir];
diff --git a/Carpet/CarpetRegrid/src/baselevel.cc b/Carpet/CarpetRegrid/src/baselevel.cc
index c380c6844..18cc8a95c 100644
--- a/Carpet/CarpetRegrid/src/baselevel.cc
+++ b/Carpet/CarpetRegrid/src/baselevel.cc
@@ -9,7 +9,7 @@
 #include "regrid.hh"
 
 extern "C" {
-  static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/src/baselevel.cc,v 1.2 2004/04/18 13:29:43 schnetter Exp $";
+  static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/src/baselevel.cc,v 1.1 2004/01/25 14:57:30 schnetter Exp $";
   CCTK_FILEVERSION(Carpet_CarpetRegrid_baselevel_cc);
 }
 
@@ -24,12 +24,22 @@ namespace CarpetRegrid {
   
   int BaseLevel (cGH const * const cctkGH,
                  gh<dim> const & hh,
+                 int const reflevel,
+                 int const map,
+                 int const size,
+                 jjvect const & nboundaryzones,
+                 jjvect const & is_internal,
+                 jjvect const & is_staggered,
+                 jjvect const & shiftout,
                  gh<dim>::rexts  & bbsss,
                  gh<dim>::rbnds  & obss,
                  gh<dim>::rprocs & pss)
   {
     DECLARE_CCTK_PARAMETERS;
     
+    assert (reflevel>=0 && reflevel<maxreflevels);
+    assert (map>=0 && map<maps);
+    
     assert (refinement_levels == 1);
     
     assert (bbsss.size() == 1);
diff --git a/Carpet/CarpetRegrid/src/centre.cc b/Carpet/CarpetRegrid/src/centre.cc
index 7599324ea..46a23f794 100644
--- a/Carpet/CarpetRegrid/src/centre.cc
+++ b/Carpet/CarpetRegrid/src/centre.cc
@@ -9,7 +9,7 @@
 #include "regrid.hh"
 
 extern "C" {
-  static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/src/centre.cc,v 1.4 2004/04/28 15:45:25 schnetter Exp $";
+  static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/src/centre.cc,v 1.1 2004/01/25 14:57:30 schnetter Exp $";
   CCTK_FILEVERSION(Carpet_CarpetRegrid_centre_cc);
 }
 
@@ -24,16 +24,26 @@ namespace CarpetRegrid {
   
   int Centre (cGH const * const cctkGH,
               gh<dim> const & hh,
+              int const reflevel,
+              int const map,
+              int const size,
+              jjvect const & nboundaryzones,
+              jjvect const & is_internal,
+              jjvect const & is_staggered,
+              jjvect const & shiftout,
               gh<dim>::rexts  & bbsss,
               gh<dim>::rbnds  & obss,
               gh<dim>::rprocs & pss)
   {
     DECLARE_CCTK_PARAMETERS;
     
+    assert (reflevel>=0 && reflevel<maxreflevels);
+    assert (map>=0 && map<maps);
+    
     assert (refinement_levels >= 1);
     
     // do nothing if the levels already exist
-    if (reflevel == refinement_levels) return 0;
+    if (bbsss.size() == refinement_levels) return 0;
     
     assert (bbsss.size() >= 1);
     
@@ -48,8 +58,6 @@ namespace CarpetRegrid {
     ivect rlb  = hh.baseextent.lower();
     ivect rub  = hh.baseextent.upper();
     
-    assert (! smart_outer_boundaries);
-    
     for (size_t rl=1; rl<bbsss.size(); ++rl) {
       
       // save old values
@@ -81,7 +89,10 @@ namespace CarpetRegrid {
       
       // make multigrid aware
       vector<vector<ibbox> > bbss;
-      MakeMultigridBoxes (cctkGH, bbs, obs, bbss);
+      MakeMultigridBoxes
+        (cctkGH,
+         size, nboundaryzones, is_internal, is_staggered, shiftout,
+         bbs, obs, bbss);
       
       bbsss.at(rl) = bbss;
       obss.at(rl) = obs;
diff --git a/Carpet/CarpetRegrid/src/make.code.defn b/Carpet/CarpetRegrid/src/make.code.defn
index ec45ed348..0b692a316 100644
--- a/Carpet/CarpetRegrid/src/make.code.defn
+++ b/Carpet/CarpetRegrid/src/make.code.defn
@@ -1,8 +1,16 @@
 # Main make.code.defn file for thorn CarpetRegrid
-# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/src/make.code.defn,v 1.2 2002/05/16 23:25:54 schnetter Exp $
+# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/src/make.code.defn,v 1.3 2004/01/25 14:57:30 schnetter Exp $
 
 # Source files in this directory
-SRCS = regrid.cc paramcheck.cc
+SRCS =  automatic.cc            \
+        baselevel.cc            \
+        centre.cc               \
+        manualcoordinatelist.cc \
+        manualcoordinates.cc    \
+        manualgridpointlist.cc  \
+        manualgridpoints.cc     \
+        regrid.cc               \
+        paramcheck.cc
 
 # Subdirectories containing source files
 SUBDIRS = 
diff --git a/Carpet/CarpetRegrid/src/manualcoordinatelist.cc b/Carpet/CarpetRegrid/src/manualcoordinatelist.cc
index 10e40ecbd..92ab6900a 100644
--- a/Carpet/CarpetRegrid/src/manualcoordinatelist.cc
+++ b/Carpet/CarpetRegrid/src/manualcoordinatelist.cc
@@ -1,7 +1,6 @@
-#include <cassert>
-#include <cmath>
-#include <cstring>
-#include <sstream>
+#include <assert.h>
+#include <string.h>
+
 #include <vector>
 
 #include "cctk.h"
@@ -13,7 +12,7 @@
 #include "regrid.hh"
 
 extern "C" {
-  static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/src/manualcoordinatelist.cc,v 1.12 2004/08/14 07:42:00 schnetter Exp $";
+  static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/src/manualcoordinatelist.cc,v 1.1 2004/01/25 14:57:30 schnetter Exp $";
   CCTK_FILEVERSION(Carpet_CarpetRegrid_manualcoordinatelist_cc);
 }
 
@@ -28,35 +27,29 @@ namespace CarpetRegrid {
   
   int ManualCoordinateList (cGH const * const cctkGH,
                             gh<dim> const & hh,
+                            int const reflevel,
+                            int const map,
+                            int const size,
+                            jjvect const & nboundaryzones,
+                            jjvect const & is_internal,
+                            jjvect const & is_staggered,
+                            jjvect const & shiftout,
                             gh<dim>::rexts  & bbsss,
                             gh<dim>::rbnds  & obss,
                             gh<dim>::rprocs & pss)
   {
     DECLARE_CCTK_PARAMETERS;
-    int ierr;
+    
+    assert (reflevel>=0 && reflevel<maxreflevels);
+    assert (map>=0 && map<maps);
     
     assert (refinement_levels >= 1);
     
     // do nothing if the levels already exist
-    if (reflevel == refinement_levels) return 0;
+    if (bbsss.size() == refinement_levels) return 0;
     
     assert (bbsss.size() >= 1);
     
-    jjvect nboundaryzones, is_internal, is_staggered, shiftout;
-    ierr = GetBoundarySpecification
-      (2*dim, &nboundaryzones[0][0], &is_internal[0][0],
-       &is_staggered[0][0], &shiftout[0][0]);
-    assert (!ierr);
-    rvect physical_min, physical_max;
-    rvect interior_min, interior_max;
-    rvect exterior_min, exterior_max;
-    rvect base_spacing;
-    ierr = GetDomainSpecification
-      (dim, &physical_min[0], &physical_max[0],
-       &interior_min[0], &interior_max[0],
-       &exterior_min[0], &exterior_max[0], &base_spacing[0]);
-    assert (!ierr);
-    
     bbsss.resize (refinement_levels);
     obss.resize (refinement_levels);
     pss.resize (refinement_levels);
@@ -70,79 +63,39 @@ namespace CarpetRegrid {
         CCTK_WARN (0, "Could not parse parameter \"coordinates\"");
       }
     }
-
+    
     vector<vector<bbvect> > newobss;
-    if (smart_outer_boundaries) {
-      // TODO:
-      // assert (domain_from_coordbase);
-      
+    if (strcmp(outerbounds, "") !=0 ) {
+      istringstream ob_str (outerbounds);
+      try {
+        ob_str >> newobss;
+      } catch (input_error) {
+        CCTK_WARN (0, "Could not parse parameter \"outerbounds\"");
+      }
+      bool good = newobss.size() == newbbss.size();
+      if (good) {
+        for (size_t rl=0; rl<newobss.size(); ++rl) {
+          good = good && newobss.at(rl).size() == newbbss.at(rl).size();
+        }
+      }
+      if (! good) {
+        cout << "coordinates: " << newbbss << endl;
+        cout << "outerbounds: " << newobss << endl;
+        CCTK_WARN (0, "The parameters \"outerbounds\" and \"coordinates\" must have the same structure");
+      }
+    } else {
       newobss.resize(newbbss.size());
       for (size_t rl=0; rl<newobss.size(); ++rl) {
         newobss.at(rl).resize(newbbss.at(rl).size());
         for (size_t c=0; c<newobss.at(rl).size(); ++c) {
-          for (int d=0; d<dim; ++d) {
-            assert (mglevel==0);
-            rvect const spacing = base_spacing * pow(CCTK_REAL(mgfact), basemglevel) / ipow(reffact, rl+1);
-            ierr = ConvertFromPhysicalBoundary
-              (dim, &physical_min[0], &physical_max[0],
-               &interior_min[0], &interior_max[0],
-               &exterior_min[0], &exterior_max[0], &spacing[0]);
-            assert (!ierr);
-            newobss.at(rl).at(c)[d][0] = abs(newbbss.at(rl).at(c).lower()[d] - physical_min[d]) < 1.0e-6 * spacing[d];
-            if (newobss.at(rl).at(c)[d][0]) {
-              rvect lo = newbbss.at(rl).at(c).lower();
-              rvect up = newbbss.at(rl).at(c).upper();
-              rvect str = newbbss.at(rl).at(c).stride();
-              lo[d] = exterior_min[d];
-              newbbss.at(rl).at(c) = rbbox(lo, up, str);
-            }
-            newobss.at(rl).at(c)[d][1] = abs(newbbss.at(rl).at(c).upper()[d] - physical_max[d]) < 1.0e-6 * base_spacing[d] / ipow(reffact, rl);
-            if (newobss.at(rl).at(c)[d][1]) {
-              rvect lo = newbbss.at(rl).at(c).lower();
-              rvect up = newbbss.at(rl).at(c).upper();
-              rvect str = newbbss.at(rl).at(c).stride();
-              up[d] = exterior_max[d];
-              newbbss.at(rl).at(c) = rbbox(lo, up, str);
-            }
-          }
+          newobss.at(rl).at(c) = bbvect(false);
         }
       }
-      
-    } else {                    // if ! smart_outer_boundaries
-      
-      if (strcmp(outerbounds, "") !=0 ) {
-        istringstream ob_str (outerbounds);
-        try {
-          ob_str >> newobss;
-        } catch (input_error) {
-          CCTK_WARN (0, "Could not parse parameter \"outerbounds\"");
-        }
-        bool good = newobss.size() == newbbss.size();
-        if (good) {
-          for (size_t rl=0; rl<newobss.size(); ++rl) {
-            good = good && newobss.at(rl).size() == newbbss.at(rl).size();
-          }
-        }
-        if (! good) {
-          cout << "coordinates: " << newbbss << endl;
-          cout << "outerbounds: " << newobss << endl;
-          CCTK_WARN (0, "The parameters \"outerbounds\" and \"coordinates\" must have the same structure");
-        }
-      } else {
-        newobss.resize(newbbss.size());
-        for (size_t rl=0; rl<newobss.size(); ++rl) {
-          newobss.at(rl).resize(newbbss.at(rl).size());
-          for (size_t c=0; c<newobss.at(rl).size(); ++c) {
-            newobss.at(rl).at(c) = bbvect(false);
-          }
-        }
-      }
-
-    } // if ! smart_outer_boundaries
+    }
     
     if (newbbss.size() < refinement_levels-1) {
       CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING,
-                  "The parameter \"coordinates\" must contain at least \"refinement_levels-1\" (here: %d) levels", int(refinement_levels-1));
+                  "The parameter \"coordinates\" must contain at least \"refinement_levels-1\" (here: %d) levels", (int)refinement_levels-1);
     }
     
     for (size_t rl=1; rl<refinement_levels; ++rl) {
@@ -156,10 +109,6 @@ namespace CarpetRegrid {
       for (size_t c=0; c<newbbss.at(rl-1).size(); ++c) {
         rbbox const & ext = newbbss.at(rl-1).at(c);
         bbvect const & ob = newobss.at(rl-1).at(c);
-        // TODO:
-        // assert (domain_from_coordbase);
-        rvect const spacing = base_spacing * pow(CCTK_REAL(mgfact), basemglevel) / ipow(reffact, rl);
-        assert (all(abs(ext.stride() - spacing) < spacing * 1.0e-10));
         ManualCoordinates_OneLevel
           (cctkGH, hh, rl, refinement_levels,
            ext.lower(), ext.upper(), ob, bbs, obs);
@@ -171,7 +120,10 @@ namespace CarpetRegrid {
       
       // make multigrid aware
       vector<vector<ibbox> > bbss;
-      MakeMultigridBoxes (cctkGH, bbs, obs, bbss);
+      MakeMultigridBoxes
+        (cctkGH,
+         size, nboundaryzones, is_internal, is_staggered, shiftout,
+         bbs, obs, bbss);
       
       bbsss.at(rl) = bbss;
       obss.at(rl) = obs;
diff --git a/Carpet/CarpetRegrid/src/manualcoordinates.cc b/Carpet/CarpetRegrid/src/manualcoordinates.cc
index 91018cf95..1e3719f18 100644
--- a/Carpet/CarpetRegrid/src/manualcoordinates.cc
+++ b/Carpet/CarpetRegrid/src/manualcoordinates.cc
@@ -11,7 +11,7 @@
 #include "regrid.hh"
 
 extern "C" {
-  static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/src/manualcoordinates.cc,v 1.5 2004/04/28 15:45:25 schnetter Exp $";
+  static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/src/manualcoordinates.cc,v 1.1 2004/01/25 14:57:30 schnetter Exp $";
   CCTK_FILEVERSION(Carpet_CarpetRegrid_manualcoordinates_cc);
 }
 
@@ -26,19 +26,29 @@ namespace CarpetRegrid {
   
   int ManualCoordinates (cGH const * const cctkGH,
                          gh<dim> const & hh,
+                         int const reflevel,
+                         int const map,
+                         int const size,
+                         jjvect const & nboundaryzones,
+                         jjvect const & is_internal,
+                         jjvect const & is_staggered,
+                         jjvect const & shiftout,
                          gh<dim>::rexts  & bbsss,
                          gh<dim>::rbnds  & obss,
                          gh<dim>::rprocs & pss)
   {
     DECLARE_CCTK_PARAMETERS;
     
+    assert (reflevel>=0 && reflevel<maxreflevels);
+    assert (map>=0 && map<maps);
+    
     if (refinement_levels > 4) {
       CCTK_WARN (0, "Cannot currently specify manual refinement regions for more than 4 refinement levels");
     }
     assert (refinement_levels >= 1 && refinement_levels <= 4);
     
     // do nothing if the levels already exist
-    if (reflevel == refinement_levels) return 0;
+    if (bbsss.size() == refinement_levels) return 0;
     
     assert (bbsss.size() >= 1);
     
@@ -54,8 +64,6 @@ namespace CarpetRegrid {
     lower.at(2) = rvect (l3xmin, l3ymin, l3zmin);
     upper.at(2) = rvect (l3xmax, l3ymax, l3zmax);
     
-    assert (! smart_outer_boundaries);
-    
     for (size_t rl=1; rl<bbsss.size(); ++rl) {
       
       bbvect const ob (false);
@@ -73,7 +81,10 @@ namespace CarpetRegrid {
       
       // make multigrid aware
       vector<vector<ibbox> > bbss;
-      MakeMultigridBoxes (cctkGH, bbs, obs, bbss);
+      MakeMultigridBoxes
+        (cctkGH,
+         size, nboundaryzones, is_internal, is_staggered, shiftout,
+         bbs, obs, bbss);
       
       bbsss.at(rl) = bbss;
       obss.at(rl) = obs;
@@ -88,27 +99,27 @@ namespace CarpetRegrid {
   
   void ManualCoordinates_OneLevel (const cGH * const cctkGH,
                                    const gh<dim> & hh,
-                                   const int rl,
-                                   const int numrl,
+                                   const int reflevel,
+                                   const int reflevels,
                                    const rvect lower,
                                    const rvect upper,
                                    const bbvect obound,
                                    vector<ibbox> & bbs,
                                    vector<bbvect> & obs)
   {
-    if (rl >= numrl) return;
+    if (reflevel >= reflevels) return;
     
-    jvect const ilower = pos2int (cctkGH, hh, lower, rl);
-    jvect const iupper = pos2int (cctkGH, hh, upper, rl);
+    jvect const ilower = pos2int (cctkGH, hh, lower, reflevel);
+    jvect const iupper = pos2int (cctkGH, hh, upper, reflevel);
     
     ManualGridpoints_OneLevel
-      (cctkGH, hh, rl, numrl, ilower, iupper, obound, bbs, obs);
+      (cctkGH, hh, reflevel, reflevels, ilower, iupper, obound, bbs, obs);
   }
   
   
   
   ivect delta2int (const cGH * const cctkGH, const gh<dim>& hh,
-                   const rvect & rpos, const int rl)
+                   const rvect & rpos, const int reflevel)
   {
     rvect global_lower, global_upper;
     for (int d=0; d<dim; ++d) {
@@ -121,13 +132,15 @@ namespace CarpetRegrid {
     }
     const ivect global_extent (hh.baseextent.upper() - hh.baseextent.lower());
     
+    CCTK_REAL (* const rfloor) (CCTK_REAL const) = floor;
+    
     const rvect scale  = rvect(global_extent) / (global_upper - global_lower);
-    const int levfac = ipow(hh.reffact, rl);
+    const int levfac = ipow(hh.reffact, reflevel);
     assert (all (hh.baseextent.stride() % levfac == 0));
     const ivect istride = hh.baseextent.stride() / levfac;
     
     const ivect ipos
-      = ivect(floor(rpos * scale / rvect(istride) + 0.5)) * istride;
+      = ivect(::map(rfloor, rpos * scale / rvect(istride) + 0.5)) * istride;
     
     const rvect apos = rpos * scale;
     assert (all(abs(apos - rvect(ipos)) < rvect(istride)*0.01));
@@ -138,7 +151,7 @@ namespace CarpetRegrid {
   
   
   ivect pos2int (const cGH* const cctkGH, const gh<dim>& hh,
-                 const rvect & rpos, const int rl)
+                 const rvect & rpos, const int reflevel)
   {
     rvect global_lower, global_upper;
     for (int d=0; d<dim; ++d) {
@@ -151,14 +164,16 @@ namespace CarpetRegrid {
     }
     const ivect global_extent (hh.baseextent.upper() - hh.baseextent.lower());
     
+    CCTK_REAL (* const rfloor) (CCTK_REAL const) = floor;
+    
     const rvect scale  = rvect(global_extent) / (global_upper - global_lower);
-    const int levfac = ipow(hh.reffact, rl);
+    const int levfac = ipow(hh.reffact, reflevel);
     assert (all (hh.baseextent.stride() % levfac == 0));
     const ivect istride = hh.baseextent.stride() / levfac;
     
     const ivect ipos
-      = (ivect(floor((rpos - global_lower) * scale / rvect(istride) + 0.5))
-         * istride);
+      = ivect(::map(rfloor, (rpos - global_lower) * scale / rvect(istride)
+                    + 0.5)) * istride;
     
     return ipos;
   }
diff --git a/Carpet/CarpetRegrid/src/manualgridpointlist.cc b/Carpet/CarpetRegrid/src/manualgridpointlist.cc
index e53866e6c..6976b982b 100644
--- a/Carpet/CarpetRegrid/src/manualgridpointlist.cc
+++ b/Carpet/CarpetRegrid/src/manualgridpointlist.cc
@@ -1,7 +1,6 @@
 #include <assert.h>
 #include <string.h>
 
-#include <sstream>
 #include <vector>
 
 #include "cctk.h"
@@ -13,7 +12,7 @@
 #include "regrid.hh"
 
 extern "C" {
-  static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/src/manualgridpointlist.cc,v 1.4 2004/07/02 10:14:51 tradke Exp $";
+  static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/src/manualgridpointlist.cc,v 1.1 2004/01/25 14:57:30 schnetter Exp $";
   CCTK_FILEVERSION(Carpet_CarpetRegrid_manualgridpointlist_cc);
 }
 
@@ -28,16 +27,26 @@ namespace CarpetRegrid {
   
   int ManualGridpointList (cGH const * const cctkGH,
                            gh<dim> const & hh,
+                           int const reflevel,
+                           int const map,
+                           int const size,
+                           jjvect const & nboundaryzones,
+                           jjvect const & is_internal,
+                           jjvect const & is_staggered,
+                           jjvect const & shiftout,
                            gh<dim>::rexts  & bbsss,
                            gh<dim>::rbnds  & obss,
                            gh<dim>::rprocs & pss)
   {
     DECLARE_CCTK_PARAMETERS;
     
+    assert (reflevel>=0 && reflevel<maxreflevels);
+    assert (map>=0 && map<maps);
+    
     assert (refinement_levels >= 1);
     
     // do nothing if the levels already exist
-    if (reflevel == refinement_levels) return 0;
+    if (bbsss.size() == refinement_levels) return 0;
     
     assert (bbsss.size() >= 1);
     
@@ -111,7 +120,10 @@ namespace CarpetRegrid {
       
       // make multigrid aware
       vector<vector<ibbox> > bbss;
-      MakeMultigridBoxes (cctkGH, bbs, obs, bbss);
+      MakeMultigridBoxes
+        (cctkGH,
+         size, nboundaryzones, is_internal, is_staggered, shiftout,
+         bbs, obs, bbss);
       
       bbsss.at(rl) = bbss;
       obss.at(rl) = obs;
diff --git a/Carpet/CarpetRegrid/src/manualgridpoints.cc b/Carpet/CarpetRegrid/src/manualgridpoints.cc
index 2a8bedb5e..e8cd6c17a 100644
--- a/Carpet/CarpetRegrid/src/manualgridpoints.cc
+++ b/Carpet/CarpetRegrid/src/manualgridpoints.cc
@@ -1,6 +1,5 @@
 #include <assert.h>
 
-#include <sstream>
 #include <vector>
 
 #include "cctk.h"
@@ -12,7 +11,7 @@
 #include "regrid.hh"
 
 extern "C" {
-  static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/src/manualgridpoints.cc,v 1.5 2004/07/02 10:14:51 tradke Exp $";
+  static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/src/manualgridpoints.cc,v 1.1 2004/01/25 14:57:30 schnetter Exp $";
   CCTK_FILEVERSION(Carpet_CarpetRegrid_manualgridpoints_cc);
 }
 
@@ -27,19 +26,29 @@ namespace CarpetRegrid {
   
   int ManualGridpoints (cGH const * const cctkGH,
                         gh<dim> const & hh,
+                        int const reflevel,
+                        int const map,
+                        int const size,
+                        jjvect const & nboundaryzones,
+                        jjvect const & is_internal,
+                        jjvect const & is_staggered,
+                        jjvect const & shiftout,
                         gh<dim>::rexts  & bbsss,
                         gh<dim>::rbnds  & obss,
                         gh<dim>::rprocs & pss)
   {
     DECLARE_CCTK_PARAMETERS;
     
+    assert (reflevel>=0 && reflevel<maxreflevels);
+    assert (map>=0 && map<maps);
+    
     if (refinement_levels > 4) {
       CCTK_WARN (0, "Cannot currently specify manual refinement regions for more than 4 refinement levels");
     }
     assert (refinement_levels >= 1 && refinement_levels <= 4);
     
     // do nothing if the levels already exist
-    if (reflevel == refinement_levels) return 0;
+    if (bbsss.size() == refinement_levels) return 0;
     
     assert (bbsss.size() >= 1);
     
@@ -55,8 +64,6 @@ namespace CarpetRegrid {
     ilower.at(2) = ivect (l3ixmin, l3iymin, l3izmin);
     iupper.at(2) = ivect (l3ixmax, l3iymax, l3izmax);
     
-    assert (! smart_outer_boundaries);
-    
     for (size_t rl=1; rl<bbsss.size(); ++rl) {
       
       bbvect const ob (false);
@@ -65,7 +72,7 @@ namespace CarpetRegrid {
       gh<dim>::cbnds obs;
       
       ManualGridpoints_OneLevel
-        (cctkGH, hh, rl,refinement_levels,
+        (cctkGH, hh, rl, refinement_levels,
          ilower.at(rl-1), iupper.at(rl-1), ob, bbs, obs);
       
       // make multiprocessor aware
@@ -74,7 +81,10 @@ namespace CarpetRegrid {
       
       // make multigrid aware
       vector<vector<ibbox> > bbss;
-      MakeMultigridBoxes (cctkGH, bbs, obs, bbss);
+      MakeMultigridBoxes
+        (cctkGH,
+         size, nboundaryzones, is_internal, is_staggered, shiftout,
+         bbs, obs, bbss);
       
       bbsss.at(rl) = bbss;
       obss.at(rl) = obs;
@@ -89,36 +99,38 @@ namespace CarpetRegrid {
   
   void ManualGridpoints_OneLevel (const cGH * const cctkGH,
                                   const gh<dim> & hh,
-                                  const int rl,
-                                  const int numrl,
+                                  const int reflevel,
+                                  const int reflevels,
                                   const ivect ilower,
                                   const ivect iupper,
                                   const bbvect obound,
                                   vector<ibbox> & bbs,
                                   vector<bbvect> & obs)
   {
+    if (reflevel >= reflevels) return;
+    
     const ivect rstr = hh.baseextent.stride();
     const ivect rlb  = hh.baseextent.lower();
     const ivect rub  = hh.baseextent.upper();
     
-    const int levfac = ipow(hh.reffact, rl);
+    const int levfac = ipow(hh.reffact, reflevel);
     assert (all (rstr % levfac == 0));
     const ivect str (rstr / levfac);
     const ivect lb  (ilower);
     const ivect ub  (iupper);
     if (! all(lb>=rlb && ub<=rub)) {
       ostringstream buf;
-      buf << "The refinement region boundaries for refinement level #" << rl << " are not within the main grid.  Allowed are the grid point boundaries " << rlb << " - " << rub << "; specified were " << lb << " - " << ub << ends;
+      buf << "The refinement region boundaries for refinement level #" << reflevel << " are not within the main grid.  Allowed are the grid point boundaries " << rlb << " - " << rub << "; specified were " << lb << " - " << ub << ends;
       CCTK_WARN (0, buf.str().c_str());
     }
     if (! all(lb<=ub)) {
       ostringstream buf;
-      buf << "The refinement region boundaries for refinement level #" << rl << " have the upper boundary (" << ub << ") less than the lower boundary (" << lb << ")" << ends;
+      buf << "The refinement region boundaries for refinement level #" << reflevel << " have the upper boundary (" << ub << ") less than the lower boundary (" << lb << ")" << ends;
       CCTK_WARN (0, buf.str().c_str());
     }
     if (! all(lb%str==0 && ub%str==0)) {
       CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING,
-		  "The refinement region boundaries for refinement level #%d are not a multiple of the stride for that level", rl);
+		  "The refinement region boundaries for refinement level #%d are not a multiple of the stride for that level", reflevel);
     }
     assert (all(lb>=rlb && ub<=rub));
     assert (all(lb<=ub));
diff --git a/Carpet/CarpetRegrid/src/regrid.cc b/Carpet/CarpetRegrid/src/regrid.cc
index 82f6e0954..3c6b8b0a7 100644
--- a/Carpet/CarpetRegrid/src/regrid.cc
+++ b/Carpet/CarpetRegrid/src/regrid.cc
@@ -1,81 +1,70 @@
 #include <assert.h>
-#include <math.h>
-#include <stdarg.h>
-#include <stdlib.h>
-#include <string.h>
 
-#include <algorithm>
-#include <list>
 #include <sstream>
 #include <string>
-#include <vector>
 
 #include "cctk.h"
 #include "cctk_Parameters.h"
 
-#include "bbox.hh"
-#include "bboxset.hh"
-#include "defs.hh"
 #include "gh.hh"
-#include "gf.hh"
 #include "vect.hh"
 
 #include "carpet.hh"
 #include "regrid.hh"
 
 extern "C" {
-  static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/src/regrid.cc,v 1.33 2004/01/15 09:45:58 cott Exp $";
+  static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/src/regrid.cc,v 1.34 2004/01/25 14:57:30 schnetter Exp $";
   CCTK_FILEVERSION(Carpet_CarpetRegrid_regrid_cc);
 }
 
 
 
-#undef AUTOMATIC_BOUNDARIES
-
-
-
 namespace CarpetRegrid {
   
   using namespace std;
   using namespace Carpet;
   
-  int initial_refinement_levels = 0; // Store the initial number of 
-                                     // levels to refine for correct
-                                     // progressive MR.
-                                     // Note that this does not need
-                                     // to be checkpointed as it is
-                                     // reset at the beginning of every
-                                     // iteration.
-  int last_regridded_iteration = -1; // The last iteration where any
-                                     // regridding was done.
-  
-  int CarpetRegridStartup ()
-  {
-    RegisterRegridRoutine (CarpetRegridRegrid);
-    return 0;
-  }
-  
-  
-  
-  int CarpetRegridRegrid (const cGH * const cctkGH,
-			  gh<dim>::rexts& bbsss,
-			  gh<dim>::rbnds& obss,
-			  gh<dim>::rprocs& pss)
+  CCTK_INT CarpetRegrid_Regrid (CCTK_POINTER_TO_CONST const cctkGH_,
+                                CCTK_INT const reflevel,
+                                CCTK_INT const map,
+                                CCTK_INT const size,
+                                CCTK_INT const * const nboundaryzones_,
+                                CCTK_INT const * const is_internal_,
+                                CCTK_INT const * const is_staggered_,
+                                CCTK_INT const * const shiftout_,
+                                CCTK_POINTER const bbsss_,
+                                CCTK_POINTER const obss_,
+                                CCTK_POINTER const pss_)
   {
     DECLARE_CCTK_PARAMETERS;
     
+    const cGH * const cctkGH = (const cGH *) cctkGH_;
+    
+    assert (reflevel>=0 && reflevel<maxreflevels);
+    assert (map>=0 && map<maps);
+    
+    assert (size == 2*dim);
+    jjvect const nboundaryzones (* (jjvect const *) nboundaryzones_);
+    jjvect const is_internal    (* (jjvect const *) is_internal_);
+    jjvect const is_staggered   (* (jjvect const *) is_staggered_);
+    jjvect const shiftout       (* (jjvect const *) shiftout_);
+    
+    gh<dim>::rexts  & bbsss = * (gh<dim>::rexts  *) bbsss_;
+    gh<dim>::rbnds  & obss  = * (gh<dim>::rbnds  *) obss_;
+    gh<dim>::rprocs & pss   = * (gh<dim>::rprocs *) pss_;
+    
+    gh<dim> const & hh = *vhh.at(map);
+    
+    assert (is_meta_mode());
+    
+    
+    
     assert (regrid_every == -1 || regrid_every == 0
 	    || regrid_every % maxmglevelfact == 0);
     
     // Return if no regridding is desired
     if (regrid_every == -1) return 0;
     
-#if 0
-    // Return if this is not the main hierarchy
-    if (mglevel != 0) return 0;
-#endif
-    assert (mglevel == -1);
-    
     // Return if we want to regrid during initial data only, and this
     // is not the time for initial data
     if (regrid_every == 0 && cctkGH->cctk_iteration != 0) return 0;
@@ -85,833 +74,105 @@ namespace CarpetRegrid {
 	&& cctkGH->cctk_iteration % regrid_every != 0) {
       return 0;
     }
-
-    // Find the number of levels that should be active.
-    int newnumlevels;
-
-    // Is this the first time that we have regridded this iteration?
-    if (cctkGH->cctk_iteration > last_regridded_iteration) {
-      last_regridded_iteration = cctkGH->cctk_iteration;
-      initial_refinement_levels = refinement_levels;
+    
+    
+    
+    // Steer parameters
+    if (activate_newlevels_on_regrid != 0) {
+      if (cctkGH->cctk_iteration >= activate_next) {
+        const int newnumlevels
+          = min(refinement_levels + activate_newlevels_on_regrid,
+                maxreflevels);
+        assert (newnumlevels>0 && newnumlevels<=maxreflevels);
+        
+        *const_cast<CCTK_INT*>(&activate_next) = cctkGH->cctk_iteration + 1;
+        ostringstream next;
+        next << activate_next;
+        CCTK_ParameterSet
+          ("activate_next", "CarpetRegrid", next.str().c_str());
+        
+        *const_cast<CCTK_INT*>(&refinement_levels) = newnumlevels;
+        ostringstream param;
+        param << refinement_levels;
+        CCTK_ParameterSet
+          ("refinement_levels", "CarpetRegrid", param.str().c_str());
+        
+      }
     }
     
-    if (cctkGH->cctk_iteration == 0) {
-      newnumlevels = refinement_levels;
-    } else {
-      if (regrid_from_function) {
-        if (CCTK_IsFunctionAliased("RegridLevel")) {
-          int tempnewnumlevels = 0;
-          newnumlevels = 0;
-          BEGIN_MGLEVEL_LOOP(cctkGH) {
-            tempnewnumlevels = 
-              RegridLevel(cctkGH, reflevel, refinement_levels,maxreflevels);
-            newnumlevels = max(tempnewnumlevels, newnumlevels);
-          } END_MGLEVEL_LOOP;
-        } else {
-          CCTK_WARN(1, "No thorn has provided the function " 
-                    "\"RegridLevel\". Regridding will not be done.");
-        }
-      } else {
-        // The standard progressive MR number of levels.
-        newnumlevels = initial_refinement_levels + 
-          activate_newlevels_on_regrid;
+    if (regrid_from_function) {
+      if (! CCTK_IsFunctionAliased("RegridLevel")) {
+        CCTK_WARN (0, "No thorn has provided the function \"RegridLevel\"");
       }
+      const int newnumlevels
+        = RegridLevel (cctkGH, reflevel, refinement_levels, maxreflevels);
+      assert (newnumlevels>0 && newnumlevels<=maxreflevels);
+      
+      *const_cast<CCTK_INT*>(&refinement_levels) = newnumlevels;
+      ostringstream param;
+      param << refinement_levels;
+      CCTK_ParameterSet
+        ("refinement_levels", "CarpetRegrid", param.str().c_str());
     }
-
-    newnumlevels = min(newnumlevels, maxreflevels);
-
-    if ( ( newnumlevels >= 1) && (newnumlevels <= maxreflevels )) {
-      char numlevelstring[10];
-      sprintf(numlevelstring,"%d",newnumlevels);
-      CCTK_ParameterSet("refinement_levels","carpetregrid",numlevelstring);
-     }
     
-    list<ibbox> bbl;
-    list<bvect> obl;
+    
+    
+    int do_recompose;
     
     if (CCTK_EQUALS(refined_regions, "none")) {
       
-      MakeRegions_BaseLevel (cctkGH, bbl, obl);
-      
+      do_recompose = BaseLevel
+        (cctkGH, hh, reflevel, map,
+         size, nboundaryzones, is_internal, is_staggered, shiftout,
+         bbsss, obss, pss);
+                 
     } else if (CCTK_EQUALS(refined_regions, "centre")) {
       
-      MakeRegions_RefineCentre (cctkGH, newnumlevels, bbl, obl);
-      
+      do_recompose = Centre
+        (cctkGH, hh, reflevel, map,
+         size, nboundaryzones, is_internal, is_staggered, shiftout,
+         bbsss, obss, pss);
+                 
     } else if (CCTK_EQUALS(refined_regions, "manual-gridpoints")) {
       
-      if (refinement_levels > 4) {
-	CCTK_WARN (0, "Cannot currently specify manual refinement regions for more than 4 refinement levels");
-      }
-      assert (refinement_levels<=4);
-      vector<ivect> lower(3), upper(3);
-      lower[0] = ivect (l1ixmin, l1iymin, l1izmin);
-      upper[0] = ivect (l1ixmax, l1iymax, l1izmax);
-      lower[1] = ivect (l2ixmin, l2iymin, l2izmin);
-      upper[1] = ivect (l2ixmax, l2iymax, l2izmax);
-      lower[2] = ivect (l3ixmin, l3iymin, l3izmin);
-      upper[2] = ivect (l3ixmax, l3iymax, l3izmax);
-      MakeRegions_AsSpecified (cctkGH, newnumlevels, lower, upper,
-			       bbl, obl);
-      
+      do_recompose = ManualGridpoints
+        (cctkGH, hh, reflevel, map,
+         size, nboundaryzones, is_internal, is_staggered, shiftout,
+         bbsss, obss, pss);
+                 
     } else if (CCTK_EQUALS(refined_regions, "manual-coordinates")) {
       
-      if (refinement_levels > 4) {
-	CCTK_WARN (0, "Cannot currently specify manual refinement regions for more than 4 refinement levels");
-      }
-      assert (refinement_levels<=4);
-      vector<rvect> lower(3), upper(3);
-      lower[0] = rvect (l1xmin, l1ymin, l1zmin);
-      upper[0] = rvect (l1xmax, l1ymax, l1zmax);
-      lower[1] = rvect (l2xmin, l2ymin, l2zmin);
-      upper[1] = rvect (l2xmax, l2ymax, l2zmax);
-      lower[2] = rvect (l3xmin, l3ymin, l3zmin);
-      upper[2] = rvect (l3xmax, l3ymax, l3zmax);
-      MakeRegions_AsSpecified (cctkGH, newnumlevels, lower, upper,
-			       bbl, obl);
-      
+      do_recompose = ManualCoordinates
+        (cctkGH, hh, reflevel, map,
+         size, nboundaryzones, is_internal, is_staggered, shiftout,
+         bbsss, obss, pss);
+                 
     } else if (CCTK_EQUALS(refined_regions, "manual-gridpoint-list")) {
       
-      vector<vector<ibbox> > bbss;
-      if (strcmp(gridpoints, "") != 0) {
-	istringstream gp_str(gridpoints);
-	gp_str >> bbss;
-      }
-      
-      vector<vector<vect<vect<bool,2>,dim> > > obss;
-      if (strcmp(outerbounds, "") !=0 ) {
-	istringstream ob_str (outerbounds);
-	ob_str >> obss;
-        bool good = obss.size() == bbss.size();
-        if (good) {
-          for (int rl=0; rl<(int)obss.size(); ++rl) {
-            good = good && obss[rl].size() == bbss[rl].size();
-          }
-        }
-        if (! good) {
-          cout << "gridpoints: " << bbss << endl;
-          cout << "outerbounds: " << obss << endl;
-          CCTK_WARN (0, "The parameters \"outerbounds\" and \"gridpoints\" must have the same structure");
-        }
-      } else {
-	obss.resize(bbss.size());
-	for (int rl=0; rl<(int)obss.size(); ++rl) {
-	  obss[rl].resize(bbss[rl].size());
-	  for (int c=0; c<(int)obss[rl].size(); ++c) {
-	    obss[rl][c] = vect<vect<bool,2>,dim>(vect<bool,2>(false));
-	  }
-	}
-      }
-      
-      MakeRegions_AsSpecified (cctkGH, newnumlevels, bbss, obss,
-			       bbl, obl);
-      
+      do_recompose = ManualGridpointList
+        (cctkGH, hh, reflevel, map,
+         size, nboundaryzones, is_internal, is_staggered, shiftout,
+         bbsss, obss, pss);
+                 
     } else if (CCTK_EQUALS(refined_regions, "manual-coordinate-list")) {
       
-      vector<vector<rbbox> > bbss;
-      if (strcmp(coordinates, "") != 0) {
-        istringstream co_str(coordinates);
-        co_str >> bbss;
-      }
-      
-      vector<vector<vect<vect<bool,2>,dim> > > obss;
-      if (strcmp(outerbounds, "") !=0 ) {
-	istringstream ob_str (outerbounds);
-	ob_str >> obss;
-        bool good = obss.size() == bbss.size();
-        if (good) {
-          for (int rl=0; rl<(int)obss.size(); ++rl) {
-            good = good && obss[rl].size() == bbss[rl].size();
-          }
-        }
-        if (! good) {
-          cout << "coordinates: " << bbss << endl;
-          cout << "outerbounds: " << obss << endl;
-          CCTK_WARN (0, "The parameters \"outerbounds\" and \"coordinates\" must have the same structure");
-        }
-      } else {
-	obss.resize(bbss.size());
-	for (int rl=0; rl<(int)obss.size(); ++rl) {
-	  obss[rl].resize(bbss[rl].size());
-	  for (int c=0; c<(int)obss[rl].size(); ++c) {
-	    obss[rl][c] = vect<vect<bool,2>,dim>(vect<bool,2>(false));
-	  }
-	}
-      }
-      
-      MakeRegions_AsSpecified (cctkGH, newnumlevels, bbss, obss,
-			       bbl, obl);
-      
+      do_recompose = ManualCoordinateList
+        (cctkGH, hh, reflevel, map,
+         size, nboundaryzones, is_internal, is_staggered, shiftout,
+         bbsss, obss, pss);
+                 
     } else if (CCTK_EQUALS(refined_regions, "automatic")) {
       
-      const int vi = CCTK_VarIndex (errorvar);
-      assert (vi>=0 && vi<CCTK_NumVars());
-      const int gi = CCTK_GroupIndexFromVarI (vi);
-      assert (gi>=0 && gi<CCTK_NumGroups());
-      const int v1 = CCTK_FirstVarIndexI(gi);
-      assert (v1>=0 && v1<=vi && v1<CCTK_NumVars());
-      
-      assert (CCTK_GroupTypeI(gi) == CCTK_GF);
-      assert (CCTK_VarTypeI(vi) == CCTK_VARIABLE_REAL);
-      assert (CCTK_GroupDimI(gi) == dim);
-      
-      assert (gi < (int)arrdata.size());
-      assert (vi-v1 < (int)arrdata[gi].data.size());
-      assert (arrdata[gi].data[vi-v1]);
-      const gf<CCTK_REAL,dim>& error
-	= *dynamic_cast<const gf<CCTK_REAL,dim>*>(arrdata[gi].data[vi-v1]);
-      
-      MakeRegions_Adaptively (cctkGH, minwidth, minfraction, maxerror, error,
-			      bbl, obl);
-      
+      do_recompose = Automatic
+        (cctkGH, hh, reflevel, map,
+         size, nboundaryzones, is_internal, is_staggered, shiftout,
+         bbsss, obss, pss);
+                 
     } else {
       assert (0);
     }
     
-#ifdef AUTOMATIC_BOUNDARIES
-    assert (obl.size() == 0);
-#else
-    assert (bbl.size() == obl.size());
-#endif
-    
-    // transform bbox list into bbox vector
-    vector<ibbox> bbs;
-    bbs.reserve (bbl.size());
-    for (list<ibbox>::const_iterator it = bbl.begin();
-	 it != bbl.end();
-	 ++it) {
-      bbs.push_back (*it);
-    }
-#ifdef AUTOMATIC_BOUNDARIES
-    vector<bvect> obs;
-    obs.resize (bbs.size());
-    for (int c=0; c<(int)bbs.size(); ++c) {
-      ivect extlo = hh->baseextent.lower();
-      ivect extup = hh->baseextent.upper();
-      ivect str = bbs[c].stride();
-      for (int d=0; d<dim; ++d) {
-	// Decide what is an outer boundary.
-	// Allow it to be one grid point to the interior.
-	assert (bbs[c].lower()[d] >= extlo[d]);
-	obs[c][d][0] = bbs[c].lower()[d] <= extlo[d] + str[d];
-	assert (bbs[c].upper()[d] <= extup[d]);
-	obs[c][d][1] = bbs[c].upper()[d] >= extup[d] - str[d];
-      }
-    }
-#else
-    vector<bvect> obs;
-    obs.reserve (obl.size());
-    for (list<bvect>::const_iterator it = obl.begin();
-	 it != obl.end();
-	 ++it) {
-      obs.push_back (*it);
-    }
-#endif
-    
-    
-    
-    // TODO: ensure nesting properties
-    
-    // make multiprocessor aware
-    SplitRegions (cctkGH, bbs, obs);
-    
-    // make multigrid aware
-    gh<dim>::cexts bbss;
-    bbss = hh->make_reflevel_multigrid_boxes(bbs, mglevels);
-    
-    // insert into grid hierarchy
-    bbsss = hh->extents;
-    obss = hh->outer_boundaries;
-    if (bbss.size() == 0) {
-      // remove all finer levels
-      // TODO: this might not work
-      bbsss.resize(reflevel+1);
-      obss.resize(reflevel+1);
-    } else {
-      assert (reflevel < (int)bbsss.size());
-      if (reflevel+1 == (int)bbsss.size()) {
-	// add a finer level
-	bbsss.push_back (bbss);
-	obss.push_back (obs);
-      } else {
-	// change a finer level
-	bbsss[reflevel+1] = bbss;
-	obss[reflevel+1] = obs;
-      }
-    }
-    
-    MakeProcessors (cctkGH, bbsss, pss);
-    
-    return 1;			// do recompose
-  }
-  
-  
-  
-  void MakeRegions_BaseLevel (const cGH* cctkGH,
-			      list<ibbox>& bbl, list<bvect>& obl)
-  {
-    assert (bbl.empty());
-    assert (obl.empty());
-  }
-  
-  
-  
-  // This is a helpful helper routine.  The user can use it to define
-  // how the hierarchy should be refined.  But the result of this
-  // routine is rather arbitrary.
-  void MakeRegions_RefineCentre (const cGH* cctkGH, const int reflevels,
-				 list<ibbox>& bbl, list<bvect>& obl)
-  {
-    assert (bbl.empty());
-    assert (obl.empty());
-    
-    if (reflevel+1 >= reflevels) return;
-      
-    ivect rstr = hh->baseextent.stride();
-    ivect rlb  = hh->baseextent.lower();
-    ivect rub  = hh->baseextent.upper();
-    
-    for (int rl=0; rl<reflevel+1; ++rl) {
-      // save old values
-      const ivect oldrlb = rlb;
-      const ivect oldrub = rub;
-      // refined boxes have smaller stride
-      assert (all(rstr%hh->reffact == 0));
-      rstr /= hh->reffact;
-      // calculate new extent
-      const ivect quarter = (rub - rlb) / 4 / rstr * rstr;
-      rlb = oldrlb + quarter;
-      rub = oldrub - quarter;
-      assert (all(rlb >= oldrlb && rub <= oldrub));
-    }
-    
-    bbl.push_back (ibbox(rlb, rub, rstr));
-    obl.push_back (bvect(vect<bool,2>(false)));
-  }
-  
-  
-  
-  static void
-  MakeRegions_AsSpecified_OneLevel (const cGH* cctkGH, const int reflevels,
-				    const ivect ilower,
-				    const ivect iupper,
-				    const bvect obound,
-				    list<ibbox>& bbl, list<bvect>& obl)
-  {
-    if (reflevel+1 >= reflevels) return;
-    
-    const ivect rstr = hh->baseextent.stride();
-    const ivect rlb  = hh->baseextent.lower();
-    const ivect rub  = hh->baseextent.upper();
-    
-    const int rl = reflevel+1;
-    
-    const int levfac = ipow(hh->reffact, rl);
-    assert (all (rstr % levfac == 0));
-    const ivect str (rstr / levfac);
-    const ivect lb  (ilower);
-    const ivect ub  (iupper);
-    if (! all(lb>=rlb && ub<=rub)) {
-      ostringstream buf;
-      buf << "The refinement region boundaries for refinement level #" << rl << " are not within the main grid.  Allowed are the grid point boundaries " << rlb << " - " << rub << "; specified were " << lb << " - " << ub << ends;
-      CCTK_WARN (0, buf.str().c_str());
-    }
-    if (! all(lb<=ub)) {
-      CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING,
-		  "The refinement region boundaries for refinement level #%d have the upper boundary less than the lower boundary", rl);
-    }
-    if (! all(lb%str==0 && ub%str==0)) {
-      CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING,
-		  "The refinement region boundaries for refinement level #%d are not a multiple of the stride for that level", rl);
-    }
-    assert (all(lb>=rlb && ub<=rub));
-    assert (all(lb<=ub));
-    assert (all(lb%str==0 && ub%str==0));
-    
-    bbl.push_back (ibbox(lb, ub, str));
-    obl.push_back (obound);
-  }
-  
-  
-  
-  void MakeRegions_AsSpecified (const cGH* cctkGH, const int reflevels,
-				const vector<ivect> lower,
-				const vector<ivect> upper,
-				list<ibbox>& bbl, list<bvect>& obl)
-  {
-    assert (lower.size() == upper.size());
-    
-    if (reflevel+1 >= reflevels) return;
-    
-    const int rl = reflevel+1;
-    
-    const ivect ilower = lower[rl-1];
-    const ivect iupper = upper[rl-1];
-    const bvect obound = bvect(vect<bool,2>(false));
-    
-    MakeRegions_AsSpecified_OneLevel (cctkGH, reflevels,
-				      ilower, iupper, obound,
-				      bbl, obl);
-  }
-  
-  
-  
-  static ivect delta2int (const cGH* cctkGH, const rvect & rpos, const int rl)
-  {
-    rvect global_lower, global_upper;
-    for (int d=0; d<dim; ++d) {
-      const int ierr = CCTK_CoordRange
-	(cctkGH, &global_lower[d], &global_upper[d], d+1, 0, "cart3d");
-      if (ierr<0) {
-	global_lower[d] = 0;
-	global_upper[d] = 1;
-      }
-    }
-    const ivect global_extent (hh->baseextent.upper() - hh->baseextent.lower());
-    
-    CCTK_REAL (* const rfloor) (CCTK_REAL const) = floor;
-    
-    const rvect scale  = rvect(global_extent) / (global_upper - global_lower);
-    const int levfac = ipow(hh->reffact, rl);
-    assert (all (hh->baseextent.stride() % levfac == 0));
-    const ivect istride = hh->baseextent.stride() / levfac;
-    
-    const ivect ipos = ivect(map(rfloor, rpos * scale / rvect(istride) + 0.5)) * istride;
-    
-    const rvect apos = rpos * scale;
-    assert (all(abs(apos - rvect(ipos)) < rvect(istride)*0.01));
-    
-    return ipos;
-  }
-  
-  
-  
-  static ivect pos2int (const cGH* cctkGH, const rvect & rpos, const int rl)
-  {
-    rvect global_lower, global_upper;
-    for (int d=0; d<dim; ++d) {
-      const int ierr = CCTK_CoordRange
-	(cctkGH, &global_lower[d], &global_upper[d], d+1, 0, "cart3d");
-      if (ierr<0) {
-	global_lower[d] = 0;
-	global_upper[d] = 1;
-      }
-    }
-    const ivect global_extent (hh->baseextent.upper() - hh->baseextent.lower());
-    
-    CCTK_REAL (* const rfloor) (CCTK_REAL const) = floor;
-    
-    const rvect scale  = rvect(global_extent) / (global_upper - global_lower);
-    const int levfac = ipow(hh->reffact, rl);
-    assert (all (hh->baseextent.stride() % levfac == 0));
-    const ivect istride = hh->baseextent.stride() / levfac;
-    
-    const ivect ipos = ivect(map(rfloor, (rpos - global_lower) * scale / rvect(istride) + 0.5)) * istride;
-    
-    return ipos;
-  }
-  
-  
-  
-#if 0
-  void MakeRegions_AsSpecified (const cGH* cctkGH, const int reflevels,
-				const vector<rvect> lower,
-				const vector<rvect> upper,
-				list<ibbox>& bbl, list<bvect>& obl)
-  {
-    assert (lower.size() == upper.size());
-    
-    if (reflevel+1 >= reflevels) return;
-    
-    rvect global_lower, global_upper;
-    for (int d=0; d<dim; ++d) {
-      const int ierr = CCTK_CoordRange
-	(cctkGH, &global_lower[d], &global_upper[d], d+1, 0, "cart3d");
-      if (ierr<0) {
-	global_lower[d] = 0;
-	global_upper[d] = 1;
-      }
-    }
-    const ivect global_extent (hh->baseextent.upper() - hh->baseextent.lower());
-    
-    const int rl = reflevel+1;
-    
-    const rvect scale  = rvect(global_extent) / (global_upper - global_lower);
-    const rvect rlower = (lower[rl-1] - global_lower) * scale;
-    const rvect rupper = (upper[rl-1] - global_lower) * scale;
-//     const ivect ilower = ivect(map((CCTK_REAL(*)(CCTK_REAL))floor, rlower + 0.5));
-//     const ivect iupper = ivect(map((CCTK_REAL(*)(CCTK_REAL))floor, rupper + 0.5));
-    CCTK_REAL (* const rfloor) (CCTK_REAL const) = floor;
-    const int levfac = ipow(hh->reffact, rl);
-    assert (all (hh->baseextent.stride() % levfac == 0));
-    const ivect istride = hh->baseextent.stride() / levfac;
-    const ivect ilower  = ivect(map(rfloor, rlower / rvect(istride) + 0.5)) * istride;
-    const ivect iupper  = ivect(map(rfloor, rupper / rvect(istride) + 0.5)) * istride;
-    const bvect obound  = bvect(vect<bool,2>(false));
-    
-    MakeRegions_AsSpecified_OneLevel (cctkGH, reflevels,
-				      ilower, iupper, obound,
-				      bbl, obl);
-  }
-#endif
-  
-  
-  
-  void MakeRegions_AsSpecified (const cGH* cctkGH, const int reflevels,
-				const vector<rvect> lower,
-				const vector<rvect> upper,
-				list<ibbox>& bbl, list<bvect>& obl)
-  {
-    assert (lower.size() == upper.size());
-    vector<ivect> ilower, iupper;
-    ilower.resize (lower.size());
-    iupper.resize (upper.size());
-    for (int rl=1; rl<lower.size()+1; ++rl) {
-      ilower[rl-1] = pos2int (cctkGH, lower[rl-1], rl);
-      iupper[rl-1] = pos2int (cctkGH, upper[rl-1], rl);
-    }
-    MakeRegions_AsSpecified (cctkGH, reflevels, ilower, iupper, bbl, obl);
-  }
-  
-  
-  
-  void MakeRegions_AsSpecified (const cGH* cctkGH, const int reflevels,
-				const vector<vector<ibbox> > bbss,
-				const vector<vector<bvect> > obss,
-				list<ibbox>& bbl, list<bvect>& obl)
-  {
-    if (reflevel+1 >= reflevels) return;
-    
-    const int rl = reflevel+1;
-    assert (rl-1 < (int)bbss.size());
-    assert (obss.size() == bbss.size());
-    
-    for (int c=0; c<(int)bbss[rl-1].size(); ++c) {
-      assert (obss[rl-1].size() == bbss[rl-1].size());
-      
-      const ivect ilower = bbss[rl-1][c].lower();
-      const ivect iupper = bbss[rl-1][c].upper();
-      const bvect obound = obss[rl-1][c];
-      
-      MakeRegions_AsSpecified_OneLevel (cctkGH, reflevels,
-					ilower, iupper, obound,
-					bbl, obl);
-      
-    }
-  }
-  
-  
-  
-  void MakeRegions_AsSpecified (const cGH* cctkGH, const int reflevels,
-				const vector<vector<rbbox> > bbss,
-				const vector<vector<bvect> > obss,
-				list<ibbox>& bbl, list<bvect>& obl)
-  {
-    vector<vector<ibbox> > ibbss;
-    ibbss.resize (bbss.size());
-    for (int rl=1; rl<bbss.size()+1; ++rl) {
-      ibbss[rl-1].resize (bbss[rl-1].size());
-      for (int c=0; c<bbss[rl-1].size(); ++c) {
-        const ivect ilower = pos2int (cctkGH, bbss[rl-1][c].lower(), rl);
-        const ivect iupper = pos2int (cctkGH, bbss[rl-1][c].upper(), rl);
-        const ivect istride = delta2int (cctkGH, bbss[rl-1][c].stride(), rl);
-        ibbss[rl-1][c] = ibbox (ilower, iupper, istride);
-      }
-    }
-    MakeRegions_AsSpecified (cctkGH, reflevels, ibbss, obss, bbl, obl);
-  }
-  
-  
-  
-  static void
-  MakeRegions_Adaptively_Recombine (list<ibbox>& bbl1,
-				    list<ibbox>& bbl2,
-				    list<ibbox>& bbl,
-				    const ibbox& iface,
-				    const int dir)
-  {
-    assert (!iface.empty());
-    assert (iface.lower()[dir] == iface.upper()[dir]);
-    
-    const int oldnumboxes = bbl.size() + bbl1.size() + bbl2.size();
-    int numcombinedboxes = 0;
-    
-    int oldnumpoints = 0;
-    for (list<ibbox>::const_iterator ibb = bbl.begin(); ibb != bbl.end(); ++ibb) {
-      oldnumpoints += ibb->size();
-    }
-    for (list<ibbox>::const_iterator ibb1 = bbl1.begin(); ibb1 != bbl1.end(); ++ibb1) {
-      oldnumpoints += ibb1->size();
-    }
-    for (list<ibbox>::const_iterator ibb2 = bbl2.begin(); ibb2 != bbl2.end(); ++ibb2) {
-      oldnumpoints += ibb2->size();
-    }
-    
-#if 0
-    // remember old bounding boxes
-    bboxset<int,dim> oldboxes;
-    for (list<ibbox>::const_iterator ibb1 = bbl1.begin(); ibb1 != bbl1.end(); ++ibb1) {
-      oldboxes += *ibb1;
-    }
-    for (list<ibbox>::const_iterator ibb2 = bbl2.begin(); ibb2 != bbl2.end(); ++ibb2) {
-      oldboxes += *ibb2;
-    }
-#endif
-#if 0
-    cout << endl;
-    cout << "MakeRegions_Adaptively_Recombine: initial list:" << endl;
-    cout << bbl << endl;
-    cout << "MakeRegions_Adaptively_Recombine: initial list 1:" << endl;
-    cout << bbl1 << endl;
-    cout << "MakeRegions_Adaptively_Recombine: initial list 2:" << endl;
-    cout << bbl2 << endl;
-#endif
-    
-    const ivect lo = iface.lower();
-    const ivect up = iface.upper();
-    const ivect str = iface.stride();
-    
-    {
-      // prune boxes on the left
-      list<ibbox>::iterator ibb1 = bbl1.begin();
-      while (ibb1 != bbl1.end()) {
-	// is this bbox just to the left of the interface?
-// 	const ivect lo1 = ibb1->lower();
-	const ivect up1 = ibb1->upper();
-	const ivect str1 = ibb1->stride();
-	assert (up1[dir]+str1[dir] <= lo[dir]);
-	assert (all(str1 == str));
-	if (up1[dir]+str1[dir] < lo[dir]) {
-	  // no: forget it
-	  bbl.push_back (*ibb1);
-	  ibb1 = bbl1.erase(ibb1);
-	  continue;
-	}
-	++ibb1;
-      }	// while
-    }
-    
-    {
-      // prune boxes on the right
-      list<ibbox>::iterator ibb2 = bbl2.begin();
-      while (ibb2 != bbl2.end()) {
-	// is this bbox just to the right of the interface?
-	const ivect lo2 = ibb2->lower();
-// 	const ivect up2 = ibb2->upper();
-	const ivect str2 = ibb2->stride();
-	assert (up[dir] <= lo2[dir]);
-	assert (all(str2 == str));
-	if (up[dir] < lo2[dir]) {
-	  // no: forget it
-	  bbl.push_back (*ibb2);
-	  ibb2 = bbl2.erase(ibb2);
-	  continue;
-	}
-	++ibb2;
-      }	// while
-    }
-    
-    {
-      // walk all boxes on the left
-      list<ibbox>::iterator ibb1 = bbl1.begin();
-      while (ibb1 != bbl1.end()) {
-	ivect lo1 = ibb1->lower();
-	ivect up1 = ibb1->upper();
-	ivect str1 = ibb1->stride();
-	assert (up1[dir]+str1[dir] == lo[dir]);
-	lo1[dir] = lo[dir];
-	up1[dir] = up[dir];
-	const ibbox iface1 (lo1,up1,str1);
-	
-	{
-	  // walk all boxes on the right
-	  list<ibbox>::iterator ibb2 = bbl2.begin();
-	  while (ibb2 != bbl2.end()) {
-	    ivect lo2 = ibb2->lower();
-	    ivect up2 = ibb2->upper();
-	    ivect str2 = ibb2->stride();
-	    assert (lo2[dir] == up[dir]);
-	    lo2[dir] = lo[dir];
-	    up2[dir] = up[dir];
-	    const ibbox iface2 (lo2,up2,str2);
-	    
-	    // check for a match
-	    if (iface1 == iface2) {
-	      const ibbox combined (ibb1->lower(), ibb2->upper(), str);
-	      bbl.push_back (combined);
-	      ibb1 = bbl1.erase(ibb1);
-	      ibb2 = bbl2.erase(ibb2);
-	      ++numcombinedboxes;
-// 	      cout << "MRA: Combining along " << "xyz"[dir] << " to " << bbl.back() << " size " << bbl.back().size() << endl;
-	      goto continue_search;
-	    }
-	    
-	    ++ibb2;
-	  } // while
-	}
-	
-	++ibb1;
-      continue_search:;
-      }	// while
-    }
-    
-    bbl.splice (bbl.end(), bbl1);
-    bbl.splice (bbl.end(), bbl2);
-    
-    assert (bbl1.empty() && bbl2.empty());
-    
-    const int newnumboxes = bbl.size();
-    assert (newnumboxes + numcombinedboxes == oldnumboxes);
-    
-    int newnumpoints = 0;
-    for (list<ibbox>::const_iterator ibb = bbl.begin(); ibb != bbl.end(); ++ibb) {
-      newnumpoints += ibb->size();
-    }
-    assert (newnumpoints == oldnumpoints);
-    
-#if 0
-    // find new bounding boxes
-    bboxset<int,dim> newboxes;
-    for (list<ibbox>::const_iterator ibb = bbl.begin(); ibb != bbl.end(); ++ibb) {
-      newboxes += *ibb;
-    }
-    
-    // Check that they are equal
-    assert (newboxes.size() <= oldboxes.size());
-    assert ((newboxes.size()==0) == (oldboxes.size()==0));
-    assert (oldboxes == newboxes);
-#endif
-#if 0
-    cout << "MakeRegions_Adaptively_Recombine: final list:" << endl;
-    cout << bbl << endl;
-    cout << endl;
-#endif
-  }
-  
-  static void
-  MakeRegions_Adaptively (const cGH* const cctkGH,
-			  const int minwidth,
-			  const CCTK_REAL minfraction,
-			  const CCTK_REAL maxerror,
-			  const data<CCTK_REAL,dim>& error,
-			  list<ibbox>& bbl,
-			  const ibbox& region)
-  {
-    // Just to be sure
-    assert (! region.empty());
-    
-    // Count grid points that need to be refined
-    // (this doesn't work yet on multiple processors)
-    assert (CCTK_nProcs(cctkGH)==1);
-    int cnt = 0;
-    for (ibbox::iterator it=region.begin(); it!=region.end(); ++it) {
-      if (error[*it] > maxerror) ++cnt;
-    }
-    const CCTK_REAL fraction = (CCTK_REAL)cnt / region.size();
-    const int width = maxval(region.shape() / region.stride());
-    
-    if (cnt == 0) {
-      // Don't refine
-    } else if (width < 2*minwidth || fraction >= minfraction) {
-      // Refine the whole region
-      const ivect lo(region.lower());
-      const ivect up(region.upper());
-      const ivect str(region.stride());
-      bbl.push_back (ibbox(lo,up+str-str/reffact,str/reffact));
-//       cout << "MRA: Refining to " << bbl.back() << " size " << bbl.back().size() << " fraction " << fraction << endl;
-    } else {
-      // Split the region and check recursively
-      const int dir = maxloc(region.shape());
-      const ivect lo(region.lower());
-      const ivect up(region.upper());
-      const ivect str(region.stride());
-      ivect lo1(lo), lo2(lo);
-      ivect up1(up), up2(up);
-      const int mgstr = ipow(hh->mgfact, mglevels); // honour multigrid factors
-      const int step = str[dir]*mgstr;
-      lo2[dir] = ((up[dir]+lo[dir])/2/step)*step;
-      up1[dir] = lo2[dir]-str[dir];
-      const ibbox region1(lo1,up1,str);
-      const ibbox region2(lo2,up2,str);
-      assert (region1.is_contained_in(region));
-      assert (region2.is_contained_in(region));
-      assert ((region1 & region2).empty());
-      assert (region1 + region2 == region);
-      list<ibbox> bbl1, bbl2;
-      MakeRegions_Adaptively (cctkGH, minwidth, minfraction, maxerror,
-			      error, bbl1, region1);
-      MakeRegions_Adaptively (cctkGH, minwidth, minfraction, maxerror,
-			      error, bbl2, region2);
-      // Combine regions if possible
-      up2 += str-str/reffact;
-      up2[dir] = lo2[dir];
-      const ibbox iface(lo2,up2,str/reffact);
-      MakeRegions_Adaptively_Recombine (bbl1, bbl2, bbl, iface, dir);
-    }
-    
-  }
-  
-  void
-  MakeRegions_Adaptively (const cGH* const cctkGH,
-			  const int minwidth,
-			  const CCTK_REAL minfraction,
-			  const CCTK_REAL maxerror,
-			  const gf<CCTK_REAL,dim>& error,
-			  list<ibbox>& bbl, list<bvect>& obl)
-  {
-    assert (bbl.empty());
-    
-    if (reflevel+1 >= maxreflevels) return;
-    
-    assert (component == -1);
-    
-    const int rl = reflevel;
-    
-    // Arbitrary
-    const int tl = 0;
-    const int ml = 0;
-    
-//     cout << endl << "MRA: Choosing regions to refine in " << hh->components(rl) << " components" << endl;
-    
-    for (int c=0; c<hh->components(rl); ++c) {
-      const ibbox region = hh->extents[rl][c][ml];
-      assert (! region.empty());
-      
-      const data<CCTK_REAL,dim>& errdata = *error(tl,rl,c,ml);
-      
-      MakeRegions_Adaptively (cctkGH, minwidth, minfraction, maxerror,
-			      errdata, bbl, region);
-    }
-    
-//     int numpoints = 0;
-//     for (list<ibbox>::const_iterator ibb = bbl.begin(); ibb != bbl.end(); ++ibb) {
-//       numpoints += ibb->size();
-//     }
-//     cout << "MRA: Chose " << bbl.size() << " regions with a total size of " << numpoints << " to refine." << endl << endl;
-    
-    // Remove grid points outside the outer boundary
-    vect<bool,dim> obp;
-    {
-      DECLARE_CCTK_PARAMETERS;
-      assert (sizeof outside_boundary_points == dim * sizeof (CCTK_INT));
-      obp = outside_boundary_points;
-    }
-    for (list<ibbox>::iterator it = bbl.begin();
-	 it != bbl.end();
-	 ++it) {
-      const ivect ub = obp.ifthen (it->upper(),
-				   min (it->upper(), hh->baseextent.upper()));
-      *it = ibbox(it->lower(), ub, it->stride());
-    }
-    
-    // Create obl from bbl
-    for (list<ibbox>::const_iterator it = bbl.begin();
-	 it != bbl.end();
-	 ++it) {
-      obl.push_back (zip ((vect<bool,2> (*) (bool, bool)) &vect<bool,2>::make,
-			  it->lower() == hh->baseextent.lower(),
-			  it->upper() == hh->baseextent.upper()));
-    }
-    
+    return do_recompose;
   }
   
 } // namespace CarpetRegrid
diff --git a/Carpet/CarpetRegrid/src/regrid.hh b/Carpet/CarpetRegrid/src/regrid.hh
index 45a390a97..27edc99c9 100644
--- a/Carpet/CarpetRegrid/src/regrid.hh
+++ b/Carpet/CarpetRegrid/src/regrid.hh
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/src/regrid.hh,v 1.10 2003/11/13 16:04:37 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/src/regrid.hh,v 1.11 2004/01/25 14:57:30 schnetter Exp $
 
 #ifndef CARPETREGRID_HH
 #define CARPETREGRID_HH
@@ -6,14 +6,15 @@
 #include <list>
 
 #include "cctk.h"
+#include "cctk_Arguments.h"
 
+#include "bbox.hh"
 #include "gf.hh"
 #include "gh.hh"
+#include "vect.hh"
 
 #include "carpet.hh"
 
-#include "regrid.h"
-
 
 
 namespace CarpetRegrid {
@@ -23,53 +24,171 @@ namespace CarpetRegrid {
   
   
   
-  typedef vect<int,dim> ivect;
-  typedef bbox<int,dim> ibbox;
-  
-  typedef vect<vect<bool,2>,dim> bvect;
-  
-  typedef vect<CCTK_REAL,dim> rvect;
-  typedef bbox<CCTK_REAL,dim> rbbox;
-  
-  
-  
-  int CarpetRegridRegrid (const cGH * const cctkGH,
-			  gh<dim>::rexts& bbsss,
-			  gh<dim>::rbnds& obss,
-			  gh<dim>::rprocs& pss);
-  
-  
-  
-  void MakeRegions_BaseLevel    (const cGH* cctkGH,
-				 list<ibbox>& bbl, list<bvect>& obl);
-  
-  void MakeRegions_RefineCentre (const cGH* cctkGH, const int reflevels,
-				 list<ibbox>& bbl, list<bvect>& obl);
-  
-  void MakeRegions_AsSpecified  (const cGH* cctkGH, const int reflevels,
-				 const vector<ivect> lower,
-				 const vector<ivect> upper,
-				 list<ibbox>& bbl, list<bvect>& obl);
-  void MakeRegions_AsSpecified  (const cGH* cctkGH, const int reflevels,
-				 const vector<rvect> lower,
-				 const vector<rvect> upper,
-				 list<ibbox>& bbl, list<bvect>& obl);
-  
-  void MakeRegions_AsSpecified  (const cGH* cctkGH, const int reflevels,
-				 const vector<vector<ibbox> > bbss,
-				 const vector<vector<bvect> > obss,
-				 list<ibbox>& bbl, list<bvect>& obl);
-  void MakeRegions_AsSpecified  (const cGH* cctkGH, const int reflevels,
-				 const vector<vector<rbbox> > bbss,
-				 const vector<vector<bvect> > obss,
-				 list<ibbox>& bbl, list<bvect>& obl);
-  
-  void MakeRegions_Adaptively   (const cGH* cctkGH,
-				 const int minwidth,
-				 const CCTK_REAL minfraction,
-				 const CCTK_REAL maxerror,
-				 const gf<CCTK_REAL,dim>& error,
-				 list<ibbox>& bbl, list<bvect>& obl);
+  extern "C" {
+    
+    /* Scheduled functions */
+    int CarpetRegridParamcheck (CCTK_ARGUMENTS);
+    
+    /* Aliased functions */
+//     CCTK_INT CarpetRegrid_Regrid (const cGH * const cctkGH,
+//                                   gh<dim>::rexts  * bbsss,
+//                                   gh<dim>::rbnds  * obss,
+//                                   gh<dim>::rprocs * pss);
+    CCTK_INT CarpetRegrid_Regrid (CCTK_POINTER_TO_CONST const cctkGH_,
+                                  CCTK_INT const reflevel,
+                                  CCTK_INT const map,
+                                  CCTK_INT const size,
+                                  CCTK_INT const * const nboundaryzones,
+                                  CCTK_INT const * const is_internal,
+                                  CCTK_INT const * const is_staggered,
+                                  CCTK_INT const * const shiftout,
+                                  CCTK_POINTER const bbsss_,
+                                  CCTK_POINTER const obss_,
+                                  CCTK_POINTER const pss_);
+  }
+  
+  
+  
+  int BaseLevel (cGH const * const cctkGH,
+                 gh<dim> const & hh,
+                 int const reflevel,
+                 int const map,
+                 int const size,
+                 jjvect const & nboundaryzones,
+                 jjvect const & is_internal,
+                 jjvect const & is_staggered,
+                 jjvect const & shiftout,
+                 gh<dim>::rexts  & bbsss,
+                 gh<dim>::rbnds  & obss,
+                 gh<dim>::rprocs & pss);
+
+  int Centre (cGH const * const cctkGH,
+              gh<dim> const & hh,
+              int const reflevel,
+              int const map,
+              int const size,
+              jjvect const & nboundaryzones,
+              jjvect const & is_internal,
+              jjvect const & is_staggered,
+              jjvect const & shiftout,
+              gh<dim>::rexts  & bbsss,
+              gh<dim>::rbnds  & obss,
+              gh<dim>::rprocs & pss);
+
+  int ManualGridpoints (cGH const * const cctkGH,
+                        gh<dim> const & hh,
+                        int const reflevel,
+                        int const map,
+                        int const size,
+                        jjvect const & nboundaryzones,
+                        jjvect const & is_internal,
+                        jjvect const & is_staggered,
+                        jjvect const & shiftout,
+                        gh<dim>::rexts  & bbsss,
+                        gh<dim>::rbnds  & obss,
+                        gh<dim>::rprocs & pss);
+  
+  void ManualGridpoints_OneLevel (const cGH * const cctkGH,
+                                  const gh<dim> & hh,
+                                  const int reflevel,
+                                  const int reflevels,
+                                  const ivect ilower,
+                                  const ivect iupper,
+                                  const bbvect obound,
+                                  vector<ibbox> & bbs,
+                                  vector<bbvect> & obs);
+
+  int ManualCoordinates (cGH const * const cctkGH,
+                         gh<dim> const & hh,
+                         int const reflevel,
+                         int const map,
+                         int const size,
+                         jjvect const & nboundaryzones,
+                         jjvect const & is_internal,
+                         jjvect const & is_staggered,
+                         jjvect const & shiftout,
+                         gh<dim>::rexts  & bbsss,
+                         gh<dim>::rbnds  & obss,
+                         gh<dim>::rprocs & pss);
+  
+  void ManualCoordinates_OneLevel (const cGH * const cctkGH,
+                                   const gh<dim> & hh,
+                                   const int reflevel,
+                                   const int reflevels,
+                                   const rvect ilower,
+                                   const rvect iupper,
+                                   const bbvect obound,
+                                   vector<ibbox> & bbs,
+                                   vector<bbvect> & obs);
+  
+  ivect delta2int (const cGH * const cctkGH, const gh<dim>& hh,
+                   const rvect & rpos, const int reflevel);
+  ivect pos2int (const cGH* const cctkGH, const gh<dim>& hh,
+                 const rvect & rpos, const int reflevel);
+
+  int ManualGridpointList (cGH const * const cctkGH,
+                           gh<dim> const & hh,
+                           int const reflevel,
+                           int const map,
+                           int const size,
+                           jjvect const & nboundaryzones,
+                           jjvect const & is_internal,
+                           jjvect const & is_staggered,
+                           jjvect const & shiftout,
+                           gh<dim>::rexts  & bbsss,
+                           gh<dim>::rbnds  & obss,
+                           gh<dim>::rprocs & pss);
+  
+  int ManualCoordinateList (cGH const * const cctkGH,
+                            gh<dim> const & hh,
+                            int const reflevel,
+                            int const map,
+                            int const size,
+                            jjvect const & nboundaryzones,
+                            jjvect const & is_internal,
+                            jjvect const & is_staggered,
+                            jjvect const & shiftout,
+                            gh<dim>::rexts  & bbsss,
+                            gh<dim>::rbnds  & obss,
+                            gh<dim>::rprocs & pss);
+  
+  int Automatic (cGH const * const cctkGH,
+                 gh<dim> const & hh,
+                 int const reflevel,
+                 int const map,
+                 int const size,
+                 jjvect const & nboundaryzones,
+                 jjvect const & is_internal,
+                 jjvect const & is_staggered,
+                 jjvect const & shiftout,
+                 gh<dim>::rexts  & bbsss,
+                 gh<dim>::rbnds  & obss,
+                 gh<dim>::rprocs & pss);
+  
+  void Automatic_OneLevel (const cGH * const cctkGH,
+                           const gh<dim> & hh,
+                           const int reflevel,
+                           const int minwidth,
+                           const CCTK_REAL minfraction,
+                           const CCTK_REAL maxerror,
+                           const gf<CCTK_REAL,dim> & errorvar,
+                           vector<ibbox> & bbs,
+                           vector<bbvect> & obs);
+  
+  void Automatic_Recursive (const cGH * const cctkGH,
+                            const gh<dim> & hh,
+                            const int minwidth,
+                            const CCTK_REAL minfraction,
+                            const CCTK_REAL maxerror,
+                            const data<CCTK_REAL,dim> & errorvar,
+                            list<ibbox> & bbl,
+                            const ibbox & region);
+  
+  void Automatic_Recombine (list<ibbox> & bbl1,
+                            list<ibbox> & bbl2,
+                            list<ibbox> & bbl,
+                            const ibbox & iface,
+                            const int dir);
   
 } // namespace CarpetRegrid