Carpet now passes the Cart3d test suite. This required a rewrite of

Carpet now passes the Cart3d test suite.  This required a rewrite of
CarpetSlabe, which now uses the data<> class instead of doing things
its own way.  This in turn required some updates to data<>.

CarpetSlab now ignores (i. e. does not use) boundaries, except outer
boundaries.  For that to work, the dh<> class has been given the
notion of an outer boundary.

In order to increase performance, the cctk_bbox[] information in
Carpet is not set correctly.  As an approximation, it just reflects
whether the current grid component lies at the outer boundary, as
calculated using cctk_lbnd, cctk_ubnd, and cctk_gsh.

darcs-hash:20010310195459-f6438-8ee874d722de4df717e5eb3002722ceb8d4a9bb7.gz

15 files changed, 475 insertions, 527 deletions

diff --git a/Carpet/Carpet/param.ccl b/Carpet/Carpet/param.ccl
index 95ce595f1..13eb6b396 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.2 2001/03/07 12:59:40 eschnett Exp $
+# $Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/param.ccl,v 1.3 2001/03/10 20:54:59 eschnett Exp $
 
 shares: Cactus
 
@@ -27,6 +27,11 @@ CCTK_INT global_nz "Grid size in z direction"
   0:* :: "must be nonnegative"
 } 10
 
+CCTK_INT global_nsize "Grid size in each spatial direction"
+{
+  -1:* :: "must be nonnegative, or -1 to use the individual parameters"
+} -1
+
 
 
 CCTK_INT ghost_size_x "Ghost zones in x direction"
@@ -44,6 +49,11 @@ CCTK_INT ghost_size_z "Ghost zones in z direction"
   0:* :: "must be nonnegative"
 } 1
 
+CCTK_INT ghost_size "Ghost zones in each spatial direction"
+{
+  -1:* :: "must be nonnegative, or -1 to use the individual parameters"
+} -1
+
 
 
 CCTK_INT max_refinement_levels "Maximum number of refinement levels (including the base level)"
diff --git a/Carpet/Carpet/src/carpet.cc b/Carpet/Carpet/src/carpet.cc
index 7edb60644..9e1ae2cfb 100644
--- a/Carpet/Carpet/src/carpet.cc
+++ b/Carpet/Carpet/src/carpet.cc
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Attic/carpet.cc,v 1.4 2001/03/07 12:59:53 eschnett Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Attic/carpet.cc,v 1.5 2001/03/10 20:55:03 eschnett 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
@@ -32,6 +32,8 @@
 
 #include "carpet.hh"
 
+static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Attic/carpet.cc,v 1.5 2001/03/10 20:55:03 eschnett Exp $";
+
 
 
 namespace Carpet {
@@ -67,6 +69,9 @@ namespace Carpet {
   // data for grid functions
   vector<gfdesc> gfdata;	// [group]
   
+  // active time level
+  int activetimelevel;		// 0 for current, 1 for next
+  
   // current position on the grid hierarchy
   int mglevel;
   int reflevel;
@@ -114,13 +119,24 @@ namespace Carpet {
     Checkpoint ("starting SetupGH...");
     
     // ghost zones
-    const vect<int,dim> lghosts(ghost_size_x, ghost_size_y, ghost_size_z);
-    const vect<int,dim> ughosts(ghost_size_x, ghost_size_y, ghost_size_z);
+    vect<int,dim> lghosts, ughosts;
+    if (ghost_size == -1) {
+      lghosts = vect<int,dim>(ghost_size_x, ghost_size_y, ghost_size_z);
+      ughosts = vect<int,dim>(ghost_size_x, ghost_size_y, ghost_size_z);
+    } else {
+      lghosts = vect<int,dim>(ghost_size, ghost_size, ghost_size);
+      ughosts = vect<int,dim>(ghost_size, ghost_size, ghost_size);
+    }
     
     // grid size
     const int stride
       = (int)floor(pow(refinement_factor, max_refinement_levels) + 0.5);
-    const vect<int,dim> npoints(global_nx, global_ny, global_nz);
+    vect<int,dim> npoints;
+    if (global_nsize == -1) {
+      npoints =  vect<int,dim>(global_nx, global_ny, global_nz);
+    } else {
+      npoints =  vect<int,dim>(global_nsize, global_nsize, global_nsize);
+    }
     
     const vect<int,dim> str(stride);
     const vect<int,dim> lb(lghosts * str);
@@ -207,6 +223,9 @@ namespace Carpet {
       }
     }
     
+    // active time level
+    activetimelevel = 0;
+    
     // current position
     reflevel  = 0;
     mglevel   = 0;
@@ -342,6 +361,10 @@ namespace Carpet {
     Checkpoint ("%*sstarting RecursiveEvolve on level %d...",
 		2*reflevel, "", reflevel);
     
+    // Move activity to next time level
+    assert (activetimelevel == 0);
+    activetimelevel = 1;
+    
     // Prestep
     CCTK_ScheduleTraverse ("CCTK_PRESTEP", cgh, CallFunction);
     
@@ -359,11 +382,15 @@ namespace Carpet {
       cgh->cctk_time += cgh->cctk_delta_time;
     }
     
+    // Restrict
+    Restrict (cgh);
+    
     // Cycle time levels
     CycleTimeLevels (cgh);
     
-    // Restrict
-    Restrict (cgh);
+    // Move activity back to current time level
+    assert (activetimelevel == 1);
+    activetimelevel = 0;
     
     // Poststep
     CCTK_ScheduleTraverse ("CCTK_POSTSTEP", cgh, CallFunction);
@@ -391,8 +418,8 @@ namespace Carpet {
     
     RecursiveShutdown (cgh);
     
-    fprintf (stdout, "--------------------------------------------------------------------------------\n"
-	     "Done.\n");
+    CCTK_PRINTSEPARATOR;
+    printf ("Done.\n");
     
     // earlier checkpoint before calling finalising MPI
     Checkpoint ("done with Shutdown.");
@@ -575,11 +602,10 @@ namespace Carpet {
 	    cgh->cctk_gsh[d]
 	      = ((base.shape() / base.stride() + dd->lghosts + dd->ughosts)[d]
 		 * cgh->cctk_levfac[d]);
-	    // TODO: set boundary information correctly
-	    cgh->cctk_bbox[2*d  ] = 1;
-	    cgh->cctk_bbox[2*d+1] = 1;
 	    cgh->cctk_lbnd[d]     = (ext.lower() / ext.stride())[d];
 	    cgh->cctk_ubnd[d]     = (ext.upper() / ext.stride())[d];
+	    cgh->cctk_bbox[2*d  ] = cgh->cctk_lbnd[d] == 0;
+	    cgh->cctk_bbox[2*d+1] = cgh->cctk_ubnd[d] == cgh->cctk_gsh[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];
@@ -750,7 +776,8 @@ namespace Carpet {
     }
     
     const int n0 = CCTK_FirstVarIndexI(group);
-    const int tl = max(0, CCTK_NumTimeLevelsFromVarI(n0) - 1);
+    const int tl = max(0, (CCTK_NumTimeLevelsFromVarI(n0)
+			   + activetimelevel - 2));
     
     switch (CCTK_GroupTypeI(group)) {
       
@@ -1001,6 +1028,7 @@ namespace Carpet {
     const int reflevels = max_refinement_levels; // arbitrary value
     const int mglevels  = 1;	// arbitrary value
     vector<vector<bbox<int,dim> > > 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;
@@ -1018,8 +1046,9 @@ namespace Carpet {
 	vect<int,dim> cstr = rstr;
 	vect<int,dim> clb = rlb;
 	vect<int,dim> cub = rub;
-	const int zstep = ((rub[dim-1] - rlb[dim-1] + rstr[dim-1])
-			   / cstr[dim-1] + nprocs) / nprocs * cstr[dim-1];
+	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];
 	clb[dim-1] = rlb[dim-1] + zstep *  c;
 	cub[dim-1] = rlb[dim-1] + zstep * (c+1);
 	if (c == nprocs-1) cub[dim-1] = rub[dim-1];
@@ -1103,7 +1132,8 @@ namespace Carpet {
       for (int group=0; group<CCTK_NumGroups(); ++group) {
 	
 	const int n0 = CCTK_FirstVarIndexI(group);
-	const int tl = max(0, CCTK_NumTimeLevelsFromVarI(n0) - 1);
+	const int tl = max(0, (CCTK_NumTimeLevelsFromVarI(n0)
+			       + activetimelevel - 2));
 	
 	switch (CCTK_GroupTypeI(group)) {
 	  
diff --git a/Carpet/Carpet/src/carpet.hh b/Carpet/Carpet/src/carpet.hh
index d94acdea4..a3cadb4c5 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.2 2001/03/05 14:30:03 eschnett Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/carpet.hh,v 1.3 2001/03/10 20:55:03 eschnett Exp $
 
 #include <vector>
 
@@ -47,6 +47,9 @@ namespace Carpet {
   };
   extern vector<gfdesc> gfdata;	// [group]
   
+  // active time level
+  extern int activetimelevel;	// 0 for current, 1 for next
+  
   // current position on the grid hierarchy
   extern int mglevel;
   extern int reflevel;
@@ -76,7 +79,6 @@ namespace Carpet {
   int EnableGroupComm (cGH *cgh, const char *groupname);
   int DisableGroupComm (cGH *cgh, const char *groupname);
   int Barrier (cGH *cgh);
-//   int ParallelInit (cGH *cgh);
   int Exit (cGH *cgh, int retval);
   int Abort (cGH *cgh, int retval);
   int myProc (cGH *cgh);
diff --git a/Carpet/CarpetLib/src/bbox.hh b/Carpet/CarpetLib/src/bbox.hh
index 6787b9240..95afb992e 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.3 2001/03/07 13:00:57 eschnett Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/bbox.hh,v 1.4 2001/03/10 20:55:06 eschnett Exp $
 
  ***************************************************************************/
 
@@ -96,6 +96,8 @@ public:
   bbox contracted_for (const bbox& b) const;
   
   // Set operations
+  // TODO: rename these; they clash with the bboxset operations
+  // (and name them & and | instead)
   // Smallest bbox containing both boxes
   bbox operator* (const bbox& b) const;
   bbox& operator*= (const bbox& b);
diff --git a/Carpet/CarpetLib/src/bboxset.cc b/Carpet/CarpetLib/src/bboxset.cc
index 1dd685178..b76dc67b6 100644
--- a/Carpet/CarpetLib/src/bboxset.cc
+++ b/Carpet/CarpetLib/src/bboxset.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/bboxset.cc,v 1.1 2001/03/01 13:40:10 eschnett Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/bboxset.cc,v 1.2 2001/03/10 20:55:06 eschnett Exp $
 
  ***************************************************************************/
 
@@ -38,7 +38,7 @@ bboxset<T,D>::bboxset () {
 
 template<class T, int D>
 bboxset<T,D>::bboxset (const box& b) {
-  bs.insert(b);
+  if (!b.empty()) bs.insert(b);
   assert (invariant());
 }
 
diff --git a/Carpet/CarpetLib/src/bboxset.hh b/Carpet/CarpetLib/src/bboxset.hh
index 009d9afdd..98aa3fb65 100644
--- a/Carpet/CarpetLib/src/bboxset.hh
+++ b/Carpet/CarpetLib/src/bboxset.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/bboxset.hh,v 1.1 2001/03/01 13:40:10 eschnett Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/bboxset.hh,v 1.2 2001/03/10 20:55:06 eschnett Exp $
 
  ***************************************************************************/
 
@@ -34,6 +34,8 @@
 // Forward definition
 template<class T, int D> class bboxset;
 
+// TODO: add more functions for the other operators
+// (but cf. bbox<T,D> first!)
 template<class T,int D>
 bboxset<T,D> operator- (const bbox<T,D>& b1, const bbox<T,D>& b2);
 template<class T,int D>
diff --git a/Carpet/CarpetLib/src/data.cc b/Carpet/CarpetLib/src/data.cc
index 3f8222f6c..b2c190865 100644
--- a/Carpet/CarpetLib/src/data.cc
+++ b/Carpet/CarpetLib/src/data.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/data.cc,v 1.3 2001/03/05 21:48:38 eschnett Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/data.cc,v 1.4 2001/03/10 20:55:06 eschnett Exp $
 
  ***************************************************************************/
 
@@ -56,13 +56,14 @@ data<T,D>::~data () {
   
 // Pseudo constructors
 template<class T, int D>
-data<T,D>* data<T,D>::make_typed (const ibbox& extent, const int proc) const {
-  return new data(extent, proc);
+data<T,D>* data<T,D>::make_typed () const {
+  return new data();
 }
 
 // Storage management
 template<class T, int D>
-void data<T,D>::allocate (const ibbox& extent, const int proc) {
+void data<T,D>::allocate (const ibbox& extent, const int proc,
+			  void* const mem=0) {
   assert (!_has_storage);
   _has_storage = true;
   // data
@@ -77,13 +78,20 @@ void data<T,D>::allocate (const ibbox& extent, const int proc) {
   int rank;
   MPI_Comm_rank (dist::comm, &rank);
   if (rank==_proc) {
-    _storage = new T[_size];
+    _owns_storage = !mem;
+    if (!mem) {
+      _storage = new T[_size];
+    } else {
+      _storage = (T*)mem;
+    }
+  } else {
+    assert (!mem);
   }
 }
 
 template<class T, int D>
 void data<T,D>::free () {
-  if (_storage) delete [] _storage;
+  if (_storage && _owns_storage) delete [] _storage;
   _storage = 0;
   _has_storage = false;
 }
@@ -99,9 +107,12 @@ void data<T,D>::transfer_from (generic_data<D>* gsrc) {
 // Processor management
 
 template<class T, int D>
-void data<T,D>::change_processor (const int newproc) {
-  if (newproc == _proc) return;
-
+void data<T,D>::change_processor (const int newproc, void* const mem=0) {
+  if (newproc == _proc) {
+    assert (!mem);
+    return;
+  }
+  
   if (_has_storage) {
     int rank;
     MPI_Comm_rank (dist::comm, &rank);
@@ -109,7 +120,12 @@ void data<T,D>::change_processor (const int newproc) {
       // copy from other processor
       
       assert (!_storage);
-      _storage = new T[_size];
+      _owns_storage = !mem;
+      if (!mem) {
+	_storage = new T[_size];
+      } else {
+	_storage = (T*)mem;
+      }
 
       T dummy;
       MPI_Status status;
@@ -118,16 +134,20 @@ void data<T,D>::change_processor (const int newproc) {
 
     } else if (rank == _proc) {
       // copy to other processor
-
+      
+      assert (!mem);
       assert (_storage);
       T dummy;
       MPI_Send (_storage, _size, dist::datatype(dummy), newproc,
 		dist::tag, dist::comm);
 
-      delete [] _storage;
+      if (_owns_storage) {
+	delete [] _storage;
+      }
       _storage = 0;
-  	
+      
     } else {
+      assert (!mem);
       assert (!_storage);
     }
   }
@@ -166,7 +186,7 @@ void data<T,D>::copy_from (const generic_data<D>* gsrc, const ibbox& box) {
   } else {
     
     // copy to different processor
-    data* tmp = new data(box, src->_proc);
+    data* const tmp = new data(box, src->_proc);
     tmp->copy_from (src, box);
     tmp->change_processor (_proc);
     copy_from (tmp, box);
@@ -228,7 +248,7 @@ void data<T,D>::interpolate_from (const generic_data<D>* gsrc,
   } else {
     // interpolate from other processor
     
-    data* tmp = new data(box, src->_proc);
+    data* const tmp = new data(box, src->_proc);
     tmp->interpolate_from (src, box);
     tmp->change_processor (_proc);
     copy_from (tmp, box);
@@ -300,9 +320,9 @@ void data<T,D>::interpolate_from (const generic_data<D>* gsrc,
   } else {
     // interpolate from other processors
 
-    data* smp = new data(box, src->_proc);
+    data* const smp = new data(box, src->_proc);
     smp->copy_from (src, box);
-    data* tmp = new data(box, trc->_proc);
+    data* const tmp = new data(box, trc->_proc);
     tmp->copy_from (trc, box);
     interpolate_from (smp, sfact, tmp, tfact, box);
     delete smp;
diff --git a/Carpet/CarpetLib/src/data.hh b/Carpet/CarpetLib/src/data.hh
index 7b4336b97..84e4c92b7 100644
--- a/Carpet/CarpetLib/src/data.hh
+++ b/Carpet/CarpetLib/src/data.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/data.hh,v 1.2 2001/03/05 14:31:03 eschnett Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/data.hh,v 1.3 2001/03/10 20:55:06 eschnett Exp $
 
  ***************************************************************************/
 
@@ -27,7 +27,6 @@
 #include "defs.hh"
 #include "dist.hh"
 #include "bbox.hh"
-#include "bboxset.hh"
 #include "gdata.hh"
 #include "vect.hh"
 
@@ -49,7 +48,6 @@ class data: public generic_data<D> {
   // Types
   typedef vect<int,D> ivect;
   typedef bbox<int,D> ibbox;
-  typedef bboxset<int,D> ibset;
 
   // Fields
   T* restrict _storage;		// the data (if located on this processor)
@@ -64,15 +62,16 @@ public:
   virtual ~data ();
 
   // Pseudo constructors
-  virtual data* make_typed (const ibbox& extent, const int proc) const;
+  virtual data* make_typed () const;
 
   // Storage management
-  virtual void allocate (const ibbox& extent, const int proc);
+  virtual void allocate (const ibbox& extent, const int proc,
+			 void* const mem=0);
   virtual void free ();
   virtual void transfer_from (generic_data<D>* gsrc);
 
   // Processor management
-  virtual void change_processor (const int newproc);
+  virtual void change_processor (const int newproc, void* const mem=0);
 
   // Accessors
   virtual const T* storage () const {
diff --git a/Carpet/CarpetLib/src/dh.cc b/Carpet/CarpetLib/src/dh.cc
index 5dc7ba0e7..756f37cad 100644
--- a/Carpet/CarpetLib/src/dh.cc
+++ b/Carpet/CarpetLib/src/dh.cc
@@ -6,7 +6,7 @@
     copyright            : (C) 2000 by Erik Schnetter
     email                : schnetter@astro.psu.edu
 
-    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/dh.cc,v 1.3 2001/03/07 13:00:57 eschnett Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/dh.cc,v 1.4 2001/03/10 20:55:06 eschnett Exp $
 
  ***************************************************************************/
 
@@ -66,18 +66,18 @@ void dh<D>::recompose () {
       boxes[rl][c].resize(h.mglevels(rl,c));
       for (int ml=0; ml<h.mglevels(rl,c); ++ml) {
        	const ibbox intr = h.extents[rl][c][ml];
-        	
+	
        	// Interior
        	// (the interior of the grid has the extent as specified by
        	// the user)
        	boxes[rl][c][ml].interior = intr;
-        	
+	
        	// Exterior (add ghost zones)
        	// (the content of the exterior is completely determined by
        	// the interior of this or other components; the content of
        	// the exterior is redundant)
        	boxes[rl][c][ml].exterior = intr.expand(lghosts, ughosts);
-        	
+	
        	// Boundaries (ghost zones only)
        	// (interior + boundaries = exterior)
        	boxes[rl][c][ml].boundaries = boxes[rl][c][ml].exterior - intr;
@@ -118,8 +118,9 @@ void dh<D>::recompose () {
       for (int ml=0; ml<h.mglevels(rl,c); ++ml) {
       	const ibbox intr = boxes[rl][c][ml].interior;
       	const ibbox extr = boxes[rl][c][ml].exterior;
+
       	const ibset bnds = boxes[rl][c][ml].boundaries;
-      	
+	
       	// Sync boxes
       	for (int cc=0; cc<h.components(rl); ++cc) {
       	  assert (ml<h.mglevels(rl,cc));
@@ -208,6 +209,31 @@ void dh<D>::recompose () {
     } // for c
   } // for rl
   
+  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) {
+	
+	// Boundaries that are neither synced nor prolonged from
+	// coarser grids (outer boundaries)
+	ibset& sync_not = boxes[rl][c][ml].sync_not;
+	
+	// The whole boundary
+	sync_not = boxes[rl][c][ml].boundaries;
+	
+	// Subtract boxes received during synchronisation
+	const iblistvect& recv_sync = boxes[rl][c][ml].recv_sync;
+	for (iblistvect::const_iterator lvi=recv_sync.begin();
+	     lvi!=recv_sync.end(); ++lvi) {
+	  for (iblist::const_iterator li=lvi->begin();
+	       li!=lvi->end(); ++li) {
+	    sync_not -= *li;
+	  }
+	}
+	
+      } // for ml
+    } // for c
+  } // for rl
+  
   bases.resize(h.reflevels());
   for (int rl=0; rl<h.reflevels(); ++rl) {
     if (h.components(rl)==0) {
@@ -221,6 +247,8 @@ void dh<D>::recompose () {
 	  bases[rl][ml].exterior *= boxes[rl][c][ml].exterior;
 	  bases[rl][ml].interior *= boxes[rl][c][ml].interior;
 	}
+	bases[rl][ml].boundaries
+	  = bases[rl][ml].exterior - bases[rl][ml].interior;
       }
     }
   }
@@ -247,6 +275,7 @@ void dh<D>::recompose () {
       	cout << "boundaries=" << boxes[rl][c][ml].boundaries << endl;
       	cout << "recv_sync=" << boxes[rl][c][ml].recv_sync << endl;
       	cout << "recv_ref_bnd_coarse=" << boxes[rl][c][ml].recv_ref_bnd_coarse << endl;
+      	cout << "sync_not=" << boxes[rl][c][ml].sync_not << endl;
       }
     }
   }
@@ -258,6 +287,7 @@ void dh<D>::recompose () {
       	cout << "rl=" << rl << " ml=" << ml << endl;
       	cout << "exterior=" << bases[rl][ml].exterior << endl;
       	cout << "interior=" << bases[rl][ml].interior << endl;
+      	cout << "boundaries=" << bases[rl][ml].boundaries << endl;
       }
     }
   }
diff --git a/Carpet/CarpetLib/src/dh.hh b/Carpet/CarpetLib/src/dh.hh
index 8960e60d1..91634d1ca 100644
--- a/Carpet/CarpetLib/src/dh.hh
+++ b/Carpet/CarpetLib/src/dh.hh
@@ -6,7 +6,7 @@
     copyright            : (C) 2000 by Erik Schnetter
     email                : schnetter@astro.psu.edu
 
-    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/dh.hh,v 1.2 2001/03/07 13:00:57 eschnett Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/dh.hh,v 1.3 2001/03/10 20:55:06 eschnett Exp $
 
  ***************************************************************************/
 
@@ -73,13 +73,15 @@ class dh {
     iblistvect send_ref_bnd_fine;
     
     ibset boundaries;		// boundaries
-    iblistvect recv_sync;	// receive while syncing
-    iblistvect recv_ref_bnd_coarse;
+    iblistvect recv_sync;	// received while syncing
+    iblistvect recv_ref_bnd_coarse; // received from coarser grids
+    ibset sync_not;		// not received while syncing (outer boundary)
   };
   
   struct dbases {
     ibbox exterior;		// whole region (including boundaries)
     ibbox interior;		// interior (without boundaries)
+    ibset boundaries;		// boundaries
   };
   
   typedef vector<dboxes> mboxes; // ... for each multigrid level
diff --git a/Carpet/CarpetLib/src/gdata.cc b/Carpet/CarpetLib/src/gdata.cc
index 36889de97..288eb9aa6 100644
--- a/Carpet/CarpetLib/src/gdata.cc
+++ b/Carpet/CarpetLib/src/gdata.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/gdata.cc,v 1.4 2001/03/07 13:00:57 eschnett Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gdata.cc,v 1.5 2001/03/10 20:55:06 eschnett Exp $
 
  ***************************************************************************/
 
@@ -20,6 +20,7 @@
 
 #include <cassert>
 #include <fstream>
+#include <iomanip>
 
 #include "bbox.hh"
 #include "defs.hh"
@@ -67,6 +68,8 @@ void generic_data<D>::write_ascii (const string name, const int time,
       ofstream file(name.c_str(), ios::app);
       assert (file.good());
       
+      file << setprecision(15);
+      
       file << "#" << endl
 	   << "# iteration " << time << endl
 	   << "# time level " << tl << "   refinement level " << rl
@@ -111,7 +114,8 @@ void generic_data<D>::write_ascii (const string name, const int time,
   } else {
     // copy to processor 0 and output there
     
-    generic_data* tmp = make_typed(_extent, 0);
+    generic_data* const tmp = make_typed();
+    tmp->allocate(_extent, 0);
     tmp->copy_from (this, _extent);
     tmp->write_ascii (name, time, org, dirs, tl, rl, c, ml);
     delete tmp;
diff --git a/Carpet/CarpetLib/src/gdata.hh b/Carpet/CarpetLib/src/gdata.hh
index 1d245be7c..4e77071d1 100644
--- a/Carpet/CarpetLib/src/gdata.hh
+++ b/Carpet/CarpetLib/src/gdata.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/gdata.hh,v 1.2 2001/03/05 14:31:03 eschnett Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gdata.hh,v 1.3 2001/03/10 20:55:06 eschnett Exp $
 
  ***************************************************************************/
 
@@ -30,7 +30,6 @@
 #include "defs.hh"
 #include "dist.hh"
 #include "bbox.hh"
-#include "bboxset.hh"
 #include "vect.hh"
 
 
@@ -51,12 +50,12 @@ class generic_data {
   // Types
   typedef vect<int,D> ivect;
   typedef bbox<int,D> ibbox;
-  typedef bboxset<int,D> ibset;
 
 protected:                      // should be readonly
 
   // Fields
   bool _has_storage;		// has storage associated (on some processor)
+  bool _owns_storage;		// owns the storage
   ivect _shape, _stride;      	// shape and index order
   int _size;			// size
 
@@ -73,21 +72,25 @@ public:
   virtual ~generic_data ();
 
   // Pseudo constructors
-  virtual generic_data* make_typed (const ibbox& extent, const int proc) const
-    = 0;
+  virtual generic_data* make_typed () const = 0;
 
   // Storage management
-  virtual void allocate (const ibbox& extent, const int proc) = 0;
+  virtual void allocate (const ibbox& extent, const int proc,
+			 void* const mem=0) = 0;
   virtual void free () = 0;
   virtual void transfer_from (generic_data* src) = 0;
 
   // Processor management
-  virtual void change_processor (const int newproc) = 0;
+  virtual void change_processor (const int newproc, void* const mem=0) = 0;
 
   // Accessors
   bool has_storage () const {
     return _has_storage;
   }
+  bool owns_storage () const {
+    assert (_has_storage);
+    return _owns_storage;
+  }
   
   virtual const void* storage () const = 0;
   
diff --git a/Carpet/CarpetSlab/src/carpetslab.cc b/Carpet/CarpetSlab/src/carpetslab.cc
index 15cb8a874..599b31bf3 100644
--- a/Carpet/CarpetSlab/src/carpetslab.cc
+++ b/Carpet/CarpetSlab/src/carpetslab.cc
@@ -1,25 +1,24 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetSlab/src/Attic/carpetslab.cc,v 1.3 2001/03/07 13:01:11 eschnett Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetSlab/src/Attic/carpetslab.cc,v 1.4 2001/03/10 20:55:09 eschnett Exp $
 
 #include <cassert>
 #include <cstdlib>
-#include <cstring>
-
-#include <mpi.h>
 
 #include "cctk.h"
 
 #include "Carpet/CarpetLib/src/bbox.hh"
-#include "Carpet/CarpetLib/src/gdata.hh"
+#include "Carpet/CarpetLib/src/bboxset.hh"
 #include "Carpet/CarpetLib/src/dh.hh"
-#include "Carpet/CarpetLib/src/dist.hh"
-#include "Carpet/CarpetLib/src/ggf.hh"
+#include "Carpet/CarpetLib/src/gdata.hh"
 #include "Carpet/CarpetLib/src/gh.hh"
+#include "Carpet/CarpetLib/src/ggf.hh"
 #include "Carpet/CarpetLib/src/vect.hh"
 
 #include "Carpet/Carpet/src/carpet.hh"
 
 #include "carpetslab.hh"
 
+static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetSlab/src/Attic/carpetslab.cc,v 1.4 2001/03/10 20:55:09 eschnett Exp $";
+
 
 
 namespace CarpetSlab {
@@ -28,81 +27,85 @@ namespace CarpetSlab {
   
   
   
-  int Hyperslab_GetLocalHyperslab (cGH* cgh,
-				   const int n,
-				   const int tl,
-				   const int hdim,
-				   const int origin [/*vdim*/],
-				   const int dir [/*vdim*/],
-				   const int len [/*hdim*/],
-				   const int downsample [/*hdim*/],
-				   void** const hdata,
-				   int hsize [/*hdim*/],
-				   int ghsize [/*hdim*/],
-				   int hoffset [/*hdim*/])
+  void* GetSlab (cGH* const cgh,
+		 const int dest_proc,
+		 const int n,
+		 const int tl,
+		 const int hdim,
+		 const int origin[/*vdim*/],
+		 const int dirs[/*hdim*/],
+		 const int stride[/*hdim*/],
+		 const int length[/*hdim*/])
   {
-    // check current status
-    assert (mglevel>=0);
+    // Check global state
     assert (reflevel>=0);
-    assert (component>=0);	// local := this component
+    assert (mglevel>=0);
+    
+    // Save global state
+    int saved_component = component;
+    component = -1;
     
-    // check arguments
+    // Check Cactus grid hierarchy
+    assert (cgh);
+    
+    // Check destination processor
+    assert (dest_proc>=-1 && dest_proc<CCTK_nProcs(cgh));
+    
+    // Check dimension
+    assert (hdim>=0 && hdim<=dim);
+    
+    // Check variable index
     assert (n>=0 && n<CCTK_NumVars());
-    assert (tl>=0);
-    assert (hdim>0 && hdim<=dim);
-    // the following assertion is too strict (better allow arbitrary values)
-    {
-      const int group = CCTK_GroupIndexFromVarI(n);
-      assert (group>=0);
-      switch (CCTK_GroupTypeFromVarI(n)) {
-      case CCTK_SCALAR:
-	abort();
-      case CCTK_ARRAY: {
-	assert (group<(int)arrdata.size());
-// 	assert (arrdata[group].hh->is_local(reflevel, component));
-// 	const bbox<int,dim> ext =
-// 	  arrdata[group].dd->boxes[reflevel][component][mglevel].exterior;
-// 	for (int d=0; d<dim; ++d) {
-// 	  assert (origin[d] >= ext.lower()[d] / ext.stride()[d]
-// 		  && origin[d] <= ext.upper()[d] / ext.stride()[d]);
-// 	}
-	break;
-      }
-      case CCTK_GF: {
-	assert (group<(int)gfdata.size());
-// 	assert (hh->is_local(reflevel, component));
-// 	const bbox<int,dim> ext =
-// 	  dd->boxes[reflevel][component][mglevel].exterior;
-// 	for (int d=0; d<dim; ++d) {
-// 	  assert (origin[d] >= ext.lower()[d] / ext.stride()[d]
-// 		  && origin[d] <= ext.upper()[d] / ext.stride()[d]);
-// 	}
-	break;
-      }
-      default:
-	abort();
-      }
-    }
-    // the following assertion is too strict (better allow arbitrary values)
-    for (int d=0; d<dim; ++d) assert (dir[d]==0 || dir[d]==1);
-    // the following assertion is too strict (better allow arbitrary values)
-    for (int d=0; d<hdim; ++d) assert (downsample[d]>0);
-    assert (hdata);
-    assert (hsize);
-    assert (ghsize);
-    assert (hoffset);
     
-    // get variable info
+    // Get info about variable
     const int group = CCTK_GroupIndexFromVarI(n);
     assert (group>=0);
-    const int var = n - CCTK_FirstVarIndexI(group);
+    const int n0 = CCTK_FirstVarIndexI(group);
+    assert (n0>=0);
+    const int var = n - n0;
     assert (var>=0);
     
-    // get grid hierarchy data hierarchy and grid function
-    gh<dim>* myhh;
-    dh<dim>* mydd;
-    generic_gf<dim>* myff;
-    switch (CCTK_GroupTypeFromVarI(n)) {
+    // Get info about group
+    cGroup gp;
+    CCTK_GroupData (group, &gp);
+    assert (gp.dim==dim);
+    assert (CCTK_QueryGroupStorageI(cgh, group));
+    const int typesize = CCTK_VarTypeSize(gp.vartype);
+    assert (typesize>0);
+    
+    // Check timelevel
+    assert (tl>=0 && tl<gp.numtimelevels);
+    
+    // Check origin
+//     for (int d=0; d<dim; ++d) {
+//       assert (origin[d]>=0 && origin[d]<=sizes[d]);
+//     }
+    
+    // Check directions
+    for (int dd=0; dd<hdim; ++dd) {
+      assert (dirs[dd]>=1 && dirs[dd]<=dim);
+    }
+    
+    // Check stride
+    for (int dd=0; dd<hdim; ++dd) {
+      assert (stride[dd]>0);
+    }
+    
+    // Check length
+    for (int dd=0; dd<hdim; ++dd) {
+      assert (length[dd]>=0);
+    }
+    
+    // Check extent
+//     for (int dd=0; dd<hdim; ++dd) {
+//       assert (origin[dirs[dd]-1] + length[dd] <= sizes[dirs[dd]]);
+//     }
+    
+    // Get insider information about variable
+    const gh<dim>* myhh;
+    const dh<dim>* mydd;
+    const generic_gf<dim>* myff;
+    switch (gp.grouptype) {
     case CCTK_SCALAR:
       abort();
     case CCTK_ARRAY:
@@ -114,7 +117,7 @@ namespace CarpetSlab {
       break;
     case CCTK_GF:
       myhh = hh;
-      mydd = dd;
+      mydd = Carpet::dd;	// necessary for pre-ANSI C++ compilers
       assert (group < (int)gfdata.size());
       assert (var < (int)gfdata[group].data.size());
       myff = gfdata[group].data[var];
@@ -122,388 +125,243 @@ namespace CarpetSlab {
     default:
       abort();
     }
+    assert (myhh);
+    assert (mydd);
+    assert (myff);
     
-    // get data
-    const generic_data<dim>* mydata
-      = (*myff)(tl, reflevel, component, mglevel);
-    
-    // get local and global bounding boxes
-    assert (reflevel < (int)mydd->boxes.size());
-    assert (component < (int)mydd->boxes[reflevel].size());
-    assert (mglevel < (int)mydd->boxes[reflevel][component].size());
-    const bbox<int,dim> locbox
-      = mydd->boxes[reflevel][component][mglevel].exterior;
-    const bbox<int,dim> globox
-      = mydd->bases[reflevel][mglevel].exterior;
-    
-    // calculate more convenient representation of the direction
-    vect<int,dim> stride[hdim];
-    // the following if statement is written according to the
-    // definition of "dir".
-    if (hdim==1) {
-      for (int d=0; d<dim; ++d) stride[0][d] = dir[d] * downsample[0];
-    } else if (hdim==dim) {
-      for (int dd=0; dd<hdim; ++dd) {
-	for (int d=0; d<dim; ++d) stride[dd][d] = d==dd ? downsample[dd] : 0;
-      }
-    } else if (hdim==2) {
-      assert (dim==3);
-      if (dir[0]==0) {
-	assert (dir[1]!=0 && dir[2]!=0);
-	stride[0] = vect<int,dim>::dir(1);
-	stride[1] = vect<int,dim>::dir(2);
-      } else if (dir[1]==0) {
-	assert (dir[0]!=0 && dir[2]!=0);
-	stride[0] = vect<int,dim>::dir(0);
-	stride[1] = vect<int,dim>::dir(2);
-      } else if (dir[2]==0) {
-	assert (dir[0]!=0 && dir[1]!=0);
-	stride[0] = vect<int,dim>::dir(0);
-	stride[1] = vect<int,dim>::dir(1);
-      } else {
-	abort();
-      }
-      for (int dd=0; dd<hdim; ++dd) stride[dd] *= downsample[dd];
-    } else {
-      abort();
-    }
-    for (int dd=0; dd<hdim; ++dd) stride[dd] *= locbox.stride();
-    
-    // local lower bound
-    vect<int,dim> lbound;
-    for (int d=0; d<dim; ++d) lbound[d] = origin[d];
-    lbound *= locbox.stride();
-    
-    // local upper bound
-    vect<int,dim> ubound = lbound;
-    for (int dd=0; dd<hdim; ++dd) {
-      if (len[dd]<0) {
-	assert (any(stride[dd]>0));
-	while (all(ubound + stride[dd] <= locbox.upper())) {
-	  ubound += stride[dd];
-	}
-      } else {
-	ubound += stride[dd] * (len[dd]-1);
-      }
-    }
-    
-    lbound = max(lbound, locbox.lower());
-    ubound = min(ubound, locbox.upper());
-    
-    // local size
-    int total_hsize = 1;
-    for (int dd=0; dd<hdim; ++dd) {
-      hsize[dd] = 0;
-      assert (any(stride[dd]>0));
-      while (all(lbound + stride[dd] * hsize[dd] <= ubound)) ++hsize[dd];
-      total_hsize *= hsize[dd];
-    }
-    assert (total_hsize>=0);
-    
-    // sanity check
-    if (total_hsize>0) {
-      vect<int,dim> tmp = lbound;
-      for (int dd=0; dd<hdim; ++dd) {
-	tmp += stride[dd] * (hsize[dd]-1);
-      }
-      assert (all(tmp == ubound));
-    }
+    // Detemine collecting processor
+    const int collect_proc = dest_proc<0 ? 0 : dest_proc;
     
-    // global lower bound
-    vect<int,dim> glbound;
-    for (int d=0; d<dim; ++d) glbound[d] = origin[d];
-    glbound *= globox.stride();
+    // Determine own rank
+    const int rank = CCTK_MyProc(cgh);
     
-    // global upper bound
-    vect<int,dim> gubound = glbound;
+    // Calculate global size
+    int totalsize = 1;
     for (int dd=0; dd<hdim; ++dd) {
-      if (len[dd]<0) {
-	assert (any(stride[dd]>0));
-	while (all(gubound + stride[dd] <= globox.upper())) {
-	  gubound += stride[dd];
-	}
-      } else {
-	gubound += stride[dd] * (len[dd]-1);
-      }
+      totalsize *= length[dd];
     }
     
-    glbound = max(glbound, globox.lower());
-    gubound = min(gubound, globox.upper());
-    
-    // global size
-    int total_ghsize = 1;
-    for (int dd=0; dd<hdim; ++dd) {
-      ghsize[dd] = 0;
-      assert (any(stride[dd]>0));
-      while (all(glbound + stride[dd] * ghsize[dd] <= gubound)) {
-	++ghsize[dd];
-      }
-      total_ghsize *= ghsize[dd];
+    // Allocate memory
+    void* hdata = 0;
+    if (dest_proc==-1 || rank==dest_proc) {
+      hdata = malloc(totalsize * typesize);
+      assert (hdata);
     }
-    assert (total_ghsize>=0);
     
-    // sanity check
-    if (total_ghsize>0) {
-      vect<int,dim> tmp = glbound;
+    if (hh->components(reflevel) > 0) {
+      
+      // Only temporarily
+      component = 0;
+      
+      // Get sample data
+      const generic_data<dim>* mydata;
+      mydata = (*myff)(tl, reflevel, component, mglevel);
+      
+      // Stride of data in memory
+      const vect<int,dim> str = mydata->extent().stride();
+      
+      // Stride of collected data
+      vect<int,dim> hstr = str;
       for (int dd=0; dd<hdim; ++dd) {
-	tmp += stride[dd] * (ghsize[dd]-1);
-      }
-      assert (all(tmp == gubound));
-    }
-    
-    // local to global offset
-    for (int dd=0; dd<hdim; ++dd) {
-      hoffset[dd] = 0;
-      assert (any(stride[dd]>0));
-      while (all(glbound + stride[dd] * (hoffset[dd]+1) <= lbound)) {
-	++hoffset[dd];
+	hstr[dirs[dd]-1] *= stride[dd];
       }
-    }
-    
-    // sanity check
-    {
-      vect<int,dim> tmp = glbound;
-      for (int dd=0; dd<hdim; ++dd) {
-	tmp += stride[dd] * hoffset[dd];
+      
+      // Lower bound of collected data
+      vect<int,dim> hlb;
+      for (int d=0; d<dim; ++d) {
+	hlb[d] = origin[d] * str[d];
       }
-      assert (all(tmp == lbound));
-    }
-    
-    // bail out if this component is on another processor
-    if (! myhh->is_local(reflevel, component)) {
-      *hdata = 0;
+      
+      // Upper bound of collected data
+      vect<int,dim> hub = hlb;
       for (int dd=0; dd<hdim; ++dd) {
-	hsize[dd] = 0;
-	hoffset[dd] = 0;
+	hub[dirs[dd]-1] += (length[dd]-1) * hstr[dirs[dd]-1];
       }
-      return -1;
-    }
-    
-    // allocate the memory
-    *hdata = malloc(total_hsize * CCTK_VarTypeSize(CCTK_VarTypeI(n)));
-    assert (*hdata);
-    
-    if (total_hsize>0) {
       
-      // copy the data to user memory
-      char* const       dest = (char*)*hdata;
-      const char* const src  = (const char*)mydata->storage();
-      const int         sz   = CCTK_VarTypeSize(CCTK_VarTypeI(n));
+      // Calculate extent to collect
+      const bbox<int,dim> hextent (hlb, hub, hstr);
+      assert (hextent.num_points() == totalsize);
       
-      int dest_index[hdim];
-      for (int dd=0; dd<hdim; ++dd) dest_index[dd] = 0;
-      for (;;) {
-	
-	vect<int,dim> src_index = lbound;
-	for (int dd=0; dd<hdim; ++dd) src_index += stride[dd] * dest_index[dd];
+      // Create collector data object
+      void* myhdata = rank==collect_proc ? hdata : 0;
+      generic_data<dim>* const alldata = mydata->make_typed();
+      alldata->allocate (hextent, collect_proc, myhdata);
+      
+      // Done with the temporary stuff
+      mydata = 0;
+      component = -1;
+      
+      // Loop over all components, copying data from them
+      assert (component == -1);
+      for (component=0; component<hh->components(reflevel); ++component) {
 	
-	int di = 0;
-	for (int dd=hdim-1; dd>=0; --dd) di = di * hsize[dd] + dest_index[dd];
+	// Get data object
+	mydata = (*myff)(tl, reflevel, component, mglevel);
 	
-	const int si = mydata->offset(src_index);
+	// Calculate overlapping extents
+	const bboxset<int,dim> myextents
+	  = ((mydd->boxes[reflevel][component][mglevel].sync_not
+	      | mydd->boxes[reflevel][component][mglevel].interior)
+	     & hextent);
 	
-	memcpy(dest + sz*di, src + sz*si, sz);
+	// Loop over overlapping extents
+	for (bboxset<int,dim>::const_iterator ext_iter = myextents.begin();
+	     ext_iter != myextents.end(); ++ext_iter) {
+	  
+	  // Copy data
+	  alldata->copy_from (mydata, *ext_iter);
+	  
+	}
 	
-	for (int dd=0; dd<hdim; ++dd) {
-	  ++dest_index[dd];
-	  if (dest_index[dd]<hsize[dd]) break;
-	  dest_index[dd]=0;
-	  if (dd==hdim-1) goto done;
+      }	// Loop over components
+      component = -1;
+      
+      // Copy result to all processors
+      if (dest_proc == -1) {
+	for (int proc=0; proc<CCTK_nProcs(cgh); ++proc) {
+	  if (proc != collect_proc) {
+	    
+	    void* myhdata = rank==proc ? hdata : 0;
+	    generic_data<dim>* const tmpdata = mydata->make_typed();
+	    tmpdata->allocate (alldata->extent(), proc, myhdata);
+	    tmpdata->copy_from (alldata, alldata->extent());
+	    delete tmpdata;
+	    
+	  }
 	}
-      }
-    done: ;
+      }	// Copy result
+      
+      delete alldata;
       
-    } // if total_hsize>0
+    } // if components>0
     
-    return 0;
+    // Restore global state
+    component = saved_component;
+    
+    // Success
+    return hdata;
   }
   
   
   
-  int Hyperslab_GetHyperslab (cGH* cgh,
+  int Hyperslab_GetHyperslab (cGH* const GH,
 			      const int target_proc,
-			      const int n,
-			      const int tl,
+			      const int vindex,
+			      const int vtimelvl,
 			      const int hdim,
-			      const int origin [/*vdim*/],
-			      const int dir [/*vdim*/],
-			      const int len [/*hdim*/],
+			      const int global_startpoint [/*vdim*/],
+			      const int directions [/*vdim*/],
+			      const int lengths [/*hdim*/],
 			      const int downsample [/*hdim*/],
 			      void** const hdata,
 			      int hsize [/*hdim*/])
   {
-    // check current status
-    assert (mglevel>=0);
-    assert (reflevel>=0);
-    
-    // in order to work, this requires that all processors have the
-    // same number of components
-    const int saved_component = component;
-    component = -1;
+    // Check some arguments
+    assert (hdim>=0 && hdim<=dim);
     
-    // check arguments
-    assert (n>=0 && n<CCTK_NumVars());
-    assert (tl>=0);
-    assert (hdim>0 && hdim<=dim);
-    // the following assertion is too strict (better allow arbitrary values)
-    {
-      // TODO: make sure that origin is within the extent of this
-      // refinement / multigrid level
-      // (but no such extent is stored in dh)
-      // (it is now; use it!)
-      const int group = CCTK_GroupIndexFromVarI(n);
-      assert (group>=0);
-      switch (CCTK_GroupTypeFromVarI(n)) {
-      case CCTK_SCALAR:
-	abort();
-      case CCTK_ARRAY:
-	assert (group<(int)arrdata.size());
-	break;
-      case CCTK_GF:
-	assert (group<(int)gfdata.size());
-	break;
-      default:
-	abort();
-      }
-    }
-    // the following assertion is too strict (better allow arbitrary values)
-    for (int d=0; d<dim; ++d) assert (dir[d]==0 || dir[d]==1);
-    // the following assertion is too strict (better allow arbitrary values)
-    for (int d=0; d<hdim; ++d) assert (downsample[d]>0);
+    // Check output arguments
     assert (hdata);
     assert (hsize);
     
-    int collect_proc = target_proc;
-    if (collect_proc<0) collect_proc = 0;
-    
-    assert (hh->components(reflevel)>0);
-    *hdata = 0;
-    for (int dd=0; dd<hdim; ++dd) hsize[dd] = 0;
-    int totalhsize = 0;
-    
-    const int sz = CCTK_VarTypeSize(CCTK_VarTypeI(n));
-    
-    MPI_Datatype type;
-    switch (CCTK_VarTypeI(n)) {
-#define TYPECASE(N,T)				\
-    case N: {					\
-      assert (sz == sizeof(T));			\
-      T dummy;					\
-      type = dist::datatype(dummy);		\
-      break;					\
-    }
-#include "Carpet/Carpet/src/typecase"
-#undef TYPECASE
-    default:
-      abort();
-    }
-    
-    int rank;
-    MPI_Comm_rank (dist::comm, &rank);
-    
-    // loop over all components
-    for (component=0; component<hh->components(reflevel); ++component) {
-      
-      void* myhdata;
-      int myhsize[hdim], ghsize[hdim], hoffset[hdim];
-      
-      const int retval = Hyperslab_GetLocalHyperslab
-	(cgh, n, tl, hdim, origin, dir, len, downsample,
-	 &myhdata, myhsize, ghsize, hoffset);
-      
-      if (hh->is_local(reflevel,component)) {
-	assert (retval == 0);
-      } else {
-	assert (retval == -1);
+    // Calculate more convenient representation of the direction
+    int dirs[hdim];
+    // The following if statement is written according to the
+    // definition of "dir".
+    if (hdim==1) {
+      // 1-dimensional hyperslab
+      int mydir = 0;
+      for (int d=0; d<dim; ++d) {
+	if (directions[d]!=0) {
+	  mydir = d+1;
+	  break;
+	}
       }
-      
-      int mytotalsize = 1;
-      for (int dd=0; dd<hdim; ++dd) mytotalsize *= myhsize[dd];
-      
-      if (component==0) {
-	if (target_proc<0 || rank == target_proc) {
-	  
-	  for (int dd=0; dd<hdim; ++dd) hsize[dd] = ghsize[dd];
-	  
-	  totalhsize = 1;
-	  for (int dd=0; dd<hdim; ++dd) totalhsize *= hsize[dd];
-	  
-	  if (rank == collect_proc) {
-	    *hdata = malloc(totalhsize * sz);
-	    assert (*hdata);
-	  } else {
-	    *hdata = 0;
-	  }
-	
+      assert (mydir>0);
+      for (int d=0; d<dim; ++d) {
+	if (d == mydir-1) {
+	  assert (directions[d]!=0);
+	} else {
+	  assert (directions[d]==0);
 	}
       }
-      
-      if (!myhdata && rank == collect_proc) {
-	MPI_Status status;
-	assert (mytotalsize==0);
-	const int src = hh->proc(reflevel, component);
-	MPI_Recv (&mytotalsize, 1, MPI_INT, src, 2001,
-		  dist::comm, &status);
-	myhdata = malloc(mytotalsize * sz);
-	assert (myhdata);
-	MPI_Recv (myhdata, mytotalsize, type, src, 2001,
-		  dist::comm, &status);
-      } else if (myhdata && rank != collect_proc) {
-	MPI_Send (&mytotalsize, 1, MPI_INT, collect_proc, 2001, dist::comm);
-	MPI_Send (myhdata, mytotalsize, type, collect_proc, 2001, dist::comm);
-	free (myhdata);
-	myhdata = 0;
-	mytotalsize = 0;
+      dirs[0] = mydir;
+    } else if (hdim==dim) {
+      // dim-dimensional hyperslab
+      for (int dd=0; dd<hdim; ++dd) {
+	dirs[dd] = dd+1;
       }
-      
-      if (myhdata>0) {
-	assert (rank == collect_proc);
-	
-	if (mytotalsize>0) {
-	  int dest_index[hdim], src_index[hdim];
-	  for (int dd=0; dd<hdim; ++dd) dest_index[dd] = hoffset[dd];
-	  for (int dd=0; dd<hdim; ++dd) src_index[dd] = 0;
-	  for (;;) {
-	    int di=0;
-	    for (int dd=hdim-1; dd>=0; --dd) {
-	      di = di * hsize[dd] + dest_index[dd];
-	    }
-	    int si=0;
-	    for (int dd=hdim-1; dd>=0; --dd) {
-	      si = si * myhsize[dd] + src_index[dd];
-	    }
-	    
-	    memcpy ((char*)*hdata + sz*di, (char*)myhdata + sz*si, sz);
-	    
-	    for (int dd=0; dd<hdim; ++dd) {
-	      ++dest_index[dd];
-	      ++src_index[dd];
-	      if (src_index[dd] < myhsize[dd]) break;
-	      dest_index[dd] = hoffset[dd];
-	      src_index[dd] = 0;
-	      if (dd==hdim-1) goto done;
-	    }
-	  }
-	done: ;
-	  
-	} // if mytotalsize>0
-	
-	free (myhdata);
-	myhdata = 0;
-      } else {
-	assert (rank != collect_proc);
+    } else if (hdim==2) {
+      // 2-dimensional hyperslab with dim==3
+      assert (dim==3);
+      int mydir = 0;
+      for (int d=0; d<dim; ++d) {
+	if (directions[d]==0) {
+	  mydir = d+1;
+	  break;
+	}
       }
-      
+      assert (mydir>0);
+      for (int d=0; d<dim; ++d) {
+	if (d == mydir-1) {
+	  assert (directions[d]==0);
+	} else {
+	  assert (directions[d]!=0);
+	}
+      }
+      int dd=0;
+      for (int d=0; d<dim; ++d) {
+	if (d != mydir-1) {
+	  dirs[dd] = d+1;
+	  ++dd;
+	}
+      }
+      assert (dd==hdim);
+    } else {
+      abort();
     }
     
-    if (target_proc<0) {
-      MPI_Bcast (hdata, totalhsize*sz, MPI_BYTE,
-		 collect_proc, CarpetMPICommunicator());
+#if 0
+    // Invert directions
+    int invdir[dim];
+    for (int d=0; d<dim; ++d) {
+      invdir[d] = 0;
     }
+    for (int dd=0; dd<hdim; ++dd) {
+      invdir[dirs[dd]-1] = dd+1;
+    }
+#endif
     
-    component = saved_component;
+    // Calculate lengths
+    for (int dd=0; dd<hdim; ++dd) {
+      if (lengths[dd]<0) {
+	const int totlen = *CCTK_ArrayGroupSizeI(GH, dirs[dd]-1, vindex);
+	assert (totlen>=0);
+	// Partial argument check
+	assert (global_startpoint[dirs[dd]-1]>=0);
+	assert (global_startpoint[dirs[dd]-1]<=totlen);
+	assert (downsample[dd]>0);
+	hsize[dd] = (totlen - global_startpoint[dirs[dd]-1]) / downsample[dd];
+      } else {
+	hsize[dd] = lengths[dd];
+      }
+      assert (hsize[dd]>=0);
+    }
     
-    return 0;
+    // Get the slab
+    *hdata = GetSlab (GH,
+		      target_proc,
+		      vindex,
+		      vtimelvl,
+		      hdim,
+		      global_startpoint,
+		      dirs,
+		      downsample,
+		      hsize);
+    
+    // Return with success
+    return 1;
   }
   
+  
+  
 } // namespace CarpetSlab
diff --git a/Carpet/CarpetSlab/src/carpetslab.h b/Carpet/CarpetSlab/src/carpetslab.h
index dd3f86136..6e3fd7f00 100644
--- a/Carpet/CarpetSlab/src/carpetslab.h
+++ b/Carpet/CarpetSlab/src/carpetslab.h
@@ -1,28 +1,28 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetSlab/src/Attic/carpetslab.h,v 1.1 2001/03/07 13:01:11 eschnett Exp $
+/* $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetSlab/src/Attic/carpetslab.h,v 1.2 2001/03/10 20:55:09 eschnett Exp $ */
 
-#include "cctk.h"
+#ifndef CARPETSLAB_H
+#define CARPETSLAB_H
 
-int Hyperslab_GetLocalHyperslab (cGH* GH,
-				 int vindex,
-				 int vtimelvl,
-				 int hdim,
-				 const int global_startpoint [/*vdim*/],
-				 const int directions [/*vdim*/],
-				 const int lengths [/*hdim*/],
-				 const int downsample [/*hdim*/],
-				 void** hdata,
-				 int hsize [/*hdim*/],
-				 int ghsize [/*hdim*/],
-				 int hoffset [/*hdim*/]);
+#ifdef __cplusplus
+namespace CarpetSlab {
+  extern "C" {
+#endif
+    
+    int Hyperslab_GetHyperslab (cGH* const GH,
+				const int target_proc,
+				const int vindex,
+				const int vtimelvl,
+				const int hdim,
+				const int global_startpoint [/*vdim*/],
+				const int directions [/*vdim*/],
+				const int lengths [/*hdim*/],
+				const int downsample [/*hdim*/],
+				void** const hdata,
+				int hsize [/*hdim*/]);
+    
+#ifdef __cplusplus
+  } /* extern "C" */
+} /* namespace CarpetSlab */
+#endif
 
-int Hyperslab_GetHyperslab (cGH* GH,
-			    int target_proc,
-			    int vindex,
-			    int vtimelvl,
-			    int hdim,
-			    const int global_startpoint [/*vdim*/],
-			    const int directions [/*vdim*/],
-			    const int lengths [/*hdim*/],
-			    const int downsample [/*hdim*/],
-			    void** hdata,
-			    int hsize [/*hdim*/]);
+#endif /* !defined(CARPETSLAB_H) */
diff --git a/Carpet/CarpetSlab/src/carpetslab.hh b/Carpet/CarpetSlab/src/carpetslab.hh
index 444d14cf1..48fa52c31 100644
--- a/Carpet/CarpetSlab/src/carpetslab.hh
+++ b/Carpet/CarpetSlab/src/carpetslab.hh
@@ -1,36 +1,22 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetSlab/src/Attic/carpetslab.hh,v 1.2 2001/03/07 13:01:11 eschnett Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetSlab/src/Attic/carpetslab.hh,v 1.3 2001/03/10 20:55:09 eschnett Exp $
 
-#include "cctk.h"
+#ifndef CARPETSLAB_HH
+#define CARPETSLAB_HH
+
+#include "carpetslab.h"
 
 namespace CarpetSlab {
   
-  extern "C" {
-    
-    int Hyperslab_GetLocalHyperslab (cGH* GH,
-				     int vindex,
-				     int vtimelvl,
-				     int hdim,
-				     const int global_startpoint [/*vdim*/],
-				     const int directions [/*vdim*/],
-				     const int lengths [/*hdim*/],
-				     const int downsample [/*hdim*/],
-				     void** hdata,
-				     int hsize [/*hdim*/],
-				     int ghsize [/*hdim*/],
-				     int hoffset [/*hdim*/]);
-    
-    int Hyperslab_GetHyperslab (cGH* GH,
-				int target_proc,
-				int vindex,
-				int vtimelvl,
-				int hdim,
-				const int global_startpoint [/*vdim*/],
-				const int directions [/*vdim*/],
-				const int lengths [/*hdim*/],
-				const int downsample [/*hdim*/],
-				void** hdata,
-				int hsize [/*hdim*/]);
-    
-  } // extern "C"
+  void* GetSlab (cGH* const cgh,
+		 const int dest_proc,
+		 const int n,
+		 const int tl,
+		 const int hdim,
+		 const int origin[/*vdim*/],
+		 const int dirs[/*hdim*/],
+		 const int stride[/*hdim*/],
+		 const int length[/*hdim*/]);
   
 } // namespace CarpetSlab
+
+#endif // !defined(CARPETSLAB_HH)