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

23 files changed, 414 insertions, 310 deletions

diff --git a/Carpet/CarpetLib/interface.ccl b/Carpet/CarpetLib/interface.ccl
index b2c8cbde0..717c5ec06 100644
--- a/Carpet/CarpetLib/interface.ccl
+++ b/Carpet/CarpetLib/interface.ccl
@@ -1,5 +1,5 @@
 # Interface definition for thorn CarpetLib
-# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/interface.ccl,v 1.3 2003/06/18 18:28:07 schnetter Exp $
+# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/interface.ccl,v 1.4 2004/01/25 14:57:29 schnetter Exp $
 
 implements: CarpetLib
 
@@ -18,3 +18,5 @@ includes header: gf.hh in gf.hh
 includes header: ggf.hh in ggf.hh
 includes header: gh.hh in gh.hh
 includes header: th.hh in th.hh
+
+includes header: operators.hh in operators.hh
diff --git a/Carpet/CarpetLib/src/bbox.cc b/Carpet/CarpetLib/src/bbox.cc
index bbfcc3df5..6decccabe 100644
--- a/Carpet/CarpetLib/src/bbox.cc
+++ b/Carpet/CarpetLib/src/bbox.cc
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/bbox.cc,v 1.18 2003/11/13 16:03:58 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/bbox.cc,v 1.19 2004/01/25 14:57:29 schnetter Exp $
 
 #include <assert.h>
 
@@ -249,20 +249,16 @@ typename bbox<T,D>::iteratorT bbox<T,D>::endT () const {
 template<class T,int D>
 void bbox<T,D>::input (istream& is) {
   skipws (is);
-  assert (is.peek() == '(');
-  is.get();
+  consume (is, '(');
   is >> _lower;
   skipws (is);
-  assert (is.peek() == ':');
-  is.get();
+  consume (is, ':');
   is >> _upper;
   skipws (is);
-  assert (is.peek() == ':');
-  is.get();
+  consume (is, ':');
   is >> _stride;
   skipws (is);
-  assert (is.peek() == ')');
-  is.get();
+  consume (is, ')');
 }
 
 
diff --git a/Carpet/CarpetLib/src/data.cc b/Carpet/CarpetLib/src/data.cc
index f6158bae8..6f65158b8 100644
--- a/Carpet/CarpetLib/src/data.cc
+++ b/Carpet/CarpetLib/src/data.cc
@@ -1,10 +1,15 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/data.cc,v 1.38 2004/01/21 16:32:04 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/data.cc,v 1.39 2004/01/25 14:57:29 schnetter Exp $
 
 #include <assert.h>
 #include <limits.h>
+#include <stdlib.h>
+#include <math.h>
 
+#include <algorithm>
 #include <iostream>
+#include <sstream>
 #include <string>
+#include <vector>
 
 #include <mpi.h>
 
@@ -34,16 +39,17 @@ static int nexttag ()
 
 // Constructors
 template<class T, int D>
-data<T,D>::data (const int varindex_)
-  : gdata<D>(varindex_),
+data<T,D>::data (const int varindex_, const operator_type transport_operator_)
+  : gdata<D>(varindex_, transport_operator_),
     _storage(0),
     comm_active(false),
     tag(nexttag())
 { }
 
 template<class T, int D>
-data<T,D>::data (const int varindex_, const ibbox& extent_, const int proc_)
-  : gdata<D>(varindex_),
+data<T,D>::data (const int varindex_, const operator_type transport_operator_,
+                 const ibbox& extent_, const int proc_)
+  : gdata<D>(varindex_, transport_operator_),
     _storage(0),
     comm_active(false),
     tag(nexttag())
@@ -59,8 +65,11 @@ data<T,D>::~data () {
   
 // Pseudo constructors
 template<class T, int D>
-data<T,D>* data<T,D>::make_typed (const int varindex_) const {
-  return new data(varindex_);
+data<T,D>* data<T,D>::make_typed (const int varindex_,
+                                  const operator_type transport_operator_)
+  const
+{
+  return new data(varindex_, transport_operator_);
 }
 
 
@@ -107,7 +116,7 @@ void data<T,D>::transfer_from (gdata<D>* gsrc) {
   data* src = (data*)gsrc;
   assert (!_storage);
   *this = *src;
-  *src = data(this->varindex);
+  *src = data(this->varindex, this->transport_operator);
 }
 
 
@@ -159,9 +168,12 @@ void data<T,D>::change_processor_recv (const int newproc, void* const mem)
 	_storage = (T*)mem;
       }
       
+      const double wtime1 = MPI_Wtime();
       T dummy;
       MPI_Irecv (_storage, this->_size, dist::datatype(dummy), this->_proc,
                  this->tag, dist::comm, &request);
+      const double wtime2 = MPI_Wtime();
+      this->wtime_irecv += wtime2 - wtime1;
       
     } else if (rank == this->_proc) {
       // copy to other processor
@@ -197,9 +209,12 @@ void data<T,D>::change_processor_send (const int newproc, void* const mem)
       assert (!mem);
       assert (_storage);
       
+      const double wtime1 = MPI_Wtime();
       T dummy;
       MPI_Isend (_storage, this->_size, dist::datatype(dummy), newproc,
                  this->tag, dist::comm, &request);
+      const double wtime2 = MPI_Wtime();
+      this->wtime_isend += wtime2 - wtime1;
       
     } else {
       assert (!mem);
@@ -227,8 +242,11 @@ void data<T,D>::change_processor_wait (const int newproc, void* const mem)
     if (rank == newproc) {
       // copy from other processor
       
+      const double wtime1 = MPI_Wtime();
       MPI_Status status;
       MPI_Wait (&request, &status);
+      const double wtime2 = MPI_Wtime();
+      this->wtime_irecvwait += wtime2 - wtime1;
       
     } else if (rank == this->_proc) {
       // copy to other processor
@@ -236,8 +254,11 @@ void data<T,D>::change_processor_wait (const int newproc, void* const mem)
       assert (!mem);
       assert (_storage);
       
+      const double wtime1 = MPI_Wtime();
       MPI_Status status;
       MPI_Wait (&request, &status);
+      const double wtime2 = MPI_Wtime();
+      this->wtime_isendwait += wtime2 - wtime1;
       
       if (this->_owns_storage) {
 	delete [] _storage;
@@ -553,6 +574,15 @@ extern "C" {
      const int srcbbox[3][3],
      const int dstbbox[3][3],
      const int regbbox[3][3]);
+  void CCTK_FCALL CCTK_FNAME(prolongate_3d_real8_2tl_rf2)
+    (const CCTK_REAL8* src1, const CCTK_REAL8& t1,
+     const CCTK_REAL8* src2, const CCTK_REAL8& t2,
+     const int& srciext, const int& srcjext, const int& srckext,
+     CCTK_REAL8* dst, const CCTK_REAL8& t,
+     const int& dstiext, const int& dstjext, const int& dstkext,
+     const int srcbbox[3][3],
+     const int dstbbox[3][3],
+     const int regbbox[3][3]);
   void CCTK_FCALL CCTK_FNAME(prolongate_3d_real8_2tl_o3)
     (const CCTK_REAL8* src1, const CCTK_REAL8& t1,
      const CCTK_REAL8* src2, const CCTK_REAL8& t2,
@@ -562,6 +592,15 @@ extern "C" {
      const int srcbbox[3][3],
      const int dstbbox[3][3],
      const int regbbox[3][3]);
+  void CCTK_FCALL CCTK_FNAME(prolongate_3d_real8_2tl_o3_rf2)
+    (const CCTK_REAL8* src1, const CCTK_REAL8& t1,
+     const CCTK_REAL8* src2, const CCTK_REAL8& t2,
+     const int& srciext, const int& srcjext, const int& srckext,
+     CCTK_REAL8* dst, const CCTK_REAL8& t,
+     const int& dstiext, const int& dstjext, const int& dstkext,
+     const int srcbbox[3][3],
+     const int dstbbox[3][3],
+     const int regbbox[3][3]);
   void CCTK_FCALL CCTK_FNAME(prolongate_3d_real8_2tl_minmod)
     (const CCTK_REAL8* src1, const CCTK_REAL8& t1,
      const CCTK_REAL8* src2, const CCTK_REAL8& t2,
@@ -651,6 +690,8 @@ void data<CCTK_REAL8,3>
 			      const int order_space,
 			      const int order_time)
 {
+  const CCTK_REAL eps = 1.0e-10;
+  
   assert (has_storage());
   assert (all(box.lower()>=extent().lower()));
   assert (all(box.upper()<=extent().upper()));
@@ -694,8 +735,46 @@ void data<CCTK_REAL8,3>
     regbbox[2][d] = box.stride()[d];
   }
   
+  // Check that the times are consistent
+  assert (times.size() > 0);
+  CCTK_REAL min_time = times[0];
+  CCTK_REAL max_time = times[0];
+  for (size_t tl=1; tl<times.size(); ++tl) {
+    min_time = min(min_time, times[tl]);
+    max_time = max(max_time, times[tl]);
+  }
+  if (time < min_time - eps || time > max_time + eps) {
+    ostringstream buf;
+    buf << "Internal error: extrapolation in time."
+        << "  time=" << time
+        << "  times=" << times;
+    CCTK_WARN (0, buf.str().c_str());
+  }
+  
+  // Is it necessary to interpolate in time?
+  if (times.size() > 1) {
+    for (size_t tl=0; tl<times.size(); ++tl) {
+      assert (abs(1.5) > 1.4);
+      if (abs(times[tl] - time) < eps) {
+        // It is not.
+        vector<const gdata<3>*> my_gsrcs(1);
+        vector<CCTK_REAL> my_times(1);
+        my_gsrcs[0] = gsrcs[tl];
+        my_times[0] = times[tl];
+        const int my_order_time = 0;
+        this->interpolate_from_innerloop
+          (my_gsrcs, my_times, box, time, order_space, my_order_time);
+        return;
+      }
+    }
+  }
+  
   assert (all(dext.stride() == box.stride()));
   if (all(sext.stride() < dext.stride())) {
+    // Restrict
+    
+    assert (times.size() == 1);
+    assert (abs(times[0] - time) < eps);
     
     switch (transport_operator) {
       
@@ -722,9 +801,10 @@ void data<CCTK_REAL8,3>
     default:
       assert (0);
       
-    } // switch (prolong_method)
+    } // switch (transport_operator)
     
   } else if (all(sext.stride() > dext.stride())) {
+    // Prolongate
     
     switch (transport_operator) {
       
@@ -732,6 +812,8 @@ void data<CCTK_REAL8,3>
       switch (order_time) {
         
       case 0:
+        assert (times.size() == 1);
+        assert (abs(times[0] - time) < eps);
         assert (srcs.size()>=1);
         switch (order_space) {
         case 0:
@@ -789,23 +871,43 @@ void data<CCTK_REAL8,3>
         switch (order_space) {
         case 0:
         case 1:
-          CCTK_FNAME(prolongate_3d_real8_2tl)
-            ((const CCTK_REAL8*)srcs[0]->storage(), times[0],
-             (const CCTK_REAL8*)srcs[1]->storage(), times[1],
-             srcshp[0], srcshp[1], srcshp[2],
-             (CCTK_REAL8*)storage(), time,
-             dstshp[0], dstshp[1], dstshp[2],
-             srcbbox, dstbbox, regbbox);
+          if (all (sext.stride() == dext.stride() * 2)) {
+            CCTK_FNAME(prolongate_3d_real8_2tl_rf2)
+              ((const CCTK_REAL8*)srcs[0]->storage(), times[0],
+               (const CCTK_REAL8*)srcs[1]->storage(), times[1],
+               srcshp[0], srcshp[1], srcshp[2],
+               (CCTK_REAL8*)storage(), time,
+               dstshp[0], dstshp[1], dstshp[2],
+               srcbbox, dstbbox, regbbox);
+          } else {
+            CCTK_FNAME(prolongate_3d_real8_2tl)
+              ((const CCTK_REAL8*)srcs[0]->storage(), times[0],
+               (const CCTK_REAL8*)srcs[1]->storage(), times[1],
+               srcshp[0], srcshp[1], srcshp[2],
+               (CCTK_REAL8*)storage(), time,
+               dstshp[0], dstshp[1], dstshp[2],
+               srcbbox, dstbbox, regbbox);
+          }
           break;
         case 2:
         case 3:
-          CCTK_FNAME(prolongate_3d_real8_2tl_o3)
-            ((const CCTK_REAL8*)srcs[0]->storage(), times[0],
-             (const CCTK_REAL8*)srcs[1]->storage(), times[1],
-             srcshp[0], srcshp[1], srcshp[2],
-             (CCTK_REAL8*)storage(), time,
-             dstshp[0], dstshp[1], dstshp[2],
-             srcbbox, dstbbox, regbbox);
+          if (all (sext.stride() == dext.stride() * 2)) {
+            CCTK_FNAME(prolongate_3d_real8_2tl_o3_rf2)
+              ((const CCTK_REAL8*)srcs[0]->storage(), times[0],
+               (const CCTK_REAL8*)srcs[1]->storage(), times[1],
+               srcshp[0], srcshp[1], srcshp[2],
+               (CCTK_REAL8*)storage(), time,
+               dstshp[0], dstshp[1], dstshp[2],
+               srcbbox, dstbbox, regbbox);
+          } else {
+            CCTK_FNAME(prolongate_3d_real8_2tl_o3)
+              ((const CCTK_REAL8*)srcs[0]->storage(), times[0],
+               (const CCTK_REAL8*)srcs[1]->storage(), times[1],
+               srcshp[0], srcshp[1], srcshp[2],
+               (CCTK_REAL8*)storage(), time,
+               dstshp[0], dstshp[1], dstshp[2],
+               srcbbox, dstbbox, regbbox);
+          }
           break;
         case 4:
         case 5:
@@ -893,40 +995,75 @@ void data<CCTK_REAL8,3>
     case op_TVD:
       switch (order_time) {
       case 0: 
-        CCTK_FNAME(prolongate_3d_real8_minmod)
-          ((const CCTK_REAL8*)srcs[0]->storage(),
-           srcshp[0], srcshp[1], srcshp[2],
-           (CCTK_REAL8*)storage(),
-           dstshp[0], dstshp[1], dstshp[2],
-           srcbbox, dstbbox, regbbox);
+        assert (times.size() == 1);
+        assert (abs(times[0] - time) < eps);
+        switch (order_space) {
+        case 0:
+        case 1:
+          CCTK_WARN (0, "There is no stencil for op=""TVD"" with order_space=1");
+          break;
+        case 2:
+        case 3:
+          CCTK_FNAME(prolongate_3d_real8_minmod)
+            ((const CCTK_REAL8*)srcs[0]->storage(),
+             srcshp[0], srcshp[1], srcshp[2],
+             (CCTK_REAL8*)storage(),
+             dstshp[0], dstshp[1], dstshp[2],
+             srcbbox, dstbbox, regbbox);
+          break;
+        default:
+          assert (0);
+        }
         break;
       case 1:
-        CCTK_FNAME(prolongate_3d_real8_2tl_minmod)
-          ((const CCTK_REAL8*)srcs[0]->storage(), times[0],
-           (const CCTK_REAL8*)srcs[1]->storage(), times[1],
-           srcshp[0], srcshp[1], srcshp[2],
-           (CCTK_REAL8*)storage(), time,
-           dstshp[0], dstshp[1], dstshp[2],
-           srcbbox, dstbbox, regbbox);
+        switch (order_space) {
+        case 0:
+        case 1:
+          CCTK_WARN (0, "There is no stencil for op=""TVD"" with order_space=1");
+          break;
+        case 2:
+        case 3:
+          CCTK_FNAME(prolongate_3d_real8_2tl_minmod)
+            ((const CCTK_REAL8*)srcs[0]->storage(), times[0],
+             (const CCTK_REAL8*)srcs[1]->storage(), times[1],
+             srcshp[0], srcshp[1], srcshp[2],
+             (CCTK_REAL8*)storage(), time,
+             dstshp[0], dstshp[1], dstshp[2],
+             srcbbox, dstbbox, regbbox);
+          break;
+        default:
+          assert (0);
+        }
         break;
       case 2: 
-        CCTK_FNAME(prolongate_3d_real8_3tl_minmod)
-          ((const CCTK_REAL8*)srcs[0]->storage(), times[0],
-           (const CCTK_REAL8*)srcs[1]->storage(), times[1],
-           (const CCTK_REAL8*)srcs[2]->storage(), times[2],
-           srcshp[0], srcshp[1], srcshp[2],
-           (CCTK_REAL8*)storage(), time,
-           dstshp[0], dstshp[1], dstshp[2],
-           srcbbox, dstbbox, regbbox);
+        switch (order_space) {
+        case 0:
+        case 1:
+          CCTK_WARN (0, "There is no stencil for op=""TVD"" with order_space=1");
+          break;
+        case 2:
+        case 3:
+          CCTK_FNAME(prolongate_3d_real8_3tl_minmod)
+            ((const CCTK_REAL8*)srcs[0]->storage(), times[0],
+             (const CCTK_REAL8*)srcs[1]->storage(), times[1],
+             (const CCTK_REAL8*)srcs[2]->storage(), times[2],
+             srcshp[0], srcshp[1], srcshp[2],
+             (CCTK_REAL8*)storage(), time,
+             dstshp[0], dstshp[1], dstshp[2],
+             srcbbox, dstbbox, regbbox);
+          break;
+        default:
+          assert (0);
+        }
         break;
       default:
         assert (0);
-      }
+      } // switch (order_time)
       break;
       
     default:
       assert(0);
-    } // switch (prolong_method)
+    } // switch (transport_operator)
     
   } else {
     assert (0);
diff --git a/Carpet/CarpetLib/src/data.hh b/Carpet/CarpetLib/src/data.hh
index bb8e5c497..016dfcd1e 100644
--- a/Carpet/CarpetLib/src/data.hh
+++ b/Carpet/CarpetLib/src/data.hh
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/data.hh,v 1.15 2003/11/05 16:18:39 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/data.hh,v 1.16 2004/01/25 14:57:29 schnetter Exp $
 
 #ifndef DATA_HH
 #define DATA_HH
@@ -39,14 +39,17 @@ class data: public gdata<D> {
 public:
   
   // Constructors
-  data (const int varindex);
-  data (const int varindex, const ibbox& extent, const int proc);
+  data (const int varindex = -1,
+        const operator_type transport_operator = op_error);
+  data (const int varindex, const operator_type transport_operator,
+        const ibbox& extent, const int proc);
 
   // Destructors
   virtual ~data ();
 
   // Pseudo constructors
-  virtual data* make_typed (const int varindex) const;
+  virtual data* make_typed (const int varindex,
+                            const operator_type transport_operator) const;
 
   // Storage management
   virtual void allocate (const ibbox& extent, const int proc,
diff --git a/Carpet/CarpetLib/src/defs.cc b/Carpet/CarpetLib/src/defs.cc
index b525578ef..28aaf0539 100644
--- a/Carpet/CarpetLib/src/defs.cc
+++ b/Carpet/CarpetLib/src/defs.cc
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/defs.cc,v 1.17 2003/11/13 16:03:58 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/defs.cc,v 1.18 2004/01/25 14:57:29 schnetter Exp $
 
 #include <assert.h>
 #include <ctype.h>
@@ -9,12 +9,31 @@
 #include <stack>
 #include <vector>
 
+#include "cctk.h"
+
 #include "defs.hh"
 
 using namespace std;
 
 
 
+template<class T>
+T ipow (T x, int y) {
+  if (y<0) {
+    y = -y;
+    x = T(1)/x;
+  }
+  T res = T(1);
+  while (y>0) {
+    if (y%2) res *= x;
+    x *= x;
+    y /= 2;
+  }
+  return res;
+}
+
+
+
 void skipws (istream& is) {
   while (is.good() && isspace(is.peek())) {
     is.get();
@@ -23,26 +42,53 @@ void skipws (istream& is) {
 
 
 
+void expect (istream& is, const char c) {
+  if (is.peek() == c) return;
+  cout << "While reading characters from a stream:" << endl
+       << "   Character '" << c << "' expected, but not found." << endl
+       << "   The next up to 100 available characters are \"";
+  for (int i=0; i<100; ++i) {
+    const int uc = is.get();
+    if (uc<0) break;
+    cout << (unsigned char)uc;
+  }
+  cout << "\"." << endl;
+  throw input_error();
+}
+
+
+
+void consume (istream& is, const char c) {
+  expect (is, c);
+  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() != ']') {
-    T elem;
-    is >> elem;
-    v.push_back (elem);
+  try {
     skipws (is);
-    if (is.peek() != ',') break;
-    is.get();
+    consume (is, '[');
+    skipws (is);
+    while (is.good() && is.peek() != ']') {
+      T elem;
+      is >> elem;
+      v.push_back (elem);
+      skipws (is);
+      if (is.peek() != ',') break;
+      is.get();
+      skipws (is);
+    }
     skipws (is);
+    consume (is, ']');
+  } catch (input_error &err) {
+    cout << "Input error while reading a vector<>" << endl
+         << "   The following elements have been read so far: " << v << endl;
+    throw err;
   }
-  skipws (is);
-  assert (is.peek() == ']');
-  is.get();
   return is;
 }
 
@@ -102,10 +148,12 @@ ostream& output (ostream& os, const vector<T>& v) {
 #include "bbox.hh"
 #include "bboxset.hh"
 
+template int ipow (int x, int y);
+
 template istream& input (istream& os, vector<bbox<int,3> >& v);
-template istream& input (istream& os, vector<bbox<double,3> >& v);
+template istream& input (istream& os, vector<bbox<CCTK_REAL,3> >& v);
 template istream& input (istream& os, vector<vector<bbox<int,3> > >& v);
-template istream& input (istream& os, vector<vector<bbox<double,3> > >& v);
+template istream& input (istream& os, vector<vector<bbox<CCTK_REAL,3> > >& v);
 template istream& input (istream& os, vector<vect<vect<bool,2>,3> >& v);
 template istream& input (istream& os, vector<vector<vect<vect<bool,2>,3> > >& v);
 
@@ -114,12 +162,14 @@ template ostream& output (ostream& os, const set<bbox<int,3> >& s);
 template ostream& output (ostream& os, const set<bboxset<int,3> >& s);
 template ostream& output (ostream& os, const stack<bbox<int,3> >& s);
 template ostream& output (ostream& os, const vector<int>& v);
+template ostream& output (ostream& os, const vector<CCTK_REAL>& v);
 template ostream& output (ostream& os, const vector<bbox<int,3> >& v);
-template ostream& output (ostream& os, const vector<bbox<double,3> >& v);
+template ostream& output (ostream& os, const vector<bbox<CCTK_REAL,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<CCTK_REAL> >& v);
 template ostream& output (ostream& os, const vector<vector<bbox<int,3> > >& v);
-template ostream& output (ostream& os, const vector<vector<bbox<double,3> > >& v);
+template ostream& output (ostream& os, const vector<vector<bbox<CCTK_REAL,3> > >& v);
 template ostream& output (ostream& os, const vector<vector<vect<vect<bool,2>,3> > >& v);
 template ostream& output (ostream& os, const vector<vect<vect<bool,2>,3> >& v);
 template ostream& output (ostream& os, const vector<vector<vector<bbox<int,3> > > >& v);
diff --git a/Carpet/CarpetLib/src/defs.hh b/Carpet/CarpetLib/src/defs.hh
index b948d7a9a..2d324032a 100644
--- a/Carpet/CarpetLib/src/defs.hh
+++ b/Carpet/CarpetLib/src/defs.hh
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/defs.hh,v 1.11 2003/03/18 17:30:25 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/defs.hh,v 1.12 2004/01/25 14:57:29 schnetter Exp $
 
 #ifndef DEFS_HH
 #define DEFS_HH
@@ -34,24 +34,15 @@ inline T square (const T& x) { return x*x; }
 
 // Another useful helper
 template<class T>
-inline T ipow (T x, int y) {
-  if (y<0) {
-    y = -y;
-    x = T(1)/x;
-  }
-  T res = T(1);
-  while (y>0) {
-    if (y%2) res *= x;
-    x *= x;
-    y /= 2;
-  }
-  return res;
-}
+T ipow (T x, int y);
 
 
 
-// Skip whitespace
+// Input streams
+struct input_error { };
 void skipws (istream& is);
+void expect (istream& is, const char c);
+void consume (istream& is, const char c);
 
 
 
diff --git a/Carpet/CarpetLib/src/dh.cc b/Carpet/CarpetLib/src/dh.cc
index 3cb9ac1c0..e23268791 100644
--- a/Carpet/CarpetLib/src/dh.cc
+++ b/Carpet/CarpetLib/src/dh.cc
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/dh.cc,v 1.47 2004/01/21 14:25:35 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/dh.cc,v 1.48 2004/01/25 14:57:29 schnetter Exp $
 
 #include <assert.h>
 
@@ -127,7 +127,8 @@ void dh<D>::recompose (const int initialise_from, const bool do_prolongate) {
       	for (int cc=0; cc<h.components(rl); ++cc) {
       	  assert (ml<h.mglevels(rl,cc));
       	  // intersect boundaries with interior of that component
-      	  const ibset ovlp = bnds & boxes.at(rl).at(cc).at(ml).interior;
+          ibset ovlp = bnds & boxes.at(rl).at(cc).at(ml).interior;
+          ovlp.normalize();
       	  for (typename ibset::const_iterator b=ovlp.begin();
 	       b!=ovlp.end(); ++b) {
 	    boxes.at(rl).at(c ).at(ml).recv_sync.at(cc).push_back(*b);
@@ -179,7 +180,6 @@ void dh<D>::recompose (const int initialise_from, const bool do_prolongate) {
     } // for c
   } // for rl
   
-  // TODO: prefer boxes from the same processor
   for (int rl=0; rl<h.reflevels(); ++rl) {
     for (int c=0; c<h.components(rl); ++c) {
       for (int ml=0; ml<h.mglevels(rl,c); ++ml) {
@@ -191,6 +191,7 @@ void dh<D>::recompose (const int initialise_from, const bool do_prolongate) {
       	  for (int cc=0; cc<h.components(rl+1); ++cc) {
       	    const ibbox intrf = boxes.at(rl+1).at(cc).at(ml).interior;
       	    // Prolongation (interior)
+            // TODO: prefer boxes from the same processor
       	    {
       	      // (the prolongation may use the exterior of the coarse
       	      // grid, and must fill all of the interior of the fine
@@ -214,8 +215,10 @@ void dh<D>::recompose (const int initialise_from, const bool do_prolongate) {
                 const ibbox send = recv.expanded_for(extr);
                 assert (! send.empty());
                 assert (send.is_contained_in(extr));
-                boxes.at(rl+1).at(cc).at(ml).recv_ref_coarse.at(c ).push_back(recv);
-                boxes.at(rl  ).at(c ).at(ml).send_ref_fine  .at(cc).push_back(send);
+                boxes.at(rl+1).at(cc).at(ml).recv_ref_coarse.at(c )
+                  .push_back(recv);
+                boxes.at(rl  ).at(c ).at(ml).send_ref_fine  .at(cc)
+                  .push_back(send);
               }
       	    }
             
@@ -239,6 +242,7 @@ void dh<D>::recompose (const int initialise_from, const bool do_prolongate) {
             const ibbox& extrf = boxes.at(rl+1).at(cc).at(ml).exterior;
             const ibset& bndsf = boxes.at(rl+1).at(cc).at(ml).boundaries;
       	    // Prolongation (boundaries)
+            // TODO: prefer boxes from the same processor
       	    {
               // (the prolongation may use the exterior of the coarse
               // grid, and must fill all of the boundary of the fine
@@ -264,14 +268,14 @@ void dh<D>::recompose (const int initialise_from, const bool do_prolongate) {
               {
                 for (typename ibset::const_iterator pbi=pbndsf.begin();
                      pbi!=pbndsf.end(); ++pbi) {
-                  buffers |= (*pbi).expand(buffer_width, buffer_width) & intrf;
+                  buffers |= (*pbi).expand(buffer_width, buffer_width) & extrf;
                 }
                 buffers.normalize();
               }
               // Add boundaries
               const ibbox maxrecvs
                 = extr.expand(-pss,-pss).contracted_for(extrf);
-              ibset recvs = (bndsf | buffers) & maxrecvs;
+              ibset recvs = buffers & maxrecvs;
               recvs.normalize();
               {
                 // Do not prolongate what is already prolongated
@@ -342,8 +346,10 @@ void dh<D>::recompose (const int initialise_from, const bool do_prolongate) {
                 if (! recv.empty()) {
                   const ibbox & send = recv.expanded_for(intrf);
                   assert (! send.empty());
-                  boxes.at(rl+1).at(cc).at(ml).send_ref_coarse.at(c ).push_back(send);
-                  boxes.at(rl  ).at(c ).at(ml).recv_ref_fine  .at(cc).push_back(recv);
+                  boxes.at(rl+1).at(cc).at(ml).send_ref_coarse.at(c )
+                    .push_back(send);
+                  boxes.at(rl  ).at(c ).at(ml).recv_ref_fine  .at(cc)
+                    .push_back(recv);
                 }
               }
       	    }
diff --git a/Carpet/CarpetLib/src/gdata.cc b/Carpet/CarpetLib/src/gdata.cc
index 86537cea6..948047e5d 100644
--- a/Carpet/CarpetLib/src/gdata.cc
+++ b/Carpet/CarpetLib/src/gdata.cc
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gdata.cc,v 1.25 2003/11/21 13:55:46 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gdata.cc,v 1.26 2004/01/25 14:57:30 schnetter Exp $
 
 #include <assert.h>
 #include <stdlib.h>
@@ -54,10 +54,11 @@ comm_state<D>::~comm_state ()
 
 // Constructors
 template<int D>
-gdata<D>::gdata (const int varindex_)
-  : varindex(varindex_),
-    _has_storage(false),
-    transport_operator (find_transport_operator(varindex_))
+gdata<D>::gdata (const int varindex_, const operator_type transport_operator_)
+  : varindex(varindex_), transport_operator(transport_operator_),
+    wtime_isend(0.0), wtime_isendwait(0.0),
+    wtime_irecv(0.0), wtime_irecvwait(0.0),
+    _has_storage(false)
 { }
 
 // Destructors
@@ -66,43 +67,6 @@ gdata<D>::~gdata () { }
 
 
 
-// Transport operator types
-template<int D>
-typename gdata<D>::operator_type
-gdata<D>::find_transport_operator (const int varindex)
-{
-  const operator_type default_operator = op_Lagrange;
-  if (varindex == -1) return op_error;
-  assert (varindex >= 0);
-  const int groupindex = CCTK_GroupIndexFromVarI (varindex);
-  assert (groupindex >= 0);
-  const int group_tags_table = CCTK_GroupTagsTableI(groupindex);
-  assert (group_tags_table >= 0);
-  char prolong_string[100];
-  const int ierr = Util_TableGetString
-    (group_tags_table, sizeof prolong_string, prolong_string, "Prolongation");
-  if (ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY) {
-    char* groupname = CCTK_GroupName(groupindex);
-    CCTK_VWarn (4, __LINE__, __FILE__, CCTK_THORNSTRING,
-                "Tags table for group \"%s\" does not contain \"Prolongation\" tag", groupname);
-    ::free (groupname);
-    return default_operator;
-  }
-  assert (ierr >= 0);
-  if (CCTK_Equals(prolong_string, "None")) {
-    return op_none;
-  } else if (CCTK_Equals(prolong_string, "Lagrange")) {
-    return op_Lagrange;
-  } else if (CCTK_Equals(prolong_string, "TVD")) {
-    return op_TVD;
-  } else {
-    assert (0);
-  }
-  return op_error;
-}
-
-
-
 // Data manipulators
 template<int D>
 void gdata<D>::copy_from (comm_state<D>& state,
@@ -174,7 +138,7 @@ void gdata<D>::copy_from_recv (comm_state<D>& state,
   } else {
     
     // copy to different processor
-    gdata<D>* const tmp = make_typed(varindex);
+    gdata<D>* const tmp = make_typed(varindex, transport_operator);
     // TODO: is this efficient?
     state.tmps.push_back (tmp);
     ++state.current;
@@ -350,7 +314,7 @@ void gdata<D>
   } else {
     // interpolate from other processor
     
-    gdata<D>* const tmp = make_typed(varindex);
+    gdata<D>* const tmp = make_typed(varindex, transport_operator);
     // TODO: is this efficient?
     state.tmps.push_back (tmp);
     ++state.current;
diff --git a/Carpet/CarpetLib/src/gdata.hh b/Carpet/CarpetLib/src/gdata.hh
index 37b544d7d..7484d75da 100644
--- a/Carpet/CarpetLib/src/gdata.hh
+++ b/Carpet/CarpetLib/src/gdata.hh
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gdata.hh,v 1.21 2003/11/21 13:55:46 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gdata.hh,v 1.22 2004/01/25 14:57:30 schnetter Exp $
 
 #ifndef GDATA_HH
 #define GDATA_HH
@@ -15,6 +15,7 @@
 #include "defs.hh"
 #include "dist.hh"
 #include "bbox.hh"
+#include "operators.hh"
 #include "vect.hh"
 
 using namespace std;
@@ -56,6 +57,10 @@ protected:                      // should be readonly
 
   // Fields
   int varindex;                 // Cactus variable index, or -1
+  operator_type transport_operator;
+  
+  double wtime_isend, wtime_isendwait;
+  double wtime_irecv, wtime_irecvwait;
   
   bool _has_storage;		// has storage associated (on some processor)
   bool _owns_storage;		// owns the storage
@@ -72,13 +77,16 @@ protected:                      // should be readonly
 public:
 
   // Constructors
-  gdata (const int varindex);
+  gdata (const int varindex,
+         const operator_type transport_operator = op_error);
 
   // Destructors
   virtual ~gdata ();
 
   // Pseudo constructors
-  virtual gdata<D>* make_typed (const int varindex) const = 0;
+  virtual gdata<D>*
+  make_typed (const int varindex,
+              const operator_type transport_operator = op_error) const = 0;
   
   // Processor management
   virtual void change_processor (comm_state<D>& state,
@@ -144,14 +152,6 @@ public:
     return dot(ind, stride());
   }
   
-  // Transport operator types
- protected:
-  enum operator_type { op_error, op_none, op_Lagrange, op_TVD };
-  // readonly
-  operator_type transport_operator;
- private:
-  static operator_type find_transport_operator (const int varindex);
-  
   // Data manipulators
  public:
     void copy_from (comm_state<D>& state, const gdata* src, const ibbox& box);
diff --git a/Carpet/CarpetLib/src/gf.cc b/Carpet/CarpetLib/src/gf.cc
index 4a913d199..2fcec0d2f 100644
--- a/Carpet/CarpetLib/src/gf.cc
+++ b/Carpet/CarpetLib/src/gf.cc
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gf.cc,v 1.14 2003/11/05 16:18:39 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gf.cc,v 1.15 2004/01/25 14:57:30 schnetter Exp $
 
 #include <assert.h>
 
@@ -15,9 +15,11 @@ using namespace std;
 // Constructors
 // VGF
 template<class T,int D>
-gf<T,D>::gf (const int varindex, th<D>& t, dh<D>& d,
+gf<T,D>::gf (const int varindex, const operator_type transport_operator,
+             th<D>& t, dh<D>& d,
 	     const int tmin, const int tmax, const int prolongation_order_time)
-  : ggf<D>(varindex, t, d, tmin, tmax, prolongation_order_time)
+  : ggf<D>(varindex, transport_operator,
+           t, d, tmin, tmax, prolongation_order_time)
 {
   // VGF
   this->recompose (0, true);
diff --git a/Carpet/CarpetLib/src/gf.hh b/Carpet/CarpetLib/src/gf.hh
index ee4ce684c..3f2027485 100644
--- a/Carpet/CarpetLib/src/gf.hh
+++ b/Carpet/CarpetLib/src/gf.hh
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gf.hh,v 1.10 2003/11/21 13:55:46 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gf.hh,v 1.11 2004/01/25 14:57:30 schnetter Exp $
 
 #ifndef GF_HH
 #define GF_HH
@@ -43,7 +43,8 @@ public:
   
   // Constructors
   // VGF
-  gf (const int varindex, th<D>& t, dh<D>& d,
+  gf (const int varindex, const operator_type transport_operator,
+      th<D>& t, dh<D>& d,
       const int tmin, const int tmax, const int prolongation_order_time);
   
   // Destructors
@@ -55,7 +56,8 @@ public:
   
 protected:
   
-  virtual gdata<D>* typed_data() { return new data<T,D>(this->varindex); }
+  virtual gdata<D>* typed_data()
+  { return new data<T,D>(this->varindex, this->transport_operator); }
   
   
   
diff --git a/Carpet/CarpetLib/src/ggf.cc b/Carpet/CarpetLib/src/ggf.cc
index 6786e5eb7..88247a6bf 100644
--- a/Carpet/CarpetLib/src/ggf.cc
+++ b/Carpet/CarpetLib/src/ggf.cc
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/ggf.cc,v 1.29 2003/11/05 16:18:39 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/ggf.cc,v 1.30 2004/01/25 14:57:30 schnetter Exp $
 
 #include <assert.h>
 #include <stdlib.h>
@@ -20,16 +20,17 @@ using namespace std;
 
 // Constructors
 template<int D>
-ggf<D>::ggf (const int varindex, th<D>& t, dh<D>& d,
+ggf<D>::ggf (const int varindex, const operator_type transport_operator,
+             th<D>& t, dh<D>& d,
              const int tmin, const int tmax,
              const int prolongation_order_time)
-  : varindex(varindex), t(t),
+  : varindex(varindex), transport_operator(transport_operator), t(t),
     tmin(tmin), tmax(tmax),
     prolongation_order_time(prolongation_order_time),
     h(d.h), d(d),
     storage(tmax-tmin+1)
 {
-  assert (t.h == &d.h);
+  assert (&t.h == &d.h);
 
   d.add(this);
 }
@@ -140,13 +141,16 @@ void ggf<D>::recompose (const int initialise_from, const bool do_prolongate) {
           for (comm_state<D> state; !state.done(); state.step()) {
             sync (state,tl,rl,c,ml);
           }
-          // TODO: assert that reflevel 0 boundaries are copied
-          if (rl>0) {
-            for (comm_state<D> state; !state.done(); state.step()) {
-              const CCTK_REAL time = t.time(tl,rl,ml);
-              ref_bnd_prolongate (state,tl,rl,c,ml,time);
-            }
-          } // if rl
+          
+          if (do_prolongate) {
+            // TODO: assert that reflevel 0 boundaries are copied
+            if (rl>0) {
+              for (comm_state<D> state; !state.done(); state.step()) {
+                const CCTK_REAL time = t.time(tl,rl,ml);
+                ref_bnd_prolongate (state,tl,rl,c,ml,time);
+              }
+            } // if rl
+          }
           
       	} // for ml
       } // for c
@@ -318,7 +322,6 @@ void ggf<D>::intercat (comm_state<D>& state,
     gsrcs[i] = storage[tl2s[i]-tmin][rl2][c2][ml2];
     times[i] = t.time(tl2s[i],rl2,ml2);
   }
-//   const CCTK_REAL time = t.time(tl1,rl1,ml1);
   
   const ibbox recv = d.boxes[rl1][c1][ml1].*recv_list;
   const ibbox send = d.boxes[rl2][c2][ml2].*send_list;
@@ -359,7 +362,6 @@ void ggf<D>::intercat (comm_state<D>& state,
     gsrcs[i] = storage[tl2s[i]-tmin][rl2][c2][ml2];
     times[i] = t.time(tl2s[i],rl2,ml2);
   }
-//   const CCTK_REAL time = t.time(tl1,rl1,ml1);
   
   const iblist recv = d.boxes[rl1][c1][ml1].*recv_list;
   const iblist send = d.boxes[rl2][c2][ml2].*send_list;
@@ -405,7 +407,6 @@ void ggf<D>::intercat (comm_state<D>& state,
       gsrcs[i] = storage[tl2s[i]-tmin][rl2][c2][ml2];
       times[i] = t.time(tl2s[i],rl2,ml2);
     }
-//     const CCTK_REAL time = t.time(tl1,rl1,ml1);
     
     const iblist recv = (d.boxes[rl1][c1][ml1].*recv_listvect)[c2];
     const iblist send = (d.boxes[rl2][c2][ml2].*send_listvect)[c1];
diff --git a/Carpet/CarpetLib/src/ggf.hh b/Carpet/CarpetLib/src/ggf.hh
index b99c7aa10..59f9bfd71 100644
--- a/Carpet/CarpetLib/src/ggf.hh
+++ b/Carpet/CarpetLib/src/ggf.hh
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/ggf.hh,v 1.17 2003/11/05 16:18:39 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/ggf.hh,v 1.18 2004/01/25 14:57:30 schnetter Exp $
 
 #ifndef GGF_HH
 #define GGF_HH
@@ -52,6 +52,7 @@ public:				// should be readonly
 
   // Fields
   int varindex;                 // Cactus variable index
+  operator_type transport_operator;
   
   th<D> &t;			// time hierarchy
   int tmin, tmax;		// timelevels
@@ -66,7 +67,8 @@ protected:
 public:
 
   // Constructors
-  ggf (const int varindex, th<D>& t, dh<D>& d,
+  ggf (const int varindex, const operator_type transport_operator,
+       th<D>& t, dh<D>& d,
        const int tmin, const int tmax,
        const int prolongation_order_time);
 
diff --git a/Carpet/CarpetLib/src/gh.cc b/Carpet/CarpetLib/src/gh.cc
index b37fb3f15..fdaa3d054 100644
--- a/Carpet/CarpetLib/src/gh.cc
+++ b/Carpet/CarpetLib/src/gh.cc
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gh.cc,v 1.23 2003/11/05 16:18:39 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gh.cc,v 1.24 2004/01/25 14:57:30 schnetter Exp $
 
 #include <assert.h>
 #include <stdlib.h>
@@ -21,7 +21,7 @@ using namespace std;
 template<int D>
 gh<D>::gh (const int reffact, const centering refcent,
 	   const int mgfact, const centering mgcent,
-	   const ibbox& baseextent)
+           const ibbox baseextent)
   : reffact(reffact), refcent(refcent),
     mgfact(mgfact), mgcent(mgcent),
     baseextent(baseextent)
@@ -169,69 +169,6 @@ void gh<D>::recompose (const rexts& exts,
   }
 }
 
-// Helpers
-template<int D>
-typename gh<D>::cexts
-gh<D>::make_reflevel_multigrid_boxes (const vector<ibbox>& exts,
-				      const int mglevels)
-  const
-{
-  assert (mglevels>0);
-  
-  cexts mexts (exts.size());
-  for (int c=0; c<(int)exts.size(); ++c) {
-    
-    mexts[c].resize(mglevels);
-    
-    ibbox ext = exts[c];
-    for (int ml=0; ml<mglevels; ++ml) {
-      
-      mexts[c][ml] = ext;
-      
-      if (ml == mglevels-1) break;
-      
-      // This level's characteristics
-      ivect str = ext.stride();
-      ivect lo = ext.lower();
-      ivect up = ext.upper();
-      
-      // Transform to next (coarser) level
-      switch (mgcent) {
-      case vertex_centered:
-	break;
-      case cell_centered:
-	for (int d=0; d<D; ++d) assert (str[d]%2 == 0);
-	lo += str/2;
-	break;
-      default:
-	assert (0);
-      }
-      str *= mgfact;
-      up = up - (up - lo) % str;
-      
-      ext = ibbox(lo,up,str);
-    } // for ml
-  } // for c
-  
-  return mexts;
-}
-
-template<int D>
-typename gh<D>::rexts
-gh<D>::make_multigrid_boxes (const vector<vector<ibbox> >& exts,
-			     const int mglevels)
-  const
-{
-  assert (mglevels>0);
-  
-  rexts mexts (exts.size());
-  for (int rl=0; rl<(int)exts.size(); ++rl) {
-    mexts[rl] = make_reflevel_multigrid_boxes (exts[rl], mglevels);
-  }
-  
-  return mexts;
-}
-
 
 
 // Accessors
@@ -277,7 +214,6 @@ ostream& gh<D>::output (ostream& os) const {
   os << "gh<" << D << ">:"
      << "reffactor=" << reffact << ",refcentering=" << refcent << ","
      << "mgfactor=" << mgfact << ",mgcentering=" << mgcent << ","
-     << "baseextent=" << baseextent << ","
      << "extents=" << extents << ","
      << "outer_boundaries=" << outer_boundaries << ","
      << "processors=" << processors << ","
diff --git a/Carpet/CarpetLib/src/gh.hh b/Carpet/CarpetLib/src/gh.hh
index db046185a..7db23b780 100644
--- a/Carpet/CarpetLib/src/gh.hh
+++ b/Carpet/CarpetLib/src/gh.hh
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gh.hh,v 1.14 2003/11/05 16:18:39 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gh.hh,v 1.15 2004/01/25 14:57:30 schnetter Exp $
 
 #ifndef GH_HH
 #define GH_HH
@@ -63,7 +63,7 @@ public:				// should be readonly
   
   list<th<D>*> ths;		// list of all time hierarchies
   
-  ibbox baseextent;		// bounds (inclusive) of base level
+  ibbox baseextent;
   vector<vector<ibbox> > bases; // [rl][ml]
   
   // TODO: invent structure for this
@@ -78,7 +78,7 @@ public:
   // Constructors
   gh (const int reffact, const centering refcent,
       const int mgfact, const centering mgcent,
-      const ibbox& baseextent);
+      const ibbox baseextent);
   
   // Destructors
   virtual ~gh ();
@@ -89,14 +89,6 @@ public:
                   const int initialise_from,
                   const bool do_prolongate);
   
-  // Helpers
-  cexts make_reflevel_multigrid_boxes (const vector<ibbox>& exts,
-				       const int mglevels)
-    const;
-  rexts make_multigrid_boxes (const vector<vector<ibbox> >& exts,
-			      const int mglevels)
-    const;
-  
   // Accessors
   int reflevels () const {
     return (int)extents.size();
diff --git a/Carpet/CarpetLib/src/make.code.defn b/Carpet/CarpetLib/src/make.code.defn
index 0b718747a..7cbecca46 100644
--- a/Carpet/CarpetLib/src/make.code.defn
+++ b/Carpet/CarpetLib/src/make.code.defn
@@ -1,5 +1,5 @@
 # Main make.code.defn file for thorn CarpetLib   -*-Makefile-*-
-# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/make.code.defn,v 1.10 2003/11/05 16:18:39 schnetter Exp $
+# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/make.code.defn,v 1.11 2004/01/25 14:57:30 schnetter Exp $
 
 # Source files in this directory
 SRCS =	bbox.cc					\
@@ -23,7 +23,9 @@ SRCS =	bbox.cc					\
 	prolongate_3d_real8_o3_rf2.F77		\
 	prolongate_3d_real8_o5.F77		\
 	prolongate_3d_real8_2tl.F77		\
+	prolongate_3d_real8_2tl_rf2.F77		\
 	prolongate_3d_real8_2tl_o3.F77		\
+	prolongate_3d_real8_2tl_o3_rf2.F77	\
 	prolongate_3d_real8_2tl_o5.F77		\
 	prolongate_3d_real8_3tl.F77		\
 	prolongate_3d_real8_3tl_rf2.F77		\
diff --git a/Carpet/CarpetLib/src/operators.hh b/Carpet/CarpetLib/src/operators.hh
index b65ff23cb..d5ea168f3 100644
--- a/Carpet/CarpetLib/src/operators.hh
+++ b/Carpet/CarpetLib/src/operators.hh
@@ -1,8 +1,8 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/operators.hh,v 1.2 2004/03/03 15:30:40 hawke Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/operators.hh,v 1.1 2004/01/25 14:57:30 schnetter Exp $
 
 #ifndef OPERATORS_HH
 #define OPERATORS_HH
 
-enum operator_type { op_error, op_none, op_Lagrange, op_TVD, op_ENO };
+enum operator_type { op_error, op_none, op_Lagrange, op_TVD };
 
 #endif // OPERATORS_HH
diff --git a/Carpet/CarpetLib/src/prolongate_3d_real8_2tl_o3_rf2.F77 b/Carpet/CarpetLib/src/prolongate_3d_real8_2tl_o3_rf2.F77
index d9d62741b..0945909bf 100644
--- a/Carpet/CarpetLib/src/prolongate_3d_real8_2tl_o3_rf2.F77
+++ b/Carpet/CarpetLib/src/prolongate_3d_real8_2tl_o3_rf2.F77
@@ -1,5 +1,5 @@
 c     -*-Fortran-*-
-c     $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_2tl_o3_rf2.F77,v 1.2 2004/03/11 12:03:09 schnetter Exp $
+c     $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_2tl_o3_rf2.F77,v 1.1 2004/01/25 14:57:30 schnetter Exp $
       
 #include "cctk.h"
 #include "cctk_Parameters.h"
@@ -85,11 +85,6 @@ c     bbox(:,3) is stride
 c     This could be handled, but is likely to point to an error elsewhere
             call CCTK_WARN (0, "Internal error: region extent is empty")
          end if
-         if (mod(srcbbox(d,2) - srcbbox(d,1), srcbbox(d,3)).ne.0
-     $        .or. mod(dstbbox(d,2) - dstbbox(d,1), dstbbox(d,3)).ne.0
-     $        .or. mod(regbbox(d,2) - regbbox(d,1), regbbox(d,3)).ne.0) then
-            call CCTK_WARN (0, "Internal error: array extents are not integer multiples of the strides")
-         end if
          regkext = (regbbox(d,2) - regbbox(d,1)) / regbbox(d,3) + 1
          srckoff = (regbbox(d,1) - srcbbox(d,1)) / dstbbox(d,3)
          offsetlo = regbbox(d,3)
diff --git a/Carpet/CarpetLib/src/prolongate_3d_real8_2tl_rf2.F77 b/Carpet/CarpetLib/src/prolongate_3d_real8_2tl_rf2.F77
index 727c2581f..0708979f6 100644
--- a/Carpet/CarpetLib/src/prolongate_3d_real8_2tl_rf2.F77
+++ b/Carpet/CarpetLib/src/prolongate_3d_real8_2tl_rf2.F77
@@ -1,5 +1,5 @@
 c     -*-Fortran-*-
-c     $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_2tl_rf2.F77,v 1.2 2004/03/11 12:03:09 schnetter Exp $
+c     $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_2tl_rf2.F77,v 1.1 2004/01/25 14:57:30 schnetter Exp $
       
 #include "cctk.h"
 #include "cctk_Parameters.h"
@@ -75,11 +75,6 @@ c     bbox(:,3) is stride
 c     This could be handled, but is likely to point to an error elsewhere
             call CCTK_WARN (0, "Internal error: region extent is empty")
          end if
-         if (mod(srcbbox(d,2) - srcbbox(d,1), srcbbox(d,3)).ne.0
-     $        .or. mod(dstbbox(d,2) - dstbbox(d,1), dstbbox(d,3)).ne.0
-     $        .or. mod(regbbox(d,2) - regbbox(d,1), regbbox(d,3)).ne.0) then
-            call CCTK_WARN (0, "Internal error: array extents are not integer multiples of the strides")
-         end if
          if (regbbox(d,1).lt.srcbbox(d,1)
      $        .or. regbbox(d,1).lt.dstbbox(d,1)
      $        .or. regbbox(d,2).gt.srcbbox(d,2)
diff --git a/Carpet/CarpetLib/src/th.cc b/Carpet/CarpetLib/src/th.cc
index d849e3c47..5fd8eabee 100644
--- a/Carpet/CarpetLib/src/th.cc
+++ b/Carpet/CarpetLib/src/th.cc
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/th.cc,v 1.12 2003/08/10 21:58:45 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/th.cc,v 1.13 2004/01/25 14:57:30 schnetter Exp $
 
 #include <assert.h>
 #include <math.h>
@@ -18,23 +18,23 @@ using namespace std;
 
 // Constructors
 template<int D>
-th<D>::th (gh<D>* h, const CCTK_REAL basedelta)
+th<D>::th (gh<D>& h, const CCTK_REAL basedelta)
   : h(h), delta(basedelta) {
-  h->add(this);
+  h.add(this);
 }
 
 // Destructors
 template<int D>
 th<D>::~th () {
-  h->remove(this);
+  h.remove(this);
 }
 
 // Modifiers
 template<int D>
 void th<D>::recompose () {
-  times.resize(h->reflevels());
-  deltas.resize(h->reflevels());
-  for (int rl=0; rl<h->reflevels(); ++rl) {
+  times.resize(h.reflevels());
+  deltas.resize(h.reflevels());
+  for (int rl=0; rl<h.reflevels(); ++rl) {
     const int old_mglevels = times[rl].size();
     CCTK_REAL mgtime;
     // Select default time
@@ -45,15 +45,15 @@ void th<D>::recompose () {
     } else {
       mgtime = times[rl][old_mglevels-1];
     }
-    times[rl].resize(h->mglevels(rl,0), mgtime);
-    deltas[rl].resize(h->mglevels(rl,0));
-    for (int ml=0; ml<h->mglevels(rl,0); ++ml) {
+    times[rl].resize(h.mglevels(rl,0), mgtime);
+    deltas[rl].resize(h.mglevels(rl,0));
+    for (int ml=0; ml<h.mglevels(rl,0); ++ml) {
       if (rl==0 && ml==0) {
 	deltas[rl][ml] = delta;
       } else if (ml==0) {
-	deltas[rl][ml] = deltas[rl-1][ml] / h->reffact;
+	deltas[rl][ml] = deltas[rl-1][ml] / h.reffact;
       } else {
-	deltas[rl][ml] = deltas[rl][ml-1] * h->mgfact;
+	deltas[rl][ml] = deltas[rl][ml-1] * h.mgfact;
       }
     }
   }
@@ -66,8 +66,8 @@ template<int D>
 void th<D>::output (ostream& os) const {
   os << "th<" << D << ">:"
      << "times={";
-  for (int rl=0; rl<h->reflevels(); ++rl) {
-    for (int ml=0; ml<h->mglevels(rl,0); ++ml) {
+  for (int rl=0; rl<h.reflevels(); ++rl) {
+    for (int ml=0; ml<h.mglevels(rl,0); ++ml) {
       if (!(rl==0 && ml==0)) os << ",";
       os << rl << ":" << ml << ":"
 	 << times[rl][ml] << "(" << deltas[rl][ml] << ")";
diff --git a/Carpet/CarpetLib/src/th.hh b/Carpet/CarpetLib/src/th.hh
index 85d028ac3..26deb154a 100644
--- a/Carpet/CarpetLib/src/th.hh
+++ b/Carpet/CarpetLib/src/th.hh
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/th.hh,v 1.9 2003/01/03 15:49:36 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/th.hh,v 1.10 2004/01/25 14:57:30 schnetter Exp $
 
 #ifndef TH_HH
 #define TH_HH
@@ -33,7 +33,7 @@ class th {
 public:				// should be readonly
   
   // Fields
-  gh<D> *h;                     // hierarchy
+  gh<D>& h;                     // hierarchy
   
 private:
   
@@ -44,7 +44,7 @@ private:
 public:
   
   // Constructors
-  th (gh<D>* h, const CCTK_REAL basedelta);
+  th (gh<D>& h, const CCTK_REAL basedelta);
   
   // Destructors
   ~th ();
@@ -54,14 +54,14 @@ public:
   
   // Time management
   CCTK_REAL get_time (const int rl, const int ml) const {
-    assert (rl>=0 && rl<h->reflevels());
-    assert (ml>=0 && ml<h->mglevels(rl,0));
+    assert (rl>=0 && rl<h.reflevels());
+    assert (ml>=0 && ml<h.mglevels(rl,0));
     return times[rl][ml];
   }
   
   void set_time (const int rl, const int ml, const CCTK_REAL t) {
-    assert (rl>=0 && rl<h->reflevels());
-    assert (ml>=0 && ml<h->mglevels(rl,0));
+    assert (rl>=0 && rl<h.reflevels());
+    assert (ml>=0 && ml<h.mglevels(rl,0));
     times[rl][ml] = t;
   }
   
@@ -70,20 +70,20 @@ public:
   }
   
   CCTK_REAL get_delta (const int rl, const int ml) const {
-    assert (rl>=0 && rl<h->reflevels());
-    assert (ml>=0 && ml<h->mglevels(rl,0));
+    assert (rl>=0 && rl<h.reflevels());
+    assert (ml>=0 && ml<h.mglevels(rl,0));
     return deltas[rl][ml];
   }
   
   void set_delta (const int rl, const int ml, const CCTK_REAL dt) {
-    assert (rl>=0 && rl<h->reflevels());
-    assert (ml>=0 && ml<h->mglevels(rl,0));
+    assert (rl>=0 && rl<h.reflevels());
+    assert (ml>=0 && ml<h.mglevels(rl,0));
     deltas[rl][ml] = dt;
   }
   
   CCTK_REAL time (const int tl, const int rl, const int ml) const {
-    assert (rl>=0 && rl<h->reflevels());
-    assert (ml>=0 && ml<h->mglevels(rl,0));
+    assert (rl>=0 && rl<h.reflevels());
+    assert (ml>=0 && ml<h.mglevels(rl,0));
     return get_time(rl, ml) + tl * get_delta(rl, ml);
   }
   
diff --git a/Carpet/CarpetLib/src/vect.cc b/Carpet/CarpetLib/src/vect.cc
index 7fe311fdb..b5773eebf 100644
--- a/Carpet/CarpetLib/src/vect.cc
+++ b/Carpet/CarpetLib/src/vect.cc
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/vect.cc,v 1.12 2003/11/13 16:03:58 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/vect.cc,v 1.13 2004/01/25 14:57:30 schnetter Exp $
 
 #include <assert.h>
 
@@ -16,19 +16,16 @@ using namespace std;
 template<class T,int D>
 void vect<T,D>::input (istream& is) {
   skipws (is);
-  assert (is.peek() == '[');
-  is.get();
+  consume (is, '[');
   for (int d=0; d<D; ++d) {
     is >> (*this)[d];
     if (d<D-1) {
       skipws (is);
-      assert (is.peek() == ',');
-      is.get();
+      consume (is, ',');
     }
   }
   skipws (is);
-  assert (is.peek() == ']');
-  is.get();
+  consume (is, ']');
 }
 
 
diff --git a/Carpet/CarpetLib/src/vect.hh b/Carpet/CarpetLib/src/vect.hh
index 5271626d2..1df8f874c 100644
--- a/Carpet/CarpetLib/src/vect.hh
+++ b/Carpet/CarpetLib/src/vect.hh
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/vect.hh,v 1.21 2003/11/13 16:03:58 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/vect.hh,v 1.22 2004/01/25 14:57:30 schnetter Exp $
 
 #ifndef VECT_HH
 #define VECT_HH
@@ -82,14 +82,17 @@ public:
     elt[0]=x; elt[1]=y; elt[2]=z; elt[3]=t;
   }
   
+#if 0
+  // This creates confusion
   /** Constructor from a pointer, i.e.\ a C array.  */
-  vect (const T* const x) {
+  explicit vect (const T* const x) {
     for (int d=0; d<D; ++d) elt[d]=x[d];
   }
+#endif
   
   /** Constructor from a vector with a different type.  */
   template<class S>
-  explicit vect (const vect<S,D>& a) {
+  /*explicit*/ vect (const vect<S,D>& a) {
     for (int d=0; d<D; ++d) elt[d]=(T)a[d];
   }
   
@@ -160,8 +163,9 @@ public:
   // Accessors
   
   /** Return a non-writable element of a vector.  */
-  // Don't return a reference; *this might be a temporary
-  T operator[] (const int d) const {
+  // (Don't return a reference; *this might be a temporary)
+  // Do return a reference, so that a vector can be accessed as array
+  const T& operator[] (const int d) const {
     assert(d>=0 && d<D);
     return elt[d];
   }
@@ -172,6 +176,14 @@ public:
     return elt[d];
   }
   
+#if 0
+  // This creates confusion
+  /** Return a pointer to a vector.  */
+  operator const T* () const {
+    return this;
+  }
+#endif
+  
   /** Return a combination of the vector elements e[a[i]].  The
       element combination is selected by another vector.  */
   template<class TT, int DD>
@@ -492,6 +504,25 @@ public:
 
 // Operators
 
+/** This corresponds to the ?: operator.  Return a vector with the
+    elements set to either b[i] or c[i], depending on whether a[i] is
+    true or not.  */
+template<class T,int D>
+inline vect<T,D> either (const vect<bool,D>& a,
+                         const vect<T,D>& b, const vect<T,D>& c) {
+  vect<T,D> r;
+  for (int d=0; d<D; ++d) r[d]=a[d]?b[d]:c[d];
+  return r;
+}
+
+/** Transpose a vector of a vector */
+template<class T, int D, int DD>
+inline vect<vect<T,D>,DD> xpose (vect<vect<T,DD>,D> const & a) {
+  vect<vect<T,D>,DD> r;
+  for (int dd=0; dd<DD; ++dd) for (int d=0; d<D; ++d) r[dd][d] = a[d][dd];
+  return r;
+}
+
 /** Return the element-wise absolute value.  */
 template<class T,int D>
 inline vect<T,D> abs (const vect<T,D>& a) {