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

3 files changed, 309 insertions, 465 deletions

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)