Added stream input routines for some CarpetLib containers.

Added stream input routines for some CarpetLib containers.
The regridder now has to explicitly say which boundaries are outer,
and which are internal.  This will make outer boundaries on fine grid
possible, and is also necessary when there are multiple grid patches.

Started to add support for arbitrariliy many user-specified refinement
regions.  Not yet finished.

The Carpet driver can now handle multiple grid patches.

Added example files for multiple grid patches.  They use initial data
that does not "fit" the boundary conditions, and they don't use
multiple refinement levels so far.

Removed old and unused example files in CarpetLib.

darcs-hash:20020311121709-07bb3-18594c42bd7a958ee0840d29e158a343208f5711.gz

23 files changed, 563 insertions, 244 deletions

diff --git a/Carpet/Carpet/param.ccl b/Carpet/Carpet/param.ccl
index 0d8456037..962f84d8e 100644
--- a/Carpet/Carpet/param.ccl
+++ b/Carpet/Carpet/param.ccl
@@ -1,5 +1,5 @@
 # Parameter definitions for thorn Carpet
-# $Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/param.ccl,v 1.17 2002/01/11 17:58:47 schnetter Exp $
+# $Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/param.ccl,v 1.18 2002/03/11 13:17:09 schnetter Exp $
 
 shares: Cactus
 
@@ -98,6 +98,40 @@ CCTK_INT prolongation_order_time "Order of prolongation operator in time"
 
 
 
+CCTK_STRING base_extents "Extents of base grid components, in grid point units of the finest level"
+{
+  "^$" :: "leave empty for one grid component covering the whole region (default)"
+# We want the string to contain a list of bboxes.  Each bbox contains
+# three vectors specifying the lower bound, upper bound, and stride.
+# (The upper bound is inclusive.  All values are nonnegative integers.)
+# The syntax for vectors, bboxes, and lists is described below.
+# All spaces are optional.
+# The () parentheses and [] brackets are literals.
+# The {} braces are used for grouping.
+# If you put everything together, you get the tapeworm below.
+# There should be a way to define abbreviations in regexps!
+# Having #defines in param.ccl files would help a lot.
+# VECT := " [ 0 , 0 , 0 ]"
+# BBOX := " (VECT :VECT :VECT )"
+# LIST := " [{{BBOX ,}*BBOX}? ]"
+  "^\[(([[:space:]]*\([[:space:]]*\[[[:space:]]*[[:digit:]]+[[:space:]]*,[[:space:]]*[[:digit:]]+[[:space:]]*,[[:space:]]*[[:digit:]]+[[:space:]]*\][[:space:]]*:[[:space:]]*\[[[:space:]]*[[:digit:]]+[[:space:]]*,[[:space:]]*[[:digit:]]+[[:space:]]*,[[:space:]]*[[:digit:]]+[[:space:]]*\][[:space:]]*:[[:space:]]*\[[[:space:]]*[[:digit:]]+[[:space:]]*,[[:space:]]*[[:digit:]]+[[:space:]]*,[[:space:]]*[[:digit:]]+[[:space:]]*\][[:space:]]*\)[[:space:]]*,)*[[:space:]]*\([[:space:]]*\[[[:space:]]*[[:digit:]]+[[:space:]]*,[[:space:]]*[[:digit:]]+[[:space:]]*,[[:space:]]*[[:digit:]]+[[:space:]]*\][[:space:]]*:[[:space:]]*\[[[:space:]]*[[:digit:]]+[[:space:]]*,[[:space:]]*[[:digit:]]+[[:space:]]*,[[:space:]]*[[:digit:]]+[[:space:]]*\][[:space:]]*:[[:space:]]*\[[[:space:]]*[[:digit:]]+[[:space:]]*,[[:space:]]*[[:digit:]]+[[:space:]]*,[[:space:]]*[[:digit:]]+[[:space:]]*\][[:space:]]*\))?[[:space:]]*\]$" :: "[ ([<imin>,<jmin>,<kmin>]:[<imax>,<jmax>,<kmax>]:[<istride>,<jstride>,<kstride>]), ... ]"
+} ""
+
+CCTK_STRING base_bboxes "Bounding boxes of base grid components"
+{
+  "^$" :: "leave empty for using the default, which depends on cctk_gsh"
+# See above for an explanation of this syntax, and of the tapeworm below.
+# Each vector element is 0 or 1,
+# where 0 is handled by synchronisation or prolongation,
+# and 1 stands for a user-supplied ("outer") boundary condition.
+# BND := " [ 0 , 0 ]"
+# VECT := " [ BND , BND , BND ]"
+# LIST := " [{{VECT ,}*VECT}? ]"
+  "^\[(([[:space:]]*\[[[:space:]]*\[[[:space:]]*[[:digit:]]+[[:space:]]*,[[:space:]]*[[:digit:]]+[[:space:]]*\][[:space:]]*,[[:space:]]*\[[[:space:]]*[[:digit:]]+[[:space:]]*,[[:space:]]*[[:digit:]]+[[:space:]]*\][[:space:]]*,[[:space:]]*\[[[:space:]]*[[:digit:]]+[[:space:]]*,[[:space:]]*[[:digit:]]+[[:space:]]*\][[:space:]]*\][[:space:]]*,)*[[:space:]]*\[[[:space:]]*\[[[:space:]]*[[:digit:]]+[[:space:]]*,[[:space:]]*[[:digit:]]+[[:space:]]*\][[:space:]]*,[[:space:]]*\[[[:space:]]*[[:digit:]]+[[:space:]]*,[[:space:]]*[[:digit:]]+[[:space:]]*\][[:space:]]*,[[:space:]]*\[[[:space:]]*[[:digit:]]+[[:space:]]*,[[:space:]]*[[:digit:]]+[[:space:]]*\][[:space:]]*\])?[[:space:]]*\]$" :: "[ [ [<ilower>,<iupper>], [<jlower>,<jupper>], [<klower>,<kupper>] ], ... ]"
+} ""
+
+
+
 BOOLEAN enable_all_storage "Enable storage for all grid functions"
 {
 } "no"
diff --git a/Carpet/Carpet/src/Recompose.cc b/Carpet/Carpet/src/Recompose.cc
index 12277aa99..ec3c65047 100644
--- a/Carpet/Carpet/src/Recompose.cc
+++ b/Carpet/Carpet/src/Recompose.cc
@@ -21,7 +21,7 @@
 
 #include "carpet.hh"
 
-static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Recompose.cc,v 1.20 2002/01/30 16:07:58 schnetter Exp $";
+static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Recompose.cc,v 1.21 2002/03/11 13:17:10 schnetter Exp $";
 
 
 
@@ -31,13 +31,22 @@ namespace Carpet {
   
   
   
+  typedef vect<int,dim> ivect;
+  typedef bbox<int,dim> ibbox;
+  
+  typedef vect<vect<bool,2>,dim> bvect;
+  
+  
+  
   static int (*regrid_routine) (const cGH * cckgGH,
 				gh<dim>::rexts& bbsss,
+				gh<dim>::rbnds& obss,
 				gh<dim>::rprocs& pss) = 0;
   
   
   
   static void CheckRegions (const gh<dim>::rexts& bbsss,
+			    const gh<dim>::rbnds& obss,
 			    const gh<dim>::rprocs& pss);
   static void Adapt (const cGH* cgh, int reflevels, gh<dim>* hh);
   
@@ -45,14 +54,17 @@ namespace Carpet {
   
   static void OutputGridStructure (const cGH *cgh,
 				   const gh<dim>::rexts& bbsss,
+				   const gh<dim>::rbnds& obss,
 				   const gh<dim>::rprocs& pss);
   
   
   
   static void SplitRegions_AlongZ       (const cGH* cgh,
-					 vector<bbox<int,dim> >& bbs);
+					 vector<ibbox>& bbs,
+					 vector<bvect>& obs);
   static void SplitRegions_AsSpecified  (const cGH* cgh,
-					 vector<bbox<int,dim> >& bbs);
+					 vector<ibbox>& bbs,
+					 vector<bvect>& obs);
   
   static void MakeProcessors_RoundRobin (const cGH* cgh,
 					 const gh<dim>::rexts& bbss,
@@ -60,7 +72,8 @@ namespace Carpet {
   
   
   
-  void CheckRegions (const gh<dim>::rexts& bbsss, const gh<dim>::rprocs& pss)
+  void CheckRegions (const gh<dim>::rexts& bbsss, const gh<dim>::rbnds& obss,
+		     const gh<dim>::rprocs& pss)
   {
     // At least one level
     assert (bbsss.size() > 0);
@@ -84,7 +97,9 @@ namespace Carpet {
     }
     
     assert (pss.size() == bbsss.size());
+    assert (obss.size() == bbsss.size());
     for (int rl=0; rl<(int)bbsss.size(); ++rl) {
+      assert (obss[rl].size() == bbsss[rl].size());
       assert (pss[rl].size() == bbsss[rl].size());
     }
   }
@@ -93,6 +108,7 @@ namespace Carpet {
   
   void RegisterRegridRoutine (int (*routine)(const cGH * cckgGH,
 					     gh<dim>::rexts& bbsss,
+					     gh<dim>::rbnds& obss,
 					     gh<dim>::rprocs& pss))
   {
     assert (!regrid_routine);
@@ -113,30 +129,32 @@ namespace Carpet {
     
     // Check whether to recompose
     gh<dim>::rexts bbsss;
+    gh<dim>::rbnds obss;
     gh<dim>::rprocs pss;
-    int do_recompose = (*regrid_routine) (cgh, bbsss, pss);
+    int do_recompose = (*regrid_routine) (cgh, bbsss, obss, pss);
     assert (do_recompose >= 0);
     if (do_recompose == 0) return;
-    Recompose (cgh, bbsss, pss);
+    Recompose (cgh, bbsss, obss, pss);
   }
   
   
   
   void Recompose (const cGH* const cgh,
 		  const gh<dim>::rexts& bbsss,
+		  const gh<dim>::rbnds& obss,
 		  const gh<dim>::rprocs& pss)
   {
     assert (mglevel == -1);
     assert (component == -1);
     
     // Check the regions
-    CheckRegions (bbsss, pss);
+    CheckRegions (bbsss, obss, pss);
     
     // Write grid structure to file
-    OutputGridStructure (cgh, bbsss, pss);
+    OutputGridStructure (cgh, bbsss, obss, pss);
     
     // Recompose
-    hh->recompose (bbsss, pss);
+    hh->recompose (bbsss, obss, pss);
     Output (cgh, hh, 0);
     
     // Adapt grid scalars
@@ -151,6 +169,7 @@ namespace Carpet {
       case CCTK_ARRAY:
 	Adapt (cgh, hh->reflevels(), arrdata[group].hh);
 	Output (cgh, arrdata[group].hh, CCTK_GroupName(group));
+	break;
       case CCTK_GF:
 	break;
       default:
@@ -167,21 +186,21 @@ namespace Carpet {
   static void Adapt (const cGH* cgh, const int reflevels, gh<dim>* hh)
   {
     const int nprocs   = CCTK_nProcs(cgh);
-    vector<vector<bbox<int,dim> > > bbss(reflevels);
+    vector<vector<ibbox> > bbss(reflevels);
     // note: what this routine calls "ub" is "ub+str" elsewhere
-    vect<int,dim> rstr = hh->baseextent.stride();
-    vect<int,dim> rlb  = hh->baseextent.lower();
-    vect<int,dim> rub  = hh->baseextent.upper() + rstr;
+    ivect rstr = hh->baseextent.stride();
+    ivect rlb  = hh->baseextent.lower();
+    ivect rub  = hh->baseextent.upper() + rstr;
     for (int rl=0; rl<reflevels; ++rl) {
       if (rl>0) {
 	// save old values
-	const vect<int,dim> oldrlb = rlb;
-	const vect<int,dim> oldrub = rub;
+	const ivect oldrlb = rlb;
+	const ivect oldrub = rub;
 	// calculate extent
-	const vect<int,dim> rextent = rub - rlb;
+	const ivect rextent = rub - rlb;
 	// calculate new extent
 	assert (all(rextent % hh->reffact == 0));
-	const vect<int,dim> newrextent = rextent / hh->reffact;
+	const ivect newrextent = rextent / hh->reffact;
 	// refined boxes have smaller stride
 	assert (all(rstr%hh->reffact == 0));
 	rstr /= hh->reffact;
@@ -191,11 +210,12 @@ namespace Carpet {
 	// require rub<oldrub because we really want rub-rstr<=oldrub-oldstr
 	assert (all(rlb >= oldrlb && rub < oldrub));
       }
-      vector<bbox<int,dim> > bbs(nprocs);
+      vector<ibbox> bbs(nprocs);
       for (int c=0; c<nprocs; ++c) {
-	vect<int,dim> cstr = rstr;
-	vect<int,dim> clb = rlb;
-	vect<int,dim> cub = rub;
+	ivect cstr = rstr;
+	ivect clb = rlb;
+	ivect cub = rub;
+	// split the components along the z axis
 	const int glonpz = (rub[dim-1] - rlb[dim-1]) / cstr[dim-1];
 	const int locnpz = (glonpz + nprocs - 1) / nprocs;
 	const int zstep = locnpz * cstr[dim-1];
@@ -205,7 +225,7 @@ namespace Carpet {
  	if (cub[dim-1] > rub[dim-1]) cub[dim-1] = rub[dim-1];
 	assert (clb[dim-1] <= cub[dim-1]);
 	assert (cub[dim-1] <= rub[dim-1]);
-	bbs[c] = bbox<int,dim>(clb, cub-cstr, cstr);
+	bbs[c] = ibbox(clb, cub-cstr, cstr);
       }
 //       bbss[rl] = bbs;
       bbss[rl].clear();
@@ -215,20 +235,28 @@ namespace Carpet {
       }
     }
     
-    vector<vector<vector<bbox<int,dim> > > > bbsss
+    vector<vector<vector<ibbox> > > bbsss
       = hh->make_multigrid_boxes(bbss, mglevels);
     
     vector<vector<int> > pss(bbss.size());
+    vector<vector<bvect> > obss(bbss.size());
     for (int rl=0; rl<reflevels; ++rl) {
       pss[rl] = vector<int>(bbss[rl].size());
+      obss[rl] = vector<bvect>(bbss[rl].size());
       // make sure all processors have the same number of components
       assert (bbss[rl].size() % nprocs == 0);
       for (int c=0; c<(int)bbss[rl].size(); ++c) {
-	pss[rl][c] = c % nprocs; // distribute among processors
+	// distribute among processors
+	pss[rl][c] = c % nprocs;
+	for (int d=0; d<dim; ++d) {
+	  // assume the components are split along the z axis
+	  obss[rl][c][d][0] = d<dim-1 || c==0;
+	  obss[rl][c][d][1] = d<dim-1 || c==(int)bbss[rl].size()-1;
+	}
       }
     }
     
-    hh->recompose(bbsss, pss);
+    hh->recompose(bbsss, obss, pss);
   }
   
   
@@ -284,6 +312,7 @@ namespace Carpet {
   
   static void OutputGridStructure (const cGH * const cgh,
 				   const gh<dim>::rexts& bbsss,
+				   const gh<dim>::rbnds& obss,
 				   const gh<dim>::rprocs& pss)
   {
     DECLARE_CCTK_PARAMETERS;
@@ -309,7 +338,7 @@ namespace Carpet {
 	file.open (filename, ios::out | ios::trunc);
 	assert (file.good());
 	file << "# grid structure" << endl
-	     << "# format: reflevel component mglevel   processor bounding-box" << endl;
+	     << "# format: reflevel component mglevel   processor bounding-box is-outer-boundary" << endl;
 	assert (file.good());
       }
     }
@@ -325,7 +354,7 @@ namespace Carpet {
       for (int c=0; c<(int)bbsss[rl].size(); ++c) {
 	file << rl << " " << c << " mglevels " << bbsss[rl][c].size() << endl;
 	for (int ml=0; ml<(int)bbsss[rl][c].size(); ++ml) {
-	  file << rl << " " << c << " " << ml << "   " << pss[rl][c] << " " << bbsss[rl][c][ml] << endl;
+	  file << rl << " " << c << " " << ml << "   " << pss[rl][c] << " " << bbsss[rl][c][ml] << obss[rl][c] << endl;
 	}
       }
     }
@@ -337,14 +366,14 @@ namespace Carpet {
   
   
   
-  void SplitRegions (const cGH* cgh, vector<bbox<int,dim> >& bbs)
+  void SplitRegions (const cGH* cgh, vector<ibbox>& bbs, vector<bvect>& obs)
   {
     DECLARE_CCTK_PARAMETERS;
     
     if (CCTK_EQUALS (processor_topology, "automatic")) {
-      SplitRegions_AlongZ (cgh, bbs);
+      SplitRegions_AlongZ (cgh, bbs, obs);
     } else if (CCTK_EQUALS (processor_topology, "manual")) {
-      SplitRegions_AsSpecified (cgh, bbs);
+      SplitRegions_AsSpecified (cgh, bbs, obs);
     } else {
       abort();
     }
@@ -352,7 +381,8 @@ namespace Carpet {
   
   
   
-  void SplitRegions_AlongZ (const cGH* cgh, vector<bbox<int,dim> >& bbs)
+  void SplitRegions_AlongZ (const cGH* cgh, vector<ibbox>& bbs,
+			    vector<bvect>& obs)
   {
     // Something to do?
     if (bbs.size() == 0) return;
@@ -363,15 +393,17 @@ namespace Carpet {
     
     assert (bbs.size() == 1);
     
-    const vect<int,dim> rstr = bbs[0].stride();
-    const vect<int,dim> rlb  = bbs[0].lower();
-    const vect<int,dim> rub  = bbs[0].upper() + rstr;
+    const ivect rstr = bbs[0].stride();
+    const ivect rlb  = bbs[0].lower();
+    const ivect rub  = bbs[0].upper() + rstr;
+    const bvect obnd = obs[0];
     
     bbs.resize(nprocs);
+    obs.resize(nprocs);
     for (int c=0; c<nprocs; ++c) {
-      vect<int,dim> cstr = rstr;
-      vect<int,dim> clb = rlb;
-      vect<int,dim> cub = rub;
+      ivect cstr = rstr;
+      ivect clb = rlb;
+      ivect cub = rub;
       const int glonpz = (rub[dim-1] - rlb[dim-1]) / cstr[dim-1];
       const int locnpz = (glonpz + nprocs - 1) / nprocs;
       const int zstep = locnpz * cstr[dim-1];
@@ -381,13 +413,17 @@ namespace Carpet {
       if (cub[dim-1] > rub[dim-1]) cub[dim-1] = rub[dim-1];
       assert (clb[dim-1] <= cub[dim-1]);
       assert (cub[dim-1] <= rub[dim-1]);
-      bbs[c] = bbox<int,dim>(clb, cub-cstr, cstr);
+      bbs[c] = ibbox(clb, cub-cstr, cstr);
+      obs[c] = obnd;
+      if (c>0)        obs[c][dim-1][0] = false;
+      if (c<nprocs-1) obs[c][dim-1][1] = false;
     }
   }
   
   
   
-  void SplitRegions_AsSpecified (const cGH* cgh, vector<bbox<int,dim> >& bbs)
+  void SplitRegions_AsSpecified (const cGH* cgh, vector<ibbox>& bbs,
+				 vector<bvect>& obs)
   {
     DECLARE_CCTK_PARAMETERS;
     
@@ -400,11 +436,12 @@ namespace Carpet {
     
     assert (bbs.size() == 1);
     
-    const vect<int,dim> rstr = bbs[0].stride();
-    const vect<int,dim> rlb  = bbs[0].lower();
-    const vect<int,dim> rub  = bbs[0].upper() + rstr;
+    const ivect rstr = bbs[0].stride();
+    const ivect rlb  = bbs[0].lower();
+    const ivect rub  = bbs[0].upper() + rstr;
+    const bvect obnd = obs[0];
     
-    const vect<int,dim> nprocs_dir
+    const ivect nprocs_dir
       (processor_topology_3d_x, processor_topology_3d_y,
        processor_topology_3d_z);
     assert (all (nprocs_dir > 0));
@@ -415,25 +452,33 @@ namespace Carpet {
     assert (prod(nprocs_dir) == nprocs);
     
     bbs.resize(nprocs);
+    obs.resize(nprocs);
     assert (dim==3);
     for (int k=0; k<nprocs_dir[2]; ++k) {
       for (int j=0; j<nprocs_dir[1]; ++j) {
 	for (int i=0; i<nprocs_dir[0]; ++i) {
 	  const int c = i + nprocs_dir[0] * (j + nprocs_dir[1] * k);
-	  const vect<int,dim> ipos (i, j, k);
-	  vect<int,dim> cstr = rstr;
-	  vect<int,dim> clb = rlb;
-	  vect<int,dim> cub = rub;
-	  const vect<int,dim> glonp = (rub - rlb) / cstr;
-	  const vect<int,dim> locnp = (glonp + nprocs_dir - 1) / nprocs_dir;
-	  const vect<int,dim> step = locnp * cstr;
+	  const ivect ipos (i, j, k);
+	  ivect cstr = rstr;
+	  ivect clb = rlb;
+	  ivect cub = rub;
+	  const ivect glonp = (rub - rlb) / cstr;
+	  const ivect locnp = (glonp + nprocs_dir - 1) / nprocs_dir;
+	  const ivect step = locnp * cstr;
 	  clb = rlb + step *  ipos;
 	  cub = rlb + step * (ipos+1);
 	  clb = min (clb, rub);
 	  cub = min (cub, rub);
 	  assert (all (clb <= cub));
 	  assert (all (cub <= rub));
-	  bbs[c] = bbox<int,dim>(clb, cub-cstr, cstr);
+	  bbs[c] = ibbox(clb, cub-cstr, cstr);
+	  obs[c] = obnd;
+	  if (i>0) obs[c][0][0] = false;
+	  if (j>0) obs[c][1][0] = false;
+	  if (k>0) obs[c][2][0] = false;
+	  if (i<nprocs_dir[0]-1) obs[c][0][1] = false;
+	  if (j<nprocs_dir[1]-1) obs[c][1][1] = false;
+	  if (k<nprocs_dir[2]-1) obs[c][2][1] = false;
 	}
       }
     }
diff --git a/Carpet/Carpet/src/SetupGH.cc b/Carpet/Carpet/src/SetupGH.cc
index 1bf19f64f..2fa0faf72 100644
--- a/Carpet/Carpet/src/SetupGH.cc
+++ b/Carpet/Carpet/src/SetupGH.cc
@@ -1,6 +1,9 @@
 #include <assert.h>
-#include <iostream.h>
 #include <stdlib.h>
+#include <string.h>
+
+#include <iostream>
+#include <sstream>
 
 #include "cctk.h"
 #include "cctk_Parameters.h"
@@ -12,7 +15,7 @@
 
 #include "carpet.hh"
 
-static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/SetupGH.cc,v 1.23 2002/01/14 16:40:45 schnetter Exp $";
+static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/SetupGH.cc,v 1.24 2002/03/11 13:17:10 schnetter Exp $";
 
 
 
@@ -246,14 +249,34 @@ namespace Carpet {
     mglevel   = -1;
     component = -1;
     
-    // Invent a refinement structure
-    vector<bbox<int,dim> > bbs(1);
-    bbs[0] = hh->baseextent;
+    // Set initial refinement structure
+    vector<bbox<int,dim> > bbs;
+    vector<vect<vect<bool,2>,dim> > obs;
+    if (strcmp(base_extents, "")==0) {
+      // default: one grid component covering everything
+      bbs.push_back (hh->baseextent);
+      obs.push_back (vect<vect<bool,2>,dim>(vect<bool,2>(true)));
+    } else {
+      // explicit grid components
+      istringstream istr (base_extents);
+      istr >> bbs;
+      CCTK_VInfo (CCTK_THORNSTRING, "Using %d grid patches", bbs.size());
+      cout << "grid-patches-are " << bbs << endl;
+      if (bbs.size()<=0) {
+	CCTK_WARN (0, "Cannot evolve with 0 grid patches");
+      }
+      istringstream istr2 (base_bboxes);
+      istr2 >> obs;
+      cout << "outer-boundaries-are " << obs << endl;
+      assert (obs.size() == bbs.size());
+    }
     
-    SplitRegions (cgh, bbs);
+    SplitRegions (cgh, bbs, obs);
     
     vector<vector<bbox<int,dim> > > bbss(1);
+    vector<vector<vect<vect<bool,2>,dim> > > obss(1);
     bbss[0] = bbs;
+    obss[0] = obs;
     
     gh<dim>::rexts bbsss;
     bbsss = hh->make_multigrid_boxes(bbss, mglevels);
@@ -262,7 +285,7 @@ namespace Carpet {
     MakeProcessors (cgh, bbsss, pss);
     
     // Recompose grid hierarchy
-    Recompose (cgh, bbsss, pss);
+    Recompose (cgh, bbsss, obss, pss);
     
     // Initialise time step on coarse grid
     base_delta_time = 0;
diff --git a/Carpet/Carpet/src/carpet.hh b/Carpet/Carpet/src/carpet.hh
index f8803ae10..6fded25f4 100644
--- a/Carpet/Carpet/src/carpet.hh
+++ b/Carpet/Carpet/src/carpet.hh
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/carpet.hh,v 1.15 2002/01/11 17:19:46 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/carpet.hh,v 1.16 2002/03/11 13:17:10 schnetter Exp $
 
 #include "Carpet/CarpetLib/src/gh.hh"
 
@@ -12,6 +12,7 @@ namespace Carpet {
   
   void Recompose (const cGH* cgh,
 		  const gh<dim>::rexts& bbsss,
+		  const gh<dim>::rbnds& obss,
 		  const gh<dim>::rprocs& pss);
   
   enum checktimes { currenttime,
diff --git a/Carpet/Carpet/src/carpet_public.hh b/Carpet/Carpet/src/carpet_public.hh
index b2bf7e8bb..b628dd78c 100644
--- a/Carpet/Carpet/src/carpet_public.hh
+++ b/Carpet/Carpet/src/carpet_public.hh
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/carpet_public.hh,v 1.16 2002/01/11 17:19:46 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/carpet_public.hh,v 1.17 2002/03/11 13:17:11 schnetter Exp $
 
 // It is assumed that the number of components of all arrays is equal
 // to the number of components of the grid functions, and that their
@@ -139,9 +139,11 @@ namespace Carpet {
   // Functions for recomposing the grid hierarchy
   void RegisterRegridRoutine (int (*routine)(const cGH * cckgGH,
 					     gh<dim>::rexts& bbsss,
+					     gh<dim>::rbnds& obss,
 					     gh<dim>::rprocs& pss));
   
-  void SplitRegions (const cGH* cgh, vector<bbox<int,dim> >& bbs);
+  void SplitRegions (const cGH* cgh, vector<bbox<int,dim> >& bbs,
+		     vector<vect<vect<bool,2>,dim> >& obs);
   
   void MakeProcessors (const cGH* cgh, const gh<dim>::rexts& bbsss,
 		       gh<dim>::rprocs& pss);
diff --git a/Carpet/Carpet/src/helpers.cc b/Carpet/Carpet/src/helpers.cc
index b412214fc..a2046a198 100644
--- a/Carpet/Carpet/src/helpers.cc
+++ b/Carpet/Carpet/src/helpers.cc
@@ -12,7 +12,7 @@
 
 #include "carpet.hh"
 
-static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/helpers.cc,v 1.18 2002/01/14 15:56:02 schnetter Exp $";
+static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/helpers.cc,v 1.19 2002/03/11 13:17:11 schnetter Exp $";
 
 
 
@@ -320,8 +320,6 @@ namespace Carpet {
 	assert (reflevel < (int)dd->boxes.size());
 	assert (component < (int)dd->boxes[reflevel].size());
 	assert (mglevel < (int)dd->boxes[reflevel][component].size());
-	const bbox<int,dim>& bext = hh->baseextent;
-	const bbox<int,dim>& iext = hh->extents[reflevel][component][mglevel];
 	const bbox<int,dim>& ext
 	  = dd->boxes[reflevel][component][mglevel].exterior;
 	for (int d=0; d<dim; ++d) {
@@ -332,16 +330,8 @@ namespace Carpet {
 // 	  assert (cgh->cctk_lbnd[d]>=0 && cgh->cctk_ubnd[d]<cgh->cctk_gsh[d]);
 	  assert (cgh->cctk_ubnd[d]-cgh->cctk_lbnd[d]+1 == cgh->cctk_lsh[d]);
 	  assert (cgh->cctk_lbnd[d]<=cgh->cctk_ubnd[d]+1);
-	  // Do not allow outer boundaries on the finer grids
-	  if (reflevel==0) {
-	    assert (iext.lower()[d] >= bext.lower()[d]);
-	    assert (iext.upper()[d] <= bext.upper()[d]);
-	    cgh->cctk_bbox[2*d  ] = iext.lower()[d] == bext.lower()[d];
-	    cgh->cctk_bbox[2*d+1] = iext.upper()[d] == bext.upper()[d];
-	  } else {
-	    cgh->cctk_bbox[2*d  ] = 0;
-	    cgh->cctk_bbox[2*d+1] = 0;
-	  }
+	  cgh->cctk_bbox[2*d  ] = hh->outer_boundaries[reflevel][component][d][0];
+	  cgh->cctk_bbox[2*d+1] = hh->outer_boundaries[reflevel][component][d][1];
 	  for (int stg=0; stg<CCTK_NSTAGGER; ++stg) {
 	    // TODO: support staggering
 	    cgh->cctk_lssh[CCTK_LSSH_IDX(stg,d)] = cgh->cctk_lsh[d];
diff --git a/Carpet/CarpetLib/src/bbox.cc b/Carpet/CarpetLib/src/bbox.cc
index 337857826..27c1e784f 100644
--- a/Carpet/CarpetLib/src/bbox.cc
+++ b/Carpet/CarpetLib/src/bbox.cc
@@ -5,7 +5,7 @@
     copyright            : (C) 2000 by Erik Schnetter
     email                : schnetter@astro.psu.edu
 
-    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/bbox.cc,v 1.7 2001/07/02 13:22:11 schnetter Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/bbox.cc,v 1.8 2002/03/11 13:17:12 schnetter Exp $
 
  ***************************************************************************/
 
@@ -225,6 +225,28 @@ bbox<T,D>::iterator bbox<T,D>::end () const {
 
 
 
+// Input
+template<class T,int D>
+void bbox<T,D>::input (istream& is) {
+  skipws (is);
+  assert (is.peek() == '(');
+  is.get();
+  is >> _lower;
+  skipws (is);
+  assert (is.peek() == ':');
+  is.get();
+  is >> _upper;
+  skipws (is);
+  assert (is.peek() == ':');
+  is.get();
+  is >> _stride;
+  skipws (is);
+  assert (is.peek() == ')');
+  is.get();
+}
+
+
+
 // Output
 template<class T,int D>
 void bbox<T,D>::output (ostream& os) const {
diff --git a/Carpet/CarpetLib/src/bbox.hh b/Carpet/CarpetLib/src/bbox.hh
index 160748f27..04db295db 100644
--- a/Carpet/CarpetLib/src/bbox.hh
+++ b/Carpet/CarpetLib/src/bbox.hh
@@ -5,7 +5,7 @@
     copyright            : (C) 2000 by Erik Schnetter
     email                : schnetter@astro.psu.edu
 
-    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/bbox.hh,v 1.8 2001/03/27 22:26:31 eschnett Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/bbox.hh,v 1.9 2002/03/11 13:17:12 schnetter Exp $
 
  ***************************************************************************/
 
@@ -33,7 +33,9 @@ using namespace std;
 // Forward declaration
 template<class T, int D> class bbox;
 
-// Output
+// Input/Output
+template<class T, int D>
+istream& operator>> (istream& is, bbox<T,D>& b);
 template<class T, int D>
 ostream& operator<< (ostream& os, const bbox<T,D>& b);
 
@@ -116,12 +118,22 @@ public:
   iterator begin () const;
   iterator end () const;
   
-  // Output
+  // Input/Output
+  void input (istream& is);
   void output (ostream& os) const;
 };
 
 
 
+// Input
+template<class T,int D>
+inline istream& operator>> (istream& is, bbox<T,D>& b) {
+  b.input(is);
+  return is;
+}
+
+
+
 // Output
 template<class T,int D>
 inline ostream& operator<< (ostream& os, const bbox<T,D>& b) {
diff --git a/Carpet/CarpetLib/src/defs.cc b/Carpet/CarpetLib/src/defs.cc
index f862195d3..5a913beb0 100644
--- a/Carpet/CarpetLib/src/defs.cc
+++ b/Carpet/CarpetLib/src/defs.cc
@@ -5,7 +5,7 @@
     copyright            : (C) 2000 by Erik Schnetter
     email                : schnetter@astro.psu.edu
 
-    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/defs.cc,v 1.8 2001/08/26 13:59:31 schnetter Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/defs.cc,v 1.9 2002/03/11 13:17:12 schnetter Exp $
 
  ***************************************************************************/
 
@@ -19,6 +19,7 @@
  ***************************************************************************/
 
 #include <assert.h>
+#include <ctype.h>
 
 #include <iostream>
 #include <list>
@@ -33,6 +34,38 @@ using namespace std;
 
 
 
+void skipws (istream& is) {
+  while (is.good() && isspace(is.peek())) {
+    is.get();
+  }
+}
+
+
+
+// Vector input
+template<class T>
+istream& input (istream& is, vector<T>& v) {
+  v.clear();
+  skipws (is);
+  assert (is.peek() == '[');
+  is.get();
+  skipws (is);
+  while (is.good() && is.peek() != ']') {
+    v.push_back ();
+    is >> v[v.size()-1];
+    skipws (is);
+    if (is.peek() != ',') break;
+    is.get();
+    skipws (is);
+  }
+  skipws (is);
+  assert (is.peek() == ']');
+  is.get();
+  return is;
+}
+
+
+
 // List output
 template<class T>
 ostream& output (ostream& os, const list<T>& l) {
@@ -76,6 +109,10 @@ ostream& output (ostream& os, const vector<T>& v) {
 #include "bbox.hh"
 #include "bboxset.hh"
 
+template istream& input (istream& os, vector<bbox<int,3> >& v);
+template istream& input (istream& os, vector<vector<bbox<int,3> > >& v);
+template istream& input (istream& os, vector<vector<vect<vect<bool,2>,3> > >& v);
+
 template ostream& output (ostream& os, const list<bbox<int,3> >& l);
 template ostream& output (ostream& os, const set<bbox<int,3> >& s);
 template ostream& output (ostream& os, const set<bboxset<int,3> >& s);
@@ -84,5 +121,6 @@ template ostream& output (ostream& os, const vector<bbox<int,3> >& v);
 template ostream& output (ostream& os, const vector<list<bbox<int,3> > >& v);
 template ostream& output (ostream& os, const vector<vector<int> >& v);
 template ostream& output (ostream& os, const vector<vector<bbox<int,3> > >& v);
+template ostream& output (ostream& os, const vector<vector<vect<vect<bool,2>,3> > >& v);
 template ostream& output (ostream& os, const vector<vector<vector<bbox<int,3> > > >& v);
 #endif
diff --git a/Carpet/CarpetLib/src/defs.hh b/Carpet/CarpetLib/src/defs.hh
index a85b164c1..8f1810bcb 100644
--- a/Carpet/CarpetLib/src/defs.hh
+++ b/Carpet/CarpetLib/src/defs.hh
@@ -5,7 +5,7 @@
     copyright            : (C) 2000 by Erik Schnetter
     email                : schnetter@astro.psu.edu
 
-    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/defs.hh,v 1.6 2002/01/09 13:56:26 schnetter Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/defs.hh,v 1.7 2002/03/11 13:17:13 schnetter Exp $
 
  ***************************************************************************/
 
@@ -62,6 +62,23 @@ inline T ipow (const T& x, const int y) {
   }
 }
 
+
+
+// Skip whitespace
+void skipws (istream& is);
+
+
+
+// Container input
+template<class T> istream& input (istream& is, vector<T>& v);
+
+template<class T>
+inline istream& operator>> (istream& is, vector<T>& v) {
+  return input(is,v);
+}
+
+
+
 // Container output
 template<class T> ostream& output (ostream& os, const list<T>& l);
 template<class T> ostream& output (ostream& os, const set<T>& s);
diff --git a/Carpet/CarpetLib/src/gh.cc b/Carpet/CarpetLib/src/gh.cc
index 95bf16350..c0cd6df08 100644
--- a/Carpet/CarpetLib/src/gh.cc
+++ b/Carpet/CarpetLib/src/gh.cc
@@ -7,7 +7,7 @@
     copyright            : (C) 2000 by Erik Schnetter
     email                : schnetter@astro.psu.edu
 
-    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gh.cc,v 1.11 2002/01/09 13:56:27 schnetter Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gh.cc,v 1.12 2002/03/11 13:17:13 schnetter Exp $
 
  ***************************************************************************/
 
@@ -52,8 +52,10 @@ gh<D>::~gh () { }
 
 // Modifiers
 template<int D>
-void gh<D>::recompose (const rexts& exts, const rprocs& procs) {
+void gh<D>::recompose (const rexts& exts, const rbnds& outer_bounds,
+		       const rprocs& procs) {
   extents = exts;
+  outer_boundaries = outer_bounds;
   processors = procs;
   
   // Consistency checks
@@ -63,8 +65,10 @@ void gh<D>::recompose (const rexts& exts, const rprocs& procs) {
   // Check processor number consistency
   for (int rl=0; rl<reflevels(); ++rl) {
     assert (processors.size() == extents.size());
+    assert (outer_boundaries.size() == extents.size());
     for (int c=0; c<components(rl); ++c) {
-      assert (procs[rl].size() == extents[rl].size());
+      assert (processors[rl].size() == extents[rl].size());
+      assert (outer_boundaries[rl].size() == extents[rl].size());
     }
   }
   
diff --git a/Carpet/CarpetLib/src/gh.hh b/Carpet/CarpetLib/src/gh.hh
index 4748fb1de..a1f917f95 100644
--- a/Carpet/CarpetLib/src/gh.hh
+++ b/Carpet/CarpetLib/src/gh.hh
@@ -7,7 +7,7 @@
     copyright            : (C) 2000 by Erik Schnetter
     email                : schnetter@astro.psu.edu
 
-    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gh.hh,v 1.8 2001/12/14 16:39:43 schnetter Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gh.hh,v 1.9 2002/03/11 13:17:13 schnetter Exp $
 
  ***************************************************************************/
 
@@ -56,19 +56,26 @@ public:
   typedef vect<int,D> ivect;
   typedef bbox<int,D> ibbox;
   
-  typedef vector<ibbox> mexts;          // ... for each multigrid level
-  typedef vector<mexts> cexts;          // ... for each component
-  typedef vector<cexts> rexts;          // ... for each refinement level
+  typedef vect<vect<bool,2>,D> bvect;
   
-  typedef vector<int>    cprocs;        // ... for each component
-  typedef vector<cprocs> rprocs;        // ... for each refinement level
+  typedef vector<ibbox> mexts;	// ... for each multigrid level
+  typedef vector<mexts> cexts;	// ... for each component
+  typedef vector<cexts> rexts;	// ... for each refinement level
+  
+  typedef vector<bvect> cbnds;	// ... for each component
+  typedef vector<cbnds> rbnds;	// ... for each refinement level
+  
+  typedef vector<int>    cprocs; // ... for each component
+  typedef vector<cprocs> rprocs; // ... for each refinement level
   
 public:				// should be readonly
   
   ibbox baseextent;		// bounds (inclusive) of base level
   vector<vector<ibbox> > bases; // [rl][ml]
   
-  rexts extents;		// bounds of all grids
+  // TODO: invent structure for this
+  rexts extents;		// extents of all grids
+  rbnds outer_boundaries;	// boundary descriptions of all grids
   rprocs processors;		// processor numbers of all grids
   
   list<dh<D>*> dhs;		// list of all data hierarchies
@@ -84,7 +91,8 @@ public:
   virtual ~gh ();
   
   // Modifiers
-  void recompose (const rexts& exts, const rprocs& procs);
+  void recompose (const rexts& exts, const rbnds& outer_bounds,
+		  const rprocs& procs);
   
   // Helpers
   cexts make_reflevel_multigrid_boxes (const vector<ibbox>& exts,
@@ -110,6 +118,12 @@ public:
     return (int)extents[rl][c].size();
   }
   
+  bvect outer_boundary (const int rl, const int c) const {
+    assert (rl>=0 && rl<reflevels());
+    assert (c>=0 && c<components(rl));
+    return outer_boundaries[rl][c];
+  }
+  
   int proc (const int rl, const int c) const {
     assert (rl>=0 && rl<reflevels());
     assert (c>=0 && c<components(rl));
diff --git a/Carpet/CarpetLib/src/vect.cc b/Carpet/CarpetLib/src/vect.cc
index 132e37afa..85a608ee9 100644
--- a/Carpet/CarpetLib/src/vect.cc
+++ b/Carpet/CarpetLib/src/vect.cc
@@ -5,7 +5,7 @@
     copyright            : (C) 2000 by Erik Schnetter
     email                : schnetter@astro.psu.edu
 
-    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/vect.cc,v 1.6 2002/01/08 12:03:55 schnetter Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/vect.cc,v 1.7 2002/03/11 13:17:13 schnetter Exp $
 
  ***************************************************************************/
 
@@ -32,13 +32,34 @@ using namespace std;
 
 
 
+// Input
+template<class T,int D>
+void vect<T,D>::input (istream& is) {
+  skipws (is);
+  assert (is.peek() == '[');
+  is.get();
+  for (int d=0; d<D; ++d) {
+    is >> (*this)[d];
+    if (d<D-1) {
+      skipws (is);
+      assert (is.peek() == ',');
+      is.get();
+    }
+  }
+  skipws (is);
+  assert (is.peek() == ']');
+  is.get();
+}
+
+
+
 // Output
 template<class T,int D>
 void vect<T,D>::output (ostream& os) const {
   os << "[";
   for (int d=0; d<D; ++d) {
-    if (d>0) os << ",";
     os << (*this)[d];
+    if (d<D-1) os << ",";
   }
   os << "]";
 }
@@ -52,6 +73,9 @@ template class vect<int,1>;
 template class vect<int,2>;
 template class vect<int,3>;
 
+template void vect<double,3>::input (istream& is);
+template void vect<vect<bool,2>,3>::input (istream& is);
 template void vect<double,3>::output (ostream& os) const;
+template void vect<vect<bool,2>,3>::output (ostream& os) const;
 
 #endif
diff --git a/Carpet/CarpetLib/src/vect.hh b/Carpet/CarpetLib/src/vect.hh
index f53038741..beee9608d 100644
--- a/Carpet/CarpetLib/src/vect.hh
+++ b/Carpet/CarpetLib/src/vect.hh
@@ -5,7 +5,7 @@
     copyright            : (C) 2000 by Erik Schnetter
     email                : schnetter@astro.psu.edu
 
-    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/vect.hh,v 1.8 2002/01/09 17:45:42 schnetter Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/vect.hh,v 1.9 2002/03/11 13:17:13 schnetter Exp $
 
  ***************************************************************************/
 
@@ -33,7 +33,9 @@ using namespace std;
 // Forward definition
 template<class T, int D> class vect;
 
-// Output
+// Input/Output
+template<class T,int D>
+istream& operator>> (istream& is, vect<T,D>& a);
 template<class T,int D>
 ostream& operator<< (ostream& os, const vect<T,D>& a);
 
@@ -444,6 +446,8 @@ public:
     return r;
   }
   
+  // Input/Output
+  void input (istream& is);
   void output (ostream& os) const;
 };
 
@@ -598,6 +602,16 @@ inline vect<TT,D> scan1 (TT (* const func)(TT val, T x), TT val,
 
 
 
+// Input
+template<class T,int D>
+inline istream& operator>> (istream& is, vect<T,D>& a) {
+  a.input(is);
+  return is;
+}
+
+
+
+// Output
 template<class T,int D>
 inline ostream& operator<< (ostream& os, const vect<T,D>& a) {
   a.output(os);
diff --git a/Carpet/CarpetReduce/src/reduce.cc b/Carpet/CarpetReduce/src/reduce.cc
index 4cff31380..16b48fd88 100644
--- a/Carpet/CarpetReduce/src/reduce.cc
+++ b/Carpet/CarpetReduce/src/reduce.cc
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetReduce/src/reduce.cc,v 1.4 2002/03/06 17:46:15 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetReduce/src/reduce.cc,v 1.5 2002/03/11 13:17:14 schnetter Exp $
 
 #include <assert.h>
 #include <limits.h>
@@ -17,7 +17,7 @@
 
 #include "reduce.hh"
 
-static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetReduce/src/reduce.cc,v 1.4 2002/03/06 17:46:15 schnetter Exp $";
+static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetReduce/src/reduce.cc,v 1.5 2002/03/11 13:17:14 schnetter Exp $";
 
 
 
@@ -41,7 +41,7 @@ namespace CarpetReduce {
   
   // Poor man's RTTI
   enum ared { do_count, do_minimum, do_maximum, do_sum, do_sum_abs,
-	      do_average, do_norm1, do_norm2, do_norm_inf };
+	      do_sum_squared, do_average, do_norm1, do_norm2, do_norm_inf };
   
   
   
@@ -174,6 +174,18 @@ namespace CarpetReduce {
     MPI_Op mpi_op () const { return MPI_SUM; }
   };
   
+  struct sum_squared : reduction {
+    ared thered () const { return do_sum_squared; }
+    bool uses_cnt () const { return false; }
+    template<class T>
+    struct op {
+      static inline void initialise (T& accum) { accum = 0; }
+      static inline void reduce (T& accum, const T& val) { accum += val*val; }
+      static inline void finalise (T& accum, const T& cnt) { }
+    };
+    MPI_Op mpi_op () const { return MPI_SUM; }
+  };
+  
   struct average : reduction {
     ared thered () const { return do_average; }
     bool uses_cnt () const { return true; }
@@ -302,6 +314,7 @@ namespace CarpetReduce {
 	  INITIALISE(maximum,T);		\
 	  INITIALISE(sum,T);			\
 	  INITIALISE(sum_abs,T);		\
+	  INITIALISE(sum_squared,T);		\
 	  INITIALISE(average,T);		\
 	  INITIALISE(norm1,T);			\
 	  INITIALISE(norm2,T);			\
@@ -373,6 +386,7 @@ namespace CarpetReduce {
 	  REDUCE(maximum,T);			\
 	  REDUCE(sum,T);			\
 	  REDUCE(sum_abs,T);			\
+	  REDUCE(sum_squared,T);		\
 	  REDUCE(average,T);			\
 	  REDUCE(norm1,T);			\
 	  REDUCE(norm2,T);			\
@@ -460,6 +474,7 @@ namespace CarpetReduce {
 	    FINALISE(maximum,T);		\
 	    FINALISE(sum,T);			\
 	    FINALISE(sum_abs,T);		\
+	    FINALISE(sum_squared,T);		\
 	    FINALISE(average,T);		\
 	    FINALISE(norm1,T);			\
 	    FINALISE(norm2,T);			\
@@ -671,6 +686,7 @@ namespace CarpetReduce {
   REDUCTION(maximum);
   REDUCTION(sum);
   REDUCTION(sum_abs);
+  REDUCTION(sum_squared);
   REDUCTION(average);
   REDUCTION(norm1);
   REDUCTION(norm2);
@@ -682,25 +698,27 @@ namespace CarpetReduce {
   
   int CarpetReduceStartup ()
   {
-    CCTK_RegisterReductionOperator (count_GVs,    "count");
-    CCTK_RegisterReductionOperator (minimum_GVs,  "minimum");
-    CCTK_RegisterReductionOperator (maximum_GVs,  "maximum");
-    CCTK_RegisterReductionOperator (sum_GVs,      "sum");
-    CCTK_RegisterReductionOperator (sum_abs_GVs,  "sum_abs");
-    CCTK_RegisterReductionOperator (average_GVs,  "average");
-    CCTK_RegisterReductionOperator (norm1_GVs,    "norm1");
-    CCTK_RegisterReductionOperator (norm2_GVs,    "norm2");
-    CCTK_RegisterReductionOperator (norm_inf_GVs, "norm_inf");
-    
-    CCTK_RegisterReductionArrayOperator (count_arrays,    "count");
-    CCTK_RegisterReductionArrayOperator (minimum_arrays,  "minimum");
-    CCTK_RegisterReductionArrayOperator (maximum_arrays,  "maximum");
-    CCTK_RegisterReductionArrayOperator (sum_arrays,      "sum");
-    CCTK_RegisterReductionArrayOperator (sum_abs_arrays,  "sum_abs");
-    CCTK_RegisterReductionArrayOperator (average_arrays,  "average");
-    CCTK_RegisterReductionArrayOperator (norm1_arrays,    "norm1");
-    CCTK_RegisterReductionArrayOperator (norm2_arrays,    "norm2");
-    CCTK_RegisterReductionArrayOperator (norm_inf_arrays, "norm_inf");
+    CCTK_RegisterReductionOperator (count_GVs,       "count");
+    CCTK_RegisterReductionOperator (minimum_GVs,     "minimum");
+    CCTK_RegisterReductionOperator (maximum_GVs,     "maximum");
+    CCTK_RegisterReductionOperator (sum_GVs,         "sum");
+    CCTK_RegisterReductionOperator (sum_abs_GVs,     "sum_abs");
+    CCTK_RegisterReductionOperator (sum_squared_GVs, "sum_squared");
+    CCTK_RegisterReductionOperator (average_GVs,     "average");
+    CCTK_RegisterReductionOperator (norm1_GVs,       "norm1");
+    CCTK_RegisterReductionOperator (norm2_GVs,       "norm2");
+    CCTK_RegisterReductionOperator (norm_inf_GVs,    "norm_inf");
+    
+    CCTK_RegisterReductionArrayOperator (count_arrays,       "count");
+    CCTK_RegisterReductionArrayOperator (minimum_arrays,     "minimum");
+    CCTK_RegisterReductionArrayOperator (maximum_arrays,     "maximum");
+    CCTK_RegisterReductionArrayOperator (sum_arrays,         "sum");
+    CCTK_RegisterReductionArrayOperator (sum_abs_arrays,     "sum_abs");
+    CCTK_RegisterReductionArrayOperator (sum_squared_arrays, "sum_squared");
+    CCTK_RegisterReductionArrayOperator (average_arrays,     "average");
+    CCTK_RegisterReductionArrayOperator (norm1_arrays,       "norm1");
+    CCTK_RegisterReductionArrayOperator (norm2_arrays,       "norm2");
+    CCTK_RegisterReductionArrayOperator (norm_inf_arrays,    "norm_inf");
     
     return 0;
   }
diff --git a/Carpet/CarpetRegrid/param.ccl b/Carpet/CarpetRegrid/param.ccl
index e6902bc80..951e7949b 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.3 2002/01/11 17:37:13 schnetter Exp $
+# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/param.ccl,v 1.4 2002/03/11 13:17:15 schnetter Exp $
 
 
 
@@ -25,6 +25,8 @@ KEYWORD refined_regions "Regions where the grid is refined"
   "centre"             :: "Refine around the centre of the grid only"
   "manual-gridpoints"  :: "Refine the regions specified by integer grid points l[123]i[xyz]{min,max}"
   "manual-coordinates" :: "Refine the regions specified by coordinates l[123][xyz]{min,max}"
+  "manual-gridpoint-list"  :: ""
+  "manual-coordinate-list" :: ""
   "automatic"          :: "Refine automatically"
 } "centre"
 
@@ -202,6 +204,24 @@ CCTK_REAL l3zmax "Upper boundary of level 3 box in z-direction"
 
 
 
+# Refinement criteria for manual-gridpoint-list
+
+CCTK_STRING gridpoints "List of bounding box gridpoints"
+{
+  .* :: "[ [ ([<imin>,<jmin>,<kmin>]:[<imax>,<jmax>,<kmax>]:[<istride>,<jstride>,<kstride>]), ... ], ... ]"
+} ""
+
+
+
+# Refinement criteria for manual-coordinate-list
+
+CCTK_STRING coordinates "List of bounding box coordinates"
+{
+  .* :: "[ [ ([<xmin>,<ymin>,<zmin>]:[<xmax>,<ymax>,<zmax>]:[<xstride>,<ystride>,<zstride>]), ... ], ... ]"
+} ""
+
+
+
 # Refinement criteria for automatic refining
 
 CCTK_INT minwidth "Minimum width of refined region"
diff --git a/Carpet/CarpetRegrid/src/regrid.cc b/Carpet/CarpetRegrid/src/regrid.cc
index 1b0b95d87..ad17ca833 100644
--- a/Carpet/CarpetRegrid/src/regrid.cc
+++ b/Carpet/CarpetRegrid/src/regrid.cc
@@ -19,7 +19,7 @@
 #include "carpet.hh"
 #include "regrid.hh"
 
-static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/src/regrid.cc,v 1.8 2002/03/06 17:47:58 schnetter Exp $";
+static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/src/regrid.cc,v 1.9 2002/03/11 13:17:15 schnetter Exp $";
 
 
 
@@ -30,6 +30,13 @@ namespace CarpetRegrid {
   
   
   
+  typedef vect<int,dim> ivect;
+  typedef bbox<int,dim> ibbox;
+  
+  typedef vect<vect<bool,2>,dim> bvect;
+  
+  
+  
   int CarpetRegridStartup ()
   {
     RegisterRegridRoutine (CarpetRegridRegrid);
@@ -40,6 +47,7 @@ namespace CarpetRegrid {
   
   int CarpetRegridRegrid (const cGH * const cctkGH,
 			  gh<dim>::rexts& bbsss,
+			  gh<dim>::rbnds& obss,
 			  gh<dim>::rprocs& pss)
   {
     DECLARE_CCTK_PARAMETERS;
@@ -56,9 +64,6 @@ namespace CarpetRegrid {
 #endif
     assert (mglevel == -1);
     
-    // Return if this is the finest possible level
-    if (reflevel == maxreflevels-1) return 0;
-    
     // 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;
@@ -69,7 +74,7 @@ namespace CarpetRegrid {
       return 0;
     }
     
-    list<bbox<int,dim> > bbl;
+    list<ibbox> bbl;
     
     if (CCTK_EQUALS(refined_regions, "none")) {
       
@@ -85,13 +90,13 @@ namespace CarpetRegrid {
 	CCTK_WARN (0, "Cannot currently specify manual refinement regions for more than 4 refinement levels");
       }
       assert (refinement_levels<=4);
-      vector<vect<int,dim> > lower(3), upper(3);
-      lower[0] = vect<int,dim> (l1ixmin, l1iymin, l1izmin);
-      upper[0] = vect<int,dim> (l1ixmax, l1iymax, l1izmax);
-      lower[1] = vect<int,dim> (l2ixmin, l2iymin, l2izmin);
-      upper[1] = vect<int,dim> (l2ixmax, l2iymax, l2izmax);
-      lower[2] = vect<int,dim> (l3ixmin, l3iymin, l3izmin);
-      upper[2] = vect<int,dim> (l3ixmax, l3iymax, l3izmax);
+      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, refinement_levels, lower, upper, bbl);
       
     } else if (CCTK_EQUALS(refined_regions, "manual-coordinates")) {
@@ -109,6 +114,14 @@ namespace CarpetRegrid {
       upper[2] = vect<CCTK_REAL,dim> (l3xmax, l3ymax, l3zmax);
       MakeRegions_AsSpecified (cctkGH, refinement_levels, lower, upper, bbl);
       
+    } else if (CCTK_EQUALS(refined_regions, "manual-gridpoint-list")) {
+      
+      abort ();
+      
+    } else if (CCTK_EQUALS(refined_regions, "manual-coordinate-list")) {
+      
+      abort ();
+      
     } else if (CCTK_EQUALS(refined_regions, "automatic")) {
       
       const int vi = CCTK_VarIndex (errorvar);
@@ -135,9 +148,10 @@ namespace CarpetRegrid {
       abort();
     }
     
-    vector<bbox<int,dim> > bbs;
+    // transform bbox list into bbox vector
+    vector<ibbox> bbs;
     bbs.reserve (bbl.size());
-    for (list<bbox<int,dim> >::const_iterator it = bbl.begin();
+    for (list<ibbox>::const_iterator it = bbl.begin();
 	 it != bbl.end();
 	 ++it) {
       bbs.push_back (*it);
@@ -145,21 +159,40 @@ namespace CarpetRegrid {
     
     // TODO: ensure nesting properties
     
-    SplitRegions (cctkGH, bbs);
+    // TODO: set outer boundaries correctly
+    vector<bvect> obs;
+    obs.resize (bbs.size());
+    for (vector<bvect>::iterator it = obs.begin();
+	 it != obs.end();
+	 ++it) {
+      *it = bvect(vect<bool,2>(false));
+    }
+    
+    // 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;
       }
     }
     
@@ -170,7 +203,7 @@ namespace CarpetRegrid {
   
   
   
-  void MakeRegions_BaseLevel (const cGH* cctkGH, list<bbox<int,dim> >& bbl)
+  void MakeRegions_BaseLevel (const cGH* cctkGH, list<ibbox>& bbl)
   {
     assert (bbl.empty());
   }
@@ -181,27 +214,27 @@ namespace CarpetRegrid {
   // 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<bbox<int,dim> >& bbl)
+				 list<ibbox>& bbl)
   {
     assert (bbl.empty());
     
     if (reflevel+1 >= reflevels) return;
       
     // note: what this routine calls "ub" is "ub+str" elsewhere
-    vect<int,dim> rstr = hh->baseextent.stride();
-    vect<int,dim> rlb  = hh->baseextent.lower();
-    vect<int,dim> rub  = hh->baseextent.upper() + rstr;
+    ivect rstr = hh->baseextent.stride();
+    ivect rlb  = hh->baseextent.lower();
+    ivect rub  = hh->baseextent.upper() + rstr;
     
     for (int rl=0; rl<reflevel+1; ++rl) {
       // save old values
-      const vect<int,dim> oldrlb = rlb;
-      const vect<int,dim> oldrub = rub;
+      const ivect oldrlb = rlb;
+      const ivect oldrub = rub;
       // calculate extent and centre
-      const vect<int,dim> rextent = rub - rlb;
-      const vect<int,dim> rcentre = rlb + (rextent / 2 / rstr) * rstr;
+      const ivect rextent = rub - rlb;
+      const ivect rcentre = rlb + (rextent / 2 / rstr) * rstr;
       // calculate new extent
       assert (all(rextent % hh->reffact == 0));
-      const vect<int,dim> newrextent = rextent / hh->reffact;
+      const ivect newrextent = rextent / hh->reffact;
       // refined boxes have smaller stride
       assert (all(rstr%hh->reffact == 0));
       rstr /= hh->reffact;
@@ -219,32 +252,32 @@ namespace CarpetRegrid {
       assert (all(rlb >= oldrlb && rub < oldrub));
     }
     
-    bbl.push_back (bbox<int,dim>(rlb, rub-rstr, rstr));
+    bbl.push_back (ibbox(rlb, rub-rstr, rstr));
   }
   
   
   
   static void
   MakeRegions_AsSpecified_OneLevel (const cGH* cctkGH, const int reflevels,
-				    const vect<int,dim> ilower,
-				    const vect<int,dim> iupper,
-				    list<bbox<int,dim> >& bbl)
+				    const ivect ilower,
+				    const ivect iupper,
+				    list<ibbox>& bbl)
   {
     assert (bbl.empty());
     
     if (reflevel+1 >= reflevels) return;
     
-    const vect<int,dim> rstr = hh->baseextent.stride();
-    const vect<int,dim> rlb  = hh->baseextent.lower();
-    const vect<int,dim> rub  = hh->baseextent.upper();
+    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 vect<int,dim> str (rstr / levfac);
-    const vect<int,dim> lb  (ilower);
-    const vect<int,dim> ub  (iupper);
+    const ivect str (rstr / levfac);
+    const ivect lb  (ilower);
+    const ivect ub  (iupper);
     if (! all(lb>=rlb && ub<=rub)) {
       CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING,
 		  "The refinement region boundaries for refinement level #%d are not within the main grid", rl);
@@ -261,22 +294,22 @@ namespace CarpetRegrid {
     assert (all(lb<=ub));
     assert (all(lb%str==0 && ub%str==0));
     
-    bbl.push_back (bbox<int,dim>(lb, ub, str));
+    bbl.push_back (ibbox(lb, ub, str));
   }
   
   
   
   void MakeRegions_AsSpecified (const cGH* cctkGH, const int reflevels,
-				const vector<vect<int,dim> > lower,
-				const vector<vect<int,dim> > upper,
-				list<bbox<int,dim> >& bbl)
+				const vector<ivect> lower,
+				const vector<ivect> upper,
+				list<ibbox>& bbl)
   {
     if (reflevel+1 >= reflevels) return;
     
     const int rl = reflevel+1;
     
-    const vect<int,dim> ilower = lower[rl-1];
-    const vect<int,dim> iupper = upper[rl-1];
+    const ivect ilower = lower[rl-1];
+    const ivect iupper = upper[rl-1];
     
     MakeRegions_AsSpecified_OneLevel (cctkGH, reflevels, ilower, iupper, bbl);
   }
@@ -286,7 +319,7 @@ namespace CarpetRegrid {
   void MakeRegions_AsSpecified (const cGH* cctkGH, const int reflevels,
 				const vector<vect<CCTK_REAL,dim> > lower,
 				const vector<vect<CCTK_REAL,dim> > upper,
-				list<bbox<int,dim> >& bbl)
+				list<ibbox>& bbl)
   {
     if (reflevel+1 >= reflevels) return;
     
@@ -299,12 +332,15 @@ namespace CarpetRegrid {
 	global_upper[d] = 1;
       }
     }
-    const vect<int,dim> global_extent (hh->baseextent.upper() - hh->baseextent.lower() + hh->baseextent.stride() * (dd->lghosts + dd->ughosts));
+    const ivect global_extent (hh->baseextent.upper() - hh->baseextent.lower() + hh->baseextent.stride() * (dd->lghosts + dd->ughosts));
     
     const int rl = reflevel+1;
     
-    const vect<int,dim> ilower = vect<int,dim>(map((CCTK_REAL(*)(CCTK_REAL))(floor), (lower[rl-1] - global_lower) / (global_upper - global_lower) * vect<CCTK_REAL,dim>(global_extent) + 0.5));
-    const vect<int,dim> iupper = vect<int,dim>(map((CCTK_REAL(*)(CCTK_REAL))(floor), (upper[rl-1] - global_lower) / (global_upper - global_lower) * vect<CCTK_REAL,dim>(global_extent) + 0.5));
+    const vect<CCTK_REAL,dim> scale = vect<CCTK_REAL,dim>(global_extent) / (global_upper - global_lower);
+    const vect<CCTK_REAL,dim> rlower = (lower[rl-1] - global_lower) * scale;
+    const vect<CCTK_REAL,dim> rupper = (upper[rl-1] - global_lower) * scale;
+    const ivect ilower = ivect(map(floor, vect<double,dim>(rlower) + 0.5));
+    const ivect iupper = ivect(map(floor, vect<double,dim>(rupper) + 0.5));
     
     MakeRegions_AsSpecified_OneLevel (cctkGH, reflevels, ilower, iupper, bbl);
   }
@@ -312,10 +348,10 @@ namespace CarpetRegrid {
   
   
   static void
-  MakeRegions_Adaptively_Recombine (list<bbox<int,dim> >& bbl1,
-				    list<bbox<int,dim> >& bbl2,
-				    list<bbox<int,dim> >& bbl,
-				    const bbox<int,dim>& iface,
+  MakeRegions_Adaptively_Recombine (list<ibbox>& bbl1,
+				    list<ibbox>& bbl2,
+				    list<ibbox>& bbl,
+				    const ibbox& iface,
 				    const int dir)
   {
     assert (!iface.empty());
@@ -325,23 +361,23 @@ namespace CarpetRegrid {
     int numcombinedboxes = 0;
     
     int oldnumpoints = 0;
-    for (list<bbox<int,dim> >::const_iterator ibb = bbl.begin(); ibb != bbl.end(); ++ibb) {
+    for (list<ibbox>::const_iterator ibb = bbl.begin(); ibb != bbl.end(); ++ibb) {
       oldnumpoints += ibb->num_points();
     }
-    for (list<bbox<int,dim> >::const_iterator ibb1 = bbl1.begin(); ibb1 != bbl1.end(); ++ibb1) {
+    for (list<ibbox>::const_iterator ibb1 = bbl1.begin(); ibb1 != bbl1.end(); ++ibb1) {
       oldnumpoints += ibb1->num_points();
     }
-    for (list<bbox<int,dim> >::const_iterator ibb2 = bbl2.begin(); ibb2 != bbl2.end(); ++ibb2) {
+    for (list<ibbox>::const_iterator ibb2 = bbl2.begin(); ibb2 != bbl2.end(); ++ibb2) {
       oldnumpoints += ibb2->num_points();
     }
     
 #if 0
     // remember old bounding boxes
     bboxset<int,dim> oldboxes;
-    for (list<bbox<int,dim> >::const_iterator ibb1 = bbl1.begin(); ibb1 != bbl1.end(); ++ibb1) {
+    for (list<ibbox>::const_iterator ibb1 = bbl1.begin(); ibb1 != bbl1.end(); ++ibb1) {
       oldboxes += *ibb1;
     }
-    for (list<bbox<int,dim> >::const_iterator ibb2 = bbl2.begin(); ibb2 != bbl2.end(); ++ibb2) {
+    for (list<ibbox>::const_iterator ibb2 = bbl2.begin(); ibb2 != bbl2.end(); ++ibb2) {
       oldboxes += *ibb2;
     }
 #endif
@@ -355,18 +391,18 @@ namespace CarpetRegrid {
     cout << bbl2 << endl;
 #endif
     
-    const vect<int,dim> lo = iface.lower();
-    const vect<int,dim> up = iface.upper();
-    const vect<int,dim> str = iface.stride();
+    const ivect lo = iface.lower();
+    const ivect up = iface.upper();
+    const ivect str = iface.stride();
     
     {
       // prune boxes on the left
-      list<bbox<int,dim> >::iterator ibb1 = bbl1.begin();
+      list<ibbox>::iterator ibb1 = bbl1.begin();
       while (ibb1 != bbl1.end()) {
 	// is this bbox just to the left of the interface?
-// 	const vect<int,dim> lo1 = ibb1->lower();
-	const vect<int,dim> up1 = ibb1->upper();
-	const vect<int,dim> str1 = ibb1->stride();
+// 	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]) {
@@ -381,12 +417,12 @@ namespace CarpetRegrid {
     
     {
       // prune boxes on the right
-      list<bbox<int,dim> >::iterator ibb2 = bbl2.begin();
+      list<ibbox>::iterator ibb2 = bbl2.begin();
       while (ibb2 != bbl2.end()) {
 	// is this bbox just to the right of the interface?
-	const vect<int,dim> lo2 = ibb2->lower();
-// 	const vect<int,dim> up2 = ibb2->upper();
-	const vect<int,dim> str2 = ibb2->stride();
+	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]) {
@@ -401,31 +437,31 @@ namespace CarpetRegrid {
     
     {
       // walk all boxes on the left
-      list<bbox<int,dim> >::iterator ibb1 = bbl1.begin();
+      list<ibbox>::iterator ibb1 = bbl1.begin();
       while (ibb1 != bbl1.end()) {
-	vect<int,dim> lo1 = ibb1->lower();
-	vect<int,dim> up1 = ibb1->upper();
-	vect<int,dim> str1 = ibb1->stride();
+	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 bbox<int,dim> iface1 (lo1,up1,str1);
+	const ibbox iface1 (lo1,up1,str1);
 	
 	{
 	  // walk all boxes on the right
-	  list<bbox<int,dim> >::iterator ibb2 = bbl2.begin();
+	  list<ibbox>::iterator ibb2 = bbl2.begin();
 	  while (ibb2 != bbl2.end()) {
-	    vect<int,dim> lo2 = ibb2->lower();
-	    vect<int,dim> up2 = ibb2->upper();
-	    vect<int,dim> str2 = ibb2->stride();
+	    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 bbox<int,dim> iface2 (lo2,up2,str2);
+	    const ibbox iface2 (lo2,up2,str2);
 	    
 	    // check for a match
 	    if (iface1 == iface2) {
-	      const bbox<int,dim> combined (ibb1->lower(), ibb2->upper(), str);
+	      const ibbox combined (ibb1->lower(), ibb2->upper(), str);
 	      bbl.push_back (combined);
 	      ibb1 = bbl1.erase(ibb1);
 	      ibb2 = bbl2.erase(ibb2);
@@ -452,7 +488,7 @@ namespace CarpetRegrid {
     assert (newnumboxes + numcombinedboxes == oldnumboxes);
     
     int newnumpoints = 0;
-    for (list<bbox<int,dim> >::const_iterator ibb = bbl.begin(); ibb != bbl.end(); ++ibb) {
+    for (list<ibbox>::const_iterator ibb = bbl.begin(); ibb != bbl.end(); ++ibb) {
       newnumpoints += ibb->num_points();
     }
     assert (newnumpoints == oldnumpoints);
@@ -460,7 +496,7 @@ namespace CarpetRegrid {
 #if 0
     // find new bounding boxes
     bboxset<int,dim> newboxes;
-    for (list<bbox<int,dim> >::const_iterator ibb = bbl.begin(); ibb != bbl.end(); ++ibb) {
+    for (list<ibbox>::const_iterator ibb = bbl.begin(); ibb != bbl.end(); ++ibb) {
       newboxes += *ibb;
     }
     
@@ -482,8 +518,8 @@ namespace CarpetRegrid {
 			  const CCTK_REAL minfraction,
 			  const CCTK_REAL maxerror,
 			  const data<CCTK_REAL,dim>& error,
-			  list<bbox<int,dim> >& bbl,
-			  const bbox<int,dim>& region)
+			  list<ibbox>& bbl,
+			  const ibbox& region)
   {
     // Just to be sure
     assert (! region.empty());
@@ -492,7 +528,7 @@ namespace CarpetRegrid {
     // (this doesn't work yet on multiple processors)
     assert (CCTK_nProcs(cctkGH)==1);
     int cnt = 0;
-    for (bbox<int,dim>::iterator it=region.begin(); it!=region.end(); ++it) {
+    for (ibbox::iterator it=region.begin(); it!=region.end(); ++it) {
       if (error[*it] > maxerror) ++cnt;
     }
     const CCTK_REAL fraction = (CCTK_REAL)cnt / region.num_points();
@@ -502,30 +538,30 @@ namespace CarpetRegrid {
       // Don't refine
     } else if (width < 2*minwidth || fraction >= minfraction) {
       // Refine the whole region
-      const vect<int,dim> lo(region.lower());
-      const vect<int,dim> up(region.upper());
-      const vect<int,dim> str(region.stride());
-      bbl.push_back (bbox<int,dim>(lo,up+str-str/reffact,str/reffact));
+      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().num_points() << " fraction " << fraction << endl;
     } else {
       // Split the region and check recursively
       const int dir = maxloc(region.shape());
-      const vect<int,dim> lo(region.lower());
-      const vect<int,dim> up(region.upper());
-      const vect<int,dim> str(region.stride());
-      vect<int,dim> lo1(lo), lo2(lo);
-      vect<int,dim> up1(up), up2(up);
+      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 bbox<int,dim> region1(lo1,up1,str);
-      const bbox<int,dim> region2(lo2,up2,str);
+      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<bbox<int,dim> > bbl1, bbl2;
+      list<ibbox> bbl1, bbl2;
       MakeRegions_Adaptively (cctkGH, minwidth, minfraction, maxerror,
 			      error, bbl1, region1);
       MakeRegions_Adaptively (cctkGH, minwidth, minfraction, maxerror,
@@ -533,7 +569,7 @@ namespace CarpetRegrid {
       // Combine regions if possible
       up2 += str-str/reffact;
       up2[dir] = lo2[dir];
-      const bbox<int,dim> iface(lo2,up2,str/reffact);
+      const ibbox iface(lo2,up2,str/reffact);
       MakeRegions_Adaptively_Recombine (bbl1, bbl2, bbl, iface, dir);
     }
     
@@ -545,7 +581,7 @@ namespace CarpetRegrid {
 			  const CCTK_REAL minfraction,
 			  const CCTK_REAL maxerror,
 			  const gf<CCTK_REAL,dim>& error,
-			  list<bbox<int,dim> >& bbl)
+			  list<ibbox>& bbl)
   {
     assert (bbl.empty());
     
@@ -562,7 +598,7 @@ namespace CarpetRegrid {
 //     cout << endl << "MRA: Choosing regions to refine in " << hh->components(rl) << " components" << endl;
     
     for (int c=0; c<hh->components(rl); ++c) {
-      const bbox<int,dim> region = hh->extents[rl][c][ml];
+      const ibbox region = hh->extents[rl][c][ml];
       assert (! region.empty());
       
       const data<CCTK_REAL,dim>& errdata = *error(tl,rl,c,ml);
@@ -571,10 +607,10 @@ namespace CarpetRegrid {
 			      errdata, bbl, region);
     }
     
-    int numpoints = 0;
-    for (list<bbox<int,dim> >::const_iterator ibb = bbl.begin(); ibb != bbl.end(); ++ibb) {
-      numpoints += ibb->num_points();
-    }
+//     int numpoints = 0;
+//     for (list<ibbox>::const_iterator ibb = bbl.begin(); ibb != bbl.end(); ++ibb) {
+//       numpoints += ibb->num_points();
+//     }
 //     cout << "MRA: Chose " << bbl.size() << " regions with a total size of " << numpoints << " to refine." << endl << endl;
   }
   
diff --git a/Carpet/CarpetRegrid/src/regrid.hh b/Carpet/CarpetRegrid/src/regrid.hh
index cef62a2fb..402e4f8df 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.3 2002/01/11 17:37:14 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/src/regrid.hh,v 1.4 2002/03/11 13:17:16 schnetter Exp $
 
 #ifndef REGRID_HH
 #define REGRID_HH
@@ -28,6 +28,7 @@ namespace CarpetRegrid {
   
   int CarpetRegridRegrid (const cGH * const cctkGH,
 			  gh<dim>::rexts& bbsss,
+			  gh<dim>::rbnds& obss,
 			  gh<dim>::rprocs& pss);
   
   
diff --git a/CarpetExtra/WaveToyF77/par/wavetoyf77_rad_full_coarse.par b/CarpetExtra/WaveToyF77/par/wavetoyf77_rad_full_coarse.par
index 9c0093d30..a882a08a8 100644
--- a/CarpetExtra/WaveToyF77/par/wavetoyf77_rad_full_coarse.par
+++ b/CarpetExtra/WaveToyF77/par/wavetoyf77_rad_full_coarse.par
@@ -7,7 +7,7 @@
 #   @enddesc
 # @@*/
 #
-# $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/WaveToyF77/par/wavetoyf77_rad_full_coarse.par,v 1.3 2002/01/30 16:07:31 schnetter Exp $
+# $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/WaveToyF77/par/wavetoyf77_rad_full_coarse.par,v 1.4 2002/03/11 13:17:16 schnetter Exp $
 
 ActiveThorns = "Boundary IOBasic IOUtil Time Cart3d Carpet CarpetIOASCII CarpetLib CarpetReduce CarpetRegrid CarpetSlab IDScalarWave WaveToyF77"
 
@@ -25,13 +25,13 @@ grid::avoid_origin		= no
 
 IO::outdir			= "wavetoyf77_rad_full_coarse"
 
-IOBasic::outinfo_every		= 1 # 10
-#IOBasic::outinfo_vars		= "wavetoy::phi"
+IOBasic::outinfo_every		= 1
+IOBasic::outinfo_vars		= "wavetoy::phi"
 
-#IOBasic::outScalar_every	= 2
-#IOBasic::outScalar_vars		= "wavetoy::phi"
+IOBasic::outScalar_every	= 2
+IOBasic::outScalar_vars		= "wavetoy::phi"
 
-IOASCII::out1D_every		= 1 # 2
+IOASCII::out1D_every		= 1
 IOASCII::out1D_vars		= "wavetoy::phi grid::coordinates"
 
 WaveToyF77::bound		= radiation
diff --git a/CarpetExtra/WaveToyF77/par/wavetoyf77_rad_full_coarse_rl2.par b/CarpetExtra/WaveToyF77/par/wavetoyf77_rad_full_coarse_rl2.par
index 61dd57386..301944cc0 100644
--- a/CarpetExtra/WaveToyF77/par/wavetoyf77_rad_full_coarse_rl2.par
+++ b/CarpetExtra/WaveToyF77/par/wavetoyf77_rad_full_coarse_rl2.par
@@ -7,11 +7,11 @@
 #   @enddesc
 # @@*/
 #
-# $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/WaveToyF77/par/wavetoyf77_rad_full_coarse_rl2.par,v 1.3 2002/01/30 16:07:31 schnetter Exp $
+# $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/WaveToyF77/par/wavetoyf77_rad_full_coarse_rl2.par,v 1.4 2002/03/11 13:17:16 schnetter Exp $
 
 ActiveThorns = "Boundary IOBasic IOUtil Time Cart3d Carpet CarpetIOASCII CarpetLib CarpetReduce CarpetRegrid CarpetSlab IDScalarWave WaveToyF77"
 
-Cactus::cctk_itlast		= 60
+Cactus::cctk_itlast		= 120
  
 Time::dtfac			= 0.5
 
@@ -21,19 +21,21 @@ driver::global_nz		= 31
 
 Carpet::max_refinement_levels	= 2
 
+CarpetRegrid::refinement_levels	= 2
+
 grid::type			= byspacing
 grid::dxyz			= 0.6
 grid::avoid_origin		= no
 
 IO::outdir			= "wavetoyf77_rad_full_coarse_rl2"
 
-IOBasic::outinfo_every		= 1 # 10
-#IOBasic::outinfo_vars		= "wavetoy::phi"
+IOBasic::outinfo_every		= 1
+IOBasic::outinfo_vars		= "wavetoy::phi"
 
-#IOBasic::outScalar_every	= 2
-#IOBasic::outScalar_vars		= "wavetoy::phi"
+IOBasic::outScalar_every	= 2
+IOBasic::outScalar_vars		= "wavetoy::phi"
 
-IOASCII::out1D_every		= 1 # 2
+IOASCII::out1D_every		= 1
 IOASCII::out1D_vars		= "wavetoy::phi grid::coordinates"
 
 WaveToyF77::bound		= radiation
diff --git a/CarpetExtra/WaveToyF77/par/wavetoyf77_rad_full_coarse_rl3.par b/CarpetExtra/WaveToyF77/par/wavetoyf77_rad_full_coarse_rl3.par
index c30bf6446..7ef9ea780 100644
--- a/CarpetExtra/WaveToyF77/par/wavetoyf77_rad_full_coarse_rl3.par
+++ b/CarpetExtra/WaveToyF77/par/wavetoyf77_rad_full_coarse_rl3.par
@@ -7,11 +7,11 @@
 #   @enddesc
 # @@*/
 #
-# $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/WaveToyF77/par/wavetoyf77_rad_full_coarse_rl3.par,v 1.3 2002/01/30 16:07:31 schnetter Exp $
+# $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/WaveToyF77/par/wavetoyf77_rad_full_coarse_rl3.par,v 1.4 2002/03/11 13:17:16 schnetter Exp $
 
 ActiveThorns = "Boundary IOBasic IOUtil Time Cart3d Carpet CarpetIOASCII CarpetLib CarpetReduce CarpetRegrid CarpetSlab IDScalarWave WaveToyF77"
 
-Cactus::cctk_itlast		= 60
+Cactus::cctk_itlast		= 240
  
 Time::dtfac			= 0.5
 
@@ -21,19 +21,21 @@ driver::global_nz		= 31
 
 Carpet::max_refinement_levels	= 3
 
+CarpetRegrid::refinement_levels	= 3
+
 grid::type			= byspacing
 grid::dxyz			= 0.6
 grid::avoid_origin		= no
 
 IO::outdir			= "wavetoyf77_rad_full_coarse_rl3"
 
-IOBasic::outinfo_every		= 1 # 10
-#IOBasic::outinfo_vars		= "wavetoy::phi"
+IOBasic::outinfo_every		= 1
+IOBasic::outinfo_vars		= "wavetoy::phi"
 
-#IOBasic::outScalar_every	= 2
-#IOBasic::outScalar_vars		= "wavetoy::phi"
+IOBasic::outScalar_every	= 2
+IOBasic::outScalar_vars		= "wavetoy::phi"
 
-IOASCII::out1D_every		= 1 # 2
+IOASCII::out1D_every		= 1
 IOASCII::out1D_vars		= "wavetoy::phi grid::coordinates"
 
 WaveToyF77::bound		= radiation
diff --git a/CarpetExtra/WaveToyF77/par/wavetoyf77_zero_ell_coarse.par b/CarpetExtra/WaveToyF77/par/wavetoyf77_zero_ell_coarse.par
index 66121b671..aa6fd3ce7 100644
--- a/CarpetExtra/WaveToyF77/par/wavetoyf77_zero_ell_coarse.par
+++ b/CarpetExtra/WaveToyF77/par/wavetoyf77_zero_ell_coarse.par
@@ -7,7 +7,7 @@
 #   @enddesc
 # @@*/
 #
-# $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/WaveToyF77/par/wavetoyf77_zero_ell_coarse.par,v 1.3 2004/03/23 12:00:37 schnetter Exp $
+# $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/WaveToyF77/par/wavetoyf77_zero_ell_coarse.par,v 1.1 2002/03/11 13:17:17 schnetter Exp $
 
 ActiveThorns = "Boundary IOBasic IOUtil Time Cart3d Carpet CarpetIOASCII CarpetLib CarpetReduce CarpetRegrid CarpetSlab IDScalarWave WaveToyF77"
 
@@ -26,7 +26,7 @@ grid::type			= byspacing
 grid::dxyz			= 0.6
 grid::avoid_origin		= no
 
-IO::out_dir			= "wavetoyf77_zero_ell_coarse"
+IO::outdir			= "wavetoyf77_zero_ell_coarse"
 
 IOBasic::outinfo_every		= 1
 IOBasic::outinfo_vars		= "wavetoy::phi"
diff --git a/CarpetExtra/WaveToyF77/par/wavetoyf77_zero_ell_coarse_rl2.par b/CarpetExtra/WaveToyF77/par/wavetoyf77_zero_ell_coarse_rl2.par
index bf73eafc9..21871b67b 100644
--- a/CarpetExtra/WaveToyF77/par/wavetoyf77_zero_ell_coarse_rl2.par
+++ b/CarpetExtra/WaveToyF77/par/wavetoyf77_zero_ell_coarse_rl2.par
@@ -7,7 +7,7 @@
 #   @enddesc
 # @@*/
 #
-# $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/WaveToyF77/par/wavetoyf77_zero_ell_coarse_rl2.par,v 1.3 2004/03/23 12:00:37 schnetter Exp $
+# $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/WaveToyF77/par/wavetoyf77_zero_ell_coarse_rl2.par,v 1.1 2002/03/11 13:17:17 schnetter Exp $
 
 ActiveThorns = "Boundary IOBasic IOUtil Time Cart3d Carpet CarpetIOASCII CarpetLib CarpetReduce CarpetRegrid CarpetSlab IDScalarWave WaveToyF77"
 
@@ -30,7 +30,7 @@ grid::type			= byspacing
 grid::dxyz			= 0.6
 grid::avoid_origin		= no
 
-IO::out_dir			= "wavetoyf77_zero_ell_coarse_rl2"
+IO::outdir			= "wavetoyf77_zero_ell_coarse_rl2"
 
 IOBasic::outinfo_every		= 1
 IOBasic::outinfo_vars		= "wavetoy::phi"