n-dim. extraction out of n-dim volume

git-svn-id: http://svn.cactuscode.org/arrangements/CactusPUGH/PUGHSlab/trunk@19 10716dce-81a3-4424-a2c8-48026a0d3035

9 files changed, 261 insertions, 322 deletions

diff --git a/src/CollectData1D.c b/src/CollectData1D.c
index 211e77b..3605423 100644
--- a/src/CollectData1D.c
+++ b/src/CollectData1D.c
@@ -8,13 +8,10 @@
 #include "cctk_Arguments.h"
 
 #include "Hyperslab.h"
+#include "Hyperslabi.h"
 #include "CactusBase/IOUtil/src/ioGH.h"
 #include "CactusPUGH/PUGH/src/include/pugh.h"
 
-/* MPI message tag */
-#define TAGBASE 1234
-#define MAX_DIM 3
-
 
 int CollectData1D(cGH *GH, int vindex, int vtimelvl, int vdim,
                   int *S_l, int *u_l, int *ds, int len_l, 
@@ -146,10 +143,10 @@ int CollectLocalData1D(cGH *GH, int vindex, int vtimelvl, int vdim,
 {
   int *locstart, *locend, nlocpoints;
   int ipnt, gridpnt[MAX_DIM];
-  int ret_val;
+  int ierr;
 
-  int gindex,vtype,vtypesize;
-  int idim,lincount;
+  int vtype,vtypesize;
+  int idim, lincount;
 
   void *data = CCTK_VarDataPtrI (GH, vtimelvl, vindex);
 
@@ -163,21 +160,22 @@ int CollectLocalData1D(cGH *GH, int vindex, int vtimelvl, int vdim,
   locstart = (int*)malloc(vdim*sizeof(int));
   locend   = (int*)malloc(vdim*sizeof(int));
 
-  ret_val  = GetLocalLine(GH, vdim, S_l, u_l, ds[0], len_l, 
+  ierr  = GetLocalLine(GH, vdim, S_l, u_l, ds[0], len_l, 
                           locstart, locend, &nlocpoints);
 
+  if (ierr<0)
+  {
+    CCTK_WARN (1, "GetLocalLine() failed");
+    return (ierr);
+  }
+
+  /* Transform from global to local index coordinate */
   for(idim=0;idim<vdim;idim++)
   {  
     locstart[idim]-=GH->cctk_lbnd[idim];
     locend[idim]  -=GH->cctk_lbnd[idim];
   }
 
-  if (ret_val<0)
-  {
-    CCTK_WARN (1, "GetLocalLine() failed");
-    return (ret_val);
-  }
-
   dsize[0]= nlocpoints;
 
   printf("CollectLocal1D: %d npoints   lb: %d %d %d\n",
@@ -256,11 +254,11 @@ int CollectLocalData1D(cGH *GH, int vindex, int vtimelvl, int vdim,
     CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
                 "Unsupported type for variable '%s'", fullname);
     free (fullname);
-    ret_val = -1;
+    ierr = -1;
   }
 
   free(locstart);
   free(locend);
 
-  return (ret_val);
+  return (ierr);
 }
diff --git a/src/CollectData2D.c b/src/CollectData2D.c
index a42ae8a..3019a1e 100644
--- a/src/CollectData2D.c
+++ b/src/CollectData2D.c
@@ -9,18 +9,15 @@
 #include "cctk_Arguments.h"
 
 #include "Hyperslab.h"
+#include "Hyperslabi.h"
 #include "CactusBase/IOUtil/src/ioGH.h"
 #include "CactusPUGH/PUGH/src/include/pugh.h"
 
 #define MIN(a,b) ((a)<(b) ? (a) : (b))
 #define MAX(a,b) ((a)>(b) ? (a) : (b))
 
-#define MAX_DIM 3
-#define BAD -42
-
-
-int CollectData2D(cGH *GH, int vindex, int vtimelvl, int vdim,
-		  int *S_s, int dir[2], int ds[2], int len[2],
+int CollectDataND(cGH *GH, int vindex, int vtimelvl, int vdim, int sdim,
+		  int *S_s, int *dir, int *ds, int *len,
 		  void **dptr, int *dsize)
 {
   return(0);
@@ -28,7 +25,7 @@ int CollectData2D(cGH *GH, int vindex, int vtimelvl, int vdim,
 
 
  /*@@
-   @routine    CollectLocalData2D
+   @routine    CollectLocalDataND
    @date       Fri May 12 17:51:30 2000
    @author     Gerd Lanfermann
    @desc 
@@ -46,30 +43,22 @@ int CollectData2D(cGH *GH, int vindex, int vtimelvl, int vdim,
 @@*/
 
 
-int CollectLocalData2D(cGH *GH, int vindex, int vtimelvl, int vdim,
-		  int *S_s, int linedir[2], int ds[2], int len[2],
-		  void **loc_dptr, int *idx1_dptr, int *idx2_dptr, int *dsize)
+int CollectLocalDataND(cGH *GH, int vindex, int vtimelvl, int vdim, int sdim,
+		       int *S_s, int *linedir, int *ds, int *len,
+		       void **loc_dptr, int *dirsize,int *totsize)
 {
 
-  int tstart[MAX_DIM];     /* temp start points */  
-  int tlocstart[MAX_DIM]; 
-  int tlocend[MAX_DIM];
-
-
-
-  int *locstart, *locend;          /* local start/end pnt in glob. index coord */
-  int nlocpoints[2], totlocpoints; /* number of loc. points in each dir, tot. pnts */
-
-  int u_l[2][MAX_DIM];             /* the 2 vectors which span the surface: 
-				      e.g. 001,100 */
+  int *locstart, *locend, nlocpoints;
+  int ipnt[MAX_DIM], gridpnt[MAX_DIM];
 
   int vtype,vtypesize;
-  int idim, idir;
+  int idim, idir, lincount;
   int ierr;
-  int phasplane=1;
-
+ 
   void *data = CCTK_VarDataPtrI (GH, vtimelvl, vindex);
-  CCTK_REAL *real_dptr = (CCTK_REAL *) *loc_dptr;
+
+  if (vdim > MAX_DIM)
+    return (-1);
 
   vtype      = CCTK_VarTypeI(vindex);
   vtypesize  = CCTK_VarTypeSize(vtype);
@@ -77,88 +66,50 @@ int CollectLocalData2D(cGH *GH, int vindex, int vtimelvl, int vdim,
   locstart = (int*)malloc(vdim*sizeof(int));
   locend   = (int*)malloc(vdim*sizeof(int));
 
-
-  /* Initialize the local start/end to BAD values and hope for the best... */
-  for (idim=0; idim<vdim; idim++) 
+  ierr = GetLocalSlab(GH, vdim, sdim, 
+		      S_s, linedir, ds, len, 
+		      locstart, locend, dirsize, &totsize);
+  if (ierr<0)
   {
-    tstart[idim]   = BAD;
-    locstart[idim] = BAD;
-    locend[idim]   = BAD;
-
-    /* line vectors are 000, 000 */
-    u_l[0][idim]  = 0;
-    u_l[1][idim]  = 0;
-
-    /* make sure that all procs above and below the plane, 
-       know that they don't have a plane */
-    if ((idim!=linedir[0]) && (idim!=linedir[1]))
-      if((GH->cctk_lbnd[idim]<S_s[idim]) && 
-	 (GH->cctk_ubnd[idim]>S_s[idim])) 
-      {
-	locstart[idim] = S_s[idim];
-	locend[idim]   = S_s[idim];
-	phasplane = 1;
-      }
-      else phasplane = 0;
+    CCTK_WARN (1, "GetLocalSlab() failed");
+    return (ierr);
+  }
+
+  /* Transform from global to local index coordinate */
+  for(idim=0;idim<vdim;idim++)
+  {  
+    locstart[idim]-=GH->cctk_lbnd[idim];
+    locend[idim]  -=GH->cctk_lbnd[idim];
   }
 
-  /* construct line vectors that span the surface:
-     e.g. linedir[0,1] = 1,2 means YZ surface:  (000,000) -> (010,001) */
-  u_l[0][linedir[0]]=1;
-  u_l[1][linedir[1]]=1;
+  *loc_dptr= malloc(nlocpoints * vtypesize);
+  lincount=0; 
 
-  /* init to 1 (neutral element for multiplication) */
-  totlocpoints=1;
-  
-  for (idir=0;idir<2;idir++) 
-  {
+  /* zero out index array so that unused dimensions are 0 */
+  memset (gridpnt, 0, sizeof (gridpnt));
+  memset (ipnt   , 0, sizeof (ipnt));
 
-    if (S_s[linedir[idir]]<GH->cctk_ubnd[linedir[idir]])
-    {
-      /* Initialize startpoints */
-      for (idim=0; idim<vdim; idim++) 
-      {
-	tstart[idim]   = GH->cctk_lbnd[idim]+GH->cctk_nghostzones[idim];
-      }
-      tstart[linedir[idir]]=S_s[linedir[idir]];
-      ierr = GetLocalLine(GH, vdim, tstart, u_l[idir], ds[idir], len[idir], 
-			  tlocstart, tlocend, &nlocpoints[idir]);
-      locstart[linedir[idir]]= tlocstart[linedir[idir]];
-      locend  [linedir[idir]]= tlocend[linedir[idir]];
-
-      totlocpoints*=nlocpoints[idir];
-
-      if (locstart[linedir[idir]]==BAD) phasplane &=0;
-    }
-    else 
-    {
-      phasplane &=0;
-      totlocpoints = 0;
-    }
-    
-  }
- 
-  if (phasplane)
-    printf("LOC START:  %d %d %d   %d %d %d   -> %d\n",
-	   locstart[0],locstart[1],locstart[2],
-	   locend[0],locend[1],locend[2],
-	   totlocpoints);
-  else 
+    /* CCTK_VARIABLE_INT */
+  if (vtype == CCTK_VARIABLE_REAL)
   {
-    printf("proc has no plane\n");
-    printf("LOC START:  %d %d %d   %d %d %d   -> %d\n",
-	   locstart[0],locstart[1],locstart[2],
-	   locend[0],locend[1],locend[2],
-	   totlocpoints);
+    CCTK_INT *int_dptr = (CCTK_REAL *) *loc_dptr;
+
+    for (ipnt[0]=0; ipnt[0]<dirsize[0] && (0 < sdim); ipnt[0]++)
+      for (ipnt[1]=0; ipnt[1]<dirsize[1] && (1 < sdim); ipnt[1]++)
+	for (ipnt[2]=0; ipnt[2]<dirsize[2] && (2 < sdim); ipnt[2]++)
+	{ 
+	  gridpnt[idim] = locstart[idim]; 
+	  for (idim=0;idim<vdim;idim++)
+	    for (idir = 0; idir < sdim; idir ++)
+	      gridpnt[idim] += ipnt[idir] * u_l[idir][idim];
+
+	}
+    int_dptr[lincount++]=((CCTK_INT *) data)
+      [CCTK_GFINDEX3D(GH,gridpnt[0],gridpnt[1],gridpnt[2])];
+    }
   }
 }
    
   
   
   
-  
-  /*$printf("LB/UB [%d , %d, %d]  - [%d, %d, %d] \n",
-	 GH->cctk_lbnd[0],GH->cctk_lbnd[1],GH->cctk_lbnd[2],
-	 GH->cctk_ubnd[0],GH->cctk_ubnd[1],GH->cctk_ubnd[2]);$*/
-  
- 
diff --git a/src/CollectDataND.c b/src/CollectDataND.c
index 0d98886..f70e8e3 100644
--- a/src/CollectDataND.c
+++ b/src/CollectDataND.c
@@ -13,19 +13,26 @@
 #include "CactusBase/IOUtil/src/ioGH.h"
 #include "CactusPUGH/PUGH/src/include/pugh.h"
 
-#define MIN(a,b) ((a)<(b) ? (a) : (b))
-#define MAX(a,b) ((a)>(b) ? (a) : (b))
+
+/* SMART returns the loop boundary: 
+   1 if dir>sdim (the hardcoded index coordinate is less than dimension of slab) 
+   sz[idir] otherwise (size of the slab in direction <dir>.
+
+   This way we make sure that we enter the loop of that particular direction
+   at least one time, even though there is no size associated along that direction.
+*/
+#define SMART(sz,dir,sd) ((dir)>(sd) ? 1 : (sz)[(dir)])
 
 int CollectDataND(cGH *GH, int vindex, int vtimelvl, int vdim, int sdim,
 		  int *S_s, int *dir, int *ds, int *len,
 		  void **dptr, int *dsize)
 {
-  return(0);
+ return(0);
 }
 
 
  /*@@
-   @routine    CollectLocalData2D
+   @routine    CollectLocalDataND
    @date       Fri May 12 17:51:30 2000
    @author     Gerd Lanfermann
    @desc 
@@ -45,116 +52,84 @@ int CollectDataND(cGH *GH, int vindex, int vtimelvl, int vdim, int sdim,
 
 int CollectLocalDataND(cGH *GH, int vindex, int vtimelvl, int vdim, int sdim,
 		       int *S_s, int *linedir, int *ds, int *len,
-		       void **loc_dptr, int *idx1_dptr, int *idx2_dptr, int *dsize)
+		       void **loc_dptr, int *dirsize, int *totsize)
 {
 
-  
-
-  int tstart[MAX_DIM];     /* temp start points */  
-  int tlocstart[MAX_DIM]; 
-  int tlocend[MAX_DIM];
-
-  int *locstart, *locend;          /* local start/end pnt in glob. index coord */
-  int nlocpoints[2], totlocpoints; /* number of loc. points in each dir, tot. pnts */
-
-  hs_vector *u_l;
+  int *locstart, *locend, nlocpoints;
+  int ipnt[MAX_DIM], gridpnt[MAX_DIM];
 
   int vtype,vtypesize;
-  int idim, idir;
-  int ierr;
-  int phasplane=1;
-
+  int idim, idir, lincount;
+  int ierr=0;
+ 
   void *data = CCTK_VarDataPtrI (GH, vtimelvl, vindex);
-  CCTK_REAL *real_dptr = (CCTK_REAL *) *loc_dptr;
+
+  if (vdim > MAX_DIM)
+    return (-1);
 
   vtype      = CCTK_VarTypeI(vindex);
   vtypesize  = CCTK_VarTypeSize(vtype);
 
   locstart = (int*)malloc(vdim*sizeof(int));
   locend   = (int*)malloc(vdim*sizeof(int));
-  u_l      = (hs_vector*)malloc(sdim*sizeof(hs_vector));
+  dirsize  = (int*)malloc(sdim*sizeof(int));
 
-  /* Initialize the local start/end to BAD values and hope for the best... */
-  for (idim=0; idim<vdim; idim++) 
+  printf("Calling GetLocalSlab \n");
+
+  ierr = GetLocalSlab(GH, vdim, sdim, 
+		      S_s, linedir, ds, len, 
+		      locstart, locend, dirsize, &nlocpoints);
+  if (ierr<0)
   {
-    tstart[idim]   = BAD;
-    locstart[idim] = BAD;
-    locend[idim]   = BAD;
-
-    /* line vectors are 000, 000 */
-    u_l[0].vec[idim]  = 0;
-    u_l[1].vec[idim]  = 0;
-
-    /* make sure that all procs above and below the plane, 
-       know that they don't have a plane */
-    if ((idim!=linedir[0]) && (idim!=linedir[1]))
-      if((GH->cctk_lbnd[idim]<S_s[idim]) && 
-	 (GH->cctk_ubnd[idim]>S_s[idim])) 
-      {
-	locstart[idim] = S_s[idim];
-	locend[idim]   = S_s[idim];
-	phasplane = 1;
-      }
-      else phasplane = 0;
+    CCTK_WARN (1, "GetLocalSlab() failed");
+    return (ierr);
   }
 
-  /* construct line vectors that span the surface:
-     e.g. linedir[0,1] = 1,2 means YZ surface:  (000,000) -> (010,001) */
-  u_l[0][linedir[0]]=1;
-  u_l[1][linedir[1]]=1;
+  /* Transform from global to local index coordinate */
+  for(idim=0;idim<vdim;idim++)
+  {  
+    locstart[idim]-=GH->cctk_lbnd[idim];
+    locend[idim]  -=GH->cctk_lbnd[idim];
+  }
 
-  /* init to 1 (neutral element for multiplication) */
-  totlocpoints=1;
-  
-  for (idir=0;idir<2;idir++) 
-  {
+  *loc_dptr= malloc(nlocpoints * vtypesize);
+  lincount=0; 
 
-    if (S_s[linedir[idir]]<GH->cctk_ubnd[linedir[idir]])
-    {
-      /* Initialize startpoints */
-      for (idim=0; idim<vdim; idim++) 
-      {
-	tstart[idim]   = GH->cctk_lbnd[idim]+GH->cctk_nghostzones[idim];
-      }
-      tstart[linedir[idir]]=S_s[linedir[idir]];
-      ierr = GetLocalLine(GH, vdim, tstart, u_l[idir], ds[idir], len[idir], 
-			  tlocstart, tlocend, &nlocpoints[idir]);
-      locstart[linedir[idir]]= tlocstart[linedir[idir]];
-      locend  [linedir[idir]]= tlocend[linedir[idir]];
+  /* zero out index array so that unused dimensions are 0 */
+  memset (gridpnt, 0, sizeof (gridpnt));
+  memset (ipnt   , 0, sizeof (ipnt));
 
-      totlocpoints*=nlocpoints[idir];
+  /* CCTK_VARIABLE_INT */
+  gridpnt[0] = locstart[0]; 
+  gridpnt[1] = locstart[1]; 
+  gridpnt[2] = locstart[2]; 
 
-      if (locstart[linedir[idir]]==BAD) phasplane &=0;
-    }
-    else 
-    {
-      phasplane &=0;
-      totlocpoints = 0;
-    }
-    
-  }
- 
-  if (phasplane)
-    printf("LOC START:  %d %d %d   %d %d %d   -> %d\n",
-	   locstart[0],locstart[1],locstart[2],
-	   locend[0],locend[1],locend[2],
-	   totlocpoints);
-  else 
+  if (vtype == CCTK_VARIABLE_REAL)
   {
-    printf("proc has no plane\n");
-    printf("LOC START:  %d %d %d   %d %d %d   -> %d\n",
-	   locstart[0],locstart[1],locstart[2],
-	   locend[0],locend[1],locend[2],
-	   totlocpoints);
+    CCTK_REAL *real_dptr = (CCTK_REAL *) *loc_dptr;
+    
+    for (ipnt[2]=0 ; ipnt[2]<SMART(dirsize, 2, sdim); ipnt[2]++) 
+
+      gridpnt[linedir[0]] = locstart[linedir[0]] + ipnt[0] * ds[0];
+
+      for (ipnt[1]=0; ipnt[1]<SMART(dirsize, 1, sdim); ipnt[1]++) {
+	
+	gridpnt[linedir[1]] = locstart[linedir[1]] + ipnt[1] * ds[1];
+	
+	for (ipnt[0]=0; ipnt[0]<SMART(dirsize, 0, sdim); ipnt[0]++) {
+	  
+	  { 
+	    gridpnt[linedir[2]] = locstart[linedir[2]] + ipnt[2] * ds[2];    
+
+	    real_dptr[lincount++]=((CCTK_REAL *) data)
+	      [CCTK_GFINDEX3D(GH,gridpnt[0],gridpnt[1],gridpnt[2])];
+	  }
+      }
+    }
+
   }
 }
    
   
   
   
-  
-  /*$printf("LB/UB [%d , %d, %d]  - [%d, %d, %d] \n",
-	 GH->cctk_lbnd[0],GH->cctk_lbnd[1],GH->cctk_lbnd[2],
-	 GH->cctk_ubnd[0],GH->cctk_ubnd[1],GH->cctk_ubnd[2]);$*/
-  
- 
diff --git a/src/GetLocalLine.c b/src/GetLocalLine.c
index f533a45..b105628 100644
--- a/src/GetLocalLine.c
+++ b/src/GetLocalLine.c
@@ -7,14 +7,9 @@
 #include "cctk_Arguments.h"
 
 #include "Hyperslab.h"
+#include "Hyperslabi.h"
 
-#define MAX_DIM 3
-#define MAX_FACE 6
-#define EPS 1e-10
-#define BAD -42
-
-#define ABS(a) ((a)<0 ? -(a) : (a))
-#define SIGN(a) ((a)<0.0 ? (-1.0) : (1.0))
+/*$#define LINE_DEBUG$*/
 
 int IsPointOnProc_INT (int dim, int *gownership, int *pnt);
 int IsPointOnProc_REAL(int dim, int *gownership, CCTK_REAL *pnt);
@@ -51,34 +46,50 @@ int GetLocalLine(cGH *GH,
 		 int vardim, int *S_l, int *u_l, int ds, int len_l, 
 		 int *locstart, int *locend, int *nlocpoints)
 {
+
+  const int n_s[MAX_DIM][MAX_FACE] = 
+  { -1, 1, 0, 0, 0, 0,
+     0, 0,-1, 1, 0, 0,
+     0, 0, 0, 0,-1, 1};
+
+
   /* variables for proc. surfaces */
   CCTK_REAL lmd_s[MAX_FACE], lmd[2], t_lmd;
   int P_s[MAX_DIM][MAX_FACE];
-  int n_s[MAX_DIM][MAX_FACE];
+  int gownership[MAX_FACE];
+
   int E_l[MAX_DIM];
-  int notcoll[MAX_FACE],isloc[2];
-  CCTK_REAL t_pnt1[MAX_DIM], t_pnt2[MAX_DIM];
-  CCTK_REAL secpnt[MAX_DIM][2];
-  int isec;
+  int notcoll[MAX_FACE], isloc[2];
+  CCTK_REAL t_pnt1[MAX_DIM];
 
   int prochasline=0;
 
   int tmp1,tmp2,t;
-  CCTK_REAL rtmp1,rtmp2;
+  CCTK_REAL rtmp1;
 
-  int idim, iface, myproc;
-  int *gownership;
+  int idim, iface, myproc, linsec;
 
-  myproc = CCTK_MyProc(GH);
 
-  gownership= (int*) malloc(2*vardim*sizeof(int));
+  myproc = CCTK_MyProc(GH);
   
-  printf("\n+++++++++\n");
+#ifdef LINE_DEBUG
+  printf("\n+++++++++ GET LOCAL LINE \n");
   printf("S_l [XYZ] : %d %d %d \n",
 	 S_l[0], S_l[1], S_l[2]);
   printf("u_l [XYZ] : %d %d %d    ds %d \n",
 	 u_l[0], u_l[1], u_l[2], ds);
 
+  printf("n_s: face 0/1: [%d %d %d], [%d %d %d] \n",
+	 n_s[0][0],n_s[1][0],n_s[2][0],
+	 n_s[0][1],n_s[1][1],n_s[2][1]);
+  printf("n_s: face 2/3: [%d %d %d], [%d %d %d] \n",
+	 n_s[0][2],n_s[1][2],n_s[2][2],
+	 n_s[0][3],n_s[1][3],n_s[2][3]);
+  printf("n_s: face 4/5: [%d %d %d], [%d %d %d] \n",
+	 n_s[0][4],n_s[1][4],n_s[2][4],
+	 n_s[0][5],n_s[1][5],n_s[2][5]);
+#endif
+
   for (idim=0;idim<vardim;idim++)
   {
     /* GLOBAL OWNERSHIP: lower/upperbnd are INCLUSIVE */
@@ -102,35 +113,12 @@ int GetLocalLine(cGH *GH,
     P_s[idim][5] = gownership[2*idim+1];  /* Aufpunkt, plus faces */
   }
   
-
-  /* norm. vector surface, yz-minus surface (face0) */
-  n_s[0][0]=-1;
-  n_s[1][0]= 0;
-  n_s[2][0]= 0;
-
-  /* norm. vector surface, xz-minus surface (face2) */
-  n_s[0][2]= 0;
-  n_s[1][2]=-1;
-  n_s[2][2]= 0;
-
-  /* norm. vector surface, xy-minus surface (face4) */
-  n_s[0][4]= 0;
-  n_s[1][4]= 0;
-  n_s[2][4]=-1;
-
-  /* plus surfaces (1,3,5) opposite: (-1)*minus_surface */
-  for(idim=0;idim<MAX_DIM;idim++)
-  {
-    n_s[idim][1]=(-1)*n_s[idim][0];
-    n_s[idim][3]=(-1)*n_s[idim][2];
-    n_s[idim][5]=(-1)*n_s[idim][4];
-  }
-
-
+#ifdef LINE_DEBUG
   printf("P_s (face 0): %d %d %d \n", 
 	 P_s[0][0],P_s[1][0],P_s[2][0]);
   printf("P_s (face 1): %d %d %d \n\n",
 	 P_s[0][1],P_s[1][1],P_s[2][1]);
+#endif
 
   /* Lamda for each surface */
   for (iface=0;iface<2*vardim;iface++)
@@ -153,7 +141,9 @@ int GetLocalLine(cGH *GH,
       rtmp1         = (CCTK_REAL)(-1.0)*(tmp1)/(tmp2);
       if (rtmp1<0.0) lmd_s[iface]=BAD;
       else           lmd_s[iface]=rtmp1;
+#ifdef LINE_DEBUG
       printf("lmd_s (face%d) %f \n",iface,rtmp1);
+#endif
     }
   }
 
@@ -163,76 +153,95 @@ int GetLocalLine(cGH *GH,
      it is getting extremly ugly in finding out the local
      points */
 
-  lmd[0]=BAD;
-  lmd[1]=BAD;
-  isec  = 0;
+  lmd[0] = BAD;
+  lmd[1] = BAD;
+  linsec = BAD;
 
+#ifdef LINE_DEBUG
   printf("testing for local Startpt: %d %d %d",
 	 S_l[0],S_l[1],S_l[2]);
+#endif
 
   if (IsPointOnProc_INT(vardim,gownership,S_l)==1)
   {  
+#ifdef LINE_DEBUG
     printf(" ... YES ... valid first\n");
+#endif
     lmd[0]      = 0;
-    isec        = 1;
-    secpnt[0][0]= S_l[0];
-    secpnt[1][0]= S_l[1];
-    secpnt[2][0]= S_l[2];
-
+    linsec      = LIN_INDEX3D(GH,S_l);
   }
   else
+  {
+#ifdef LINE_DEBUG
     printf(" ... NO \n");
+#endif
+  }
   
 
   for (iface=0;iface<MAX_FACE; iface++)
   {
+#ifdef LINE_DEBUG
+    printf("checking face %d: coll: %d  lmd_s %d ", 
+	   iface, notcoll[iface], lmd_s[iface]);
+#endif
     if (notcoll[iface]==1 && lmd_s[iface]!=BAD) 
     {
       for (idim=0;idim<vardim;idim++)
 	t_pnt1[idim]=S_l[idim]+lmd_s[iface]*ds*u_l[idim];
 
-      printf("testing for local intersec. (face %d lamda %+4.2f):  %4.2f %4.2f %4.2f\n",
+#ifdef LINE_DEBUG
+      printf("  testing for local intersec. (face %d lamda %+4.2f):  %4.2f %4.2f %4.2f\n",
 	     iface, lmd_s[iface], t_pnt1[0],t_pnt1[1],t_pnt1[2]);
+#endif
       
       /* check if point is local */
       if (IsPointOnProc_REAL(vardim,gownership,t_pnt1)==1)
       {
+#ifdef LINE_DEBUG
 	printf("    ... YES ... ");
+#endif
 
-	/* first valid face: this point is good */
-	if (isec==0)
+	/* CASE1: first valid face: this point is good */
+	if (linsec<0)
 	{
-	  printf(" valid first %d\n",isec);
-	  lmd[0]      =lmd_s[iface];
-	  secpnt[0][0]=t_pnt1[0];
-	  secpnt[1][0]=t_pnt1[1];
-	  secpnt[2][0]=t_pnt1[2];
-	  isec++;
+#ifdef LINE_DEBUG
+	  printf(" valid first \n");
+#endif
+	  lmd[0]      = lmd_s[iface];
+	  linsec      = LIN_INDEX3D(GH,t_pnt1);
 	}            
-	/* we already have intersec. check for duplicates (make sure 
-	   that we allow the endpoint to be located on the startpoint lmd=0 )*/
-	else if ((secpnt[0][0]==t_pnt1[0]) &&
-		 (secpnt[1][0]==t_pnt1[1]) &&
-		 (secpnt[2][0]==t_pnt1[2]))
-	  {
-	    printf(" but duplicate, skip \n");
-	  }
+	/* CASE2: we already have intersec. 
+	   check for duplicates: eg. a diagonal through the n-space, with 
+	   one exception: make sure that we allow the endpoint to be located 
+	   on the startpoint if lmd == 0 (which is the case if we start at the last
+	   point of a proc) */
+	else if (linsec==LIN_INDEX3D(GH,t_pnt1) && lmd_s[iface]>EPS)
+	{
+#ifdef LINE_DEBUG
+	  printf(" but duplicate, skip \n");
+#endif
+	}
 	else
-	  {
-	    printf(" valid second %d\n",isec);
-	    lmd[1]      =lmd_s[iface];
-	    secpnt[0][1]=t_pnt1[0];
-	    secpnt[1][1]=t_pnt1[1];
-	    secpnt[2][1]=t_pnt1[2]; 
-	    isec++;
-	  }
+	{
+#ifdef LINE_DEBUG
+	  printf(" valid second \n");
+#endif
+	  lmd[1]      =lmd_s[iface];
+	}
       }
       else
       {
+#ifdef LINE_DEBUG
 	printf("     ... NO \n");
+#endif
       }
     }
-    else printf(" ... skipping face: %d %4.2f \n",iface,lmd_s[iface]);
+    else
+    {
+#ifdef LINE_DEBUG
+    printf(" ... skipping face: %d %4.2f \n",iface,lmd_s[iface]);
+#endif
+    }
   }
  
   /* flag if a proc has intersections, order the intersections 
@@ -247,30 +256,28 @@ int GetLocalLine(cGH *GH,
       t_lmd  = lmd[0];
       lmd[0] = lmd[1];
       lmd[1] = t_lmd;
-      for (idim=0;idim<vardim;idim++)
-      {
-	t_pnt1[idim]   = secpnt[idim][0];
-	secpnt[idim][0]= secpnt[idim][1];
-	secpnt[idim][1]= t_pnt1[idim];
-      }
     }
     lmd[0]=ceil (lmd[0]);
     lmd[1]=floor(lmd[1]);
+#ifdef LINE_DEBUG
     printf("LOCAL INTERSECTIONS: (ordered in lmd) %4.2f %4.2f \n",lmd[0],lmd[1]);
+#endif
   }
   else if (lmd[0]==BAD && lmd[1]==BAD)
   {   
     prochasline = 0;
+#ifdef LINE_DEBUG
     printf("NO INTERSECTION ON THIS PROC\n");
+#endif
   }
   else if ((lmd[0]==BAD && lmd[1]!=BAD) ||
 	   (lmd[1]==BAD && lmd[0]!=BAD))
   {
-    printf("ERROR: only one local intersection: that is wrong %4.2f %4.2f\n",
+    printf("ERROR: only one local intersection: you are screwed. lmd 0/1  %4.2f %4.2f\n",
 	   lmd[0],lmd[1]);
   }
   
-  /* patch if lambdas are greater line length */
+  /* patch if lambdas are greater specified line length */
   if (len_l>=0) 
   {
     if (len_l>lmd[0] && len_l<lmd[1])
@@ -283,7 +290,7 @@ int GetLocalLine(cGH *GH,
     }
   }
 
-  /* order intersections in lambda-distance from aufpunkt */  
+  /* Calculate the intersections points in global index space*/
   if (prochasline) 
   {
     nlocpoints[0]= ABS(lmd[0]-lmd[1])+1;
@@ -303,19 +310,12 @@ int GetLocalLine(cGH *GH,
     }
   }
 
+#ifdef LINE_DEBUG
   printf("\n>>>>> GLOBAL START/END: (%d %d %d) -> (%d %d %d)  = %d\n",
 	 locstart[0],locstart[1],locstart[2],
 	 locend[0],locend[1],locend[2],nlocpoints[0]);
+#endif
 
-  /*$for (idim=0;idim<vardim;idim++)
-  {
-    locstart[idim] -= GH->cctk_lbnd[idim];
-    locend[idim]   -= GH->cctk_lbnd[idim];
-  }
-  printf(">>>>> LOCAL START/END : (%d %d %d) -> (%d %d %d)  = %d\n",
-	 locstart[0],locstart[1],locstart[2],
-	 locend[0],locend[1],locend[2],nlocpoints[0]);$*/
-	
   return(0);
   
 }
diff --git a/src/Hyperslab.h b/src/Hyperslab.h
index e0be51f..7376341 100644
--- a/src/Hyperslab.h
+++ b/src/Hyperslab.h
@@ -16,8 +16,19 @@ int CollectLocalData2D(cGH *GH, int vindex, int vtimelvl, int vdim,
 		       int *S_s, int dir[2], int ds[2], int len[2],
 		       void **dptr, int *idx1_ptr, int *idx2_ptr, int *dsize);
 
+int CollectDataND(cGH *GH, int vindex, int vtimelvl, int vdim, int sdim,
+		  int *S_s, int *dir, int *ds, int *len,
+		  void **dptr, int *dsize);
+
+int CollectLocalDataND(cGH *GH, int vindex, int vtimelvl, int vdim, int sdim,
+		       int *S_s, int *dir, int *ds, int *len,
+		       void **dptr,  int *dirsize,int *totsize);
 int GetLocalLine(cGH *GH, 
 		 int vardim, int *S_l, int *u_l, int ds, int len_l, 
 		 int *locstart, int *locend, int *nlocpoints);
 
+int GetLocalSlab(cGH *GH, int vdim, int sdim,
+		 int *S_s, int *linedir, int *ds, int *len,
+		 int *locstart, int *locend, int *dirsize, int *totsize);
+
 int IsOnGlobGrid(cGH *GH, int vdim, int *grdpnt);
diff --git a/src/Hyperslabi.h b/src/Hyperslabi.h
index 9ddfa55..2adae72 100644
--- a/src/Hyperslabi.h
+++ b/src/Hyperslabi.h
@@ -1,9 +1,15 @@
 
 
 #define MAX_DIM 3
+#define MAX_FACE 6
 #define BAD -42
+#define EPS 1e-10
+
+#define ABS(a) ((a)<0 ? -(a) : (a))
+#define SIGN(a) ((a)<0.0 ? (-1.0) : (1.0))
+#define LIN_INDEX3D(GH,P) ((P[0])+ (GH)->cctk_gsh[0]*((P[1]) + (GH)->cctk_gsh[1]*(P[2])))
+
+#define MIN(a,b) ((a)<(b) ? (a) : (b))
+#define MAX(a,b) ((a)>(b) ? (a) : (b))
 
 
-typedef struct t_vector {
-  int vec[MAX_DIM];
-}  hs_vector;
diff --git a/src/TestSlab.c b/src/TestSlab.c
index ac6e72f..9f48e3f 100644
--- a/src/TestSlab.c
+++ b/src/TestSlab.c
@@ -58,15 +58,6 @@ void TestSlab(CCTK_ARGUMENTS) {
   vindex   = CCTK_VarIndex("grid::z");
   vtimelvl = CCTK_NumTimeLevelsFromVarI (vindex) - 1;
   vdim     = CCTK_GroupDimFromVarI(vindex);
-  
-
-  /*$gs[0]=1;
-  gs[1]=5;
-  gs[2]=9;
- 
-  di[0]= 0;
-  di[1]=-1;
-  di[2]=-2;$*/
 
   gs[0]=1;
   gs[1]=1;
diff --git a/src/TestSlab2D.c b/src/TestSlab2D.c
index 48d6b86..f8cf6fa 100644
--- a/src/TestSlab2D.c
+++ b/src/TestSlab2D.c
@@ -1,4 +1,5 @@
 
+
 #include <stdio.h>
 #include <stdlib.h>
 #include <math.h>  
@@ -22,10 +23,11 @@
 void TestSlab2D(CCTK_ARGUMENTS) {  
   DECLARE_CCTK_ARGUMENTS
 
-    int gs[3], di[2];    /* global start/direction (2times) */
+  int gs[3], di[2];    /* global start/direction (2times) */
   int ds[2], len[2];           /* downsampling/length */
   int ls[3], le[3];            /* local start/end */
-  int dsize[1];                /* number of points */
+  int dirsize[2],  totsize;      /* number of points (dir./tot)*/
+  int locstart[3], locend[3];
   int nprocs,iproc;            /* number of procs */
 
   int l;
@@ -41,20 +43,25 @@ void TestSlab2D(CCTK_ARGUMENTS) {
 
   gs[0]=6;
   gs[1]=6;
-  gs[2]=5;
+  gs[2]=3;
  
   di[0]= 2;
   di[1]= 1;
 
-  ds[0]=2;
-  ds[1]=2;
+  ds[0]=1;
+  ds[1]=1;
 
   len[0] = -1;
   len[1] = -1;
 
-  printf("Calling CollectData2D \n");
-  CollectLocalData2D(cctkGH, vindex, vtimelvl, vdim, 
-		     gs, di, ds, len,
-		     &dptr, index1, index2, dsize);
+  
+
+  printf("Calling CollectLocalDataND \n");
+  CollectLocalDataND(cctkGH, vindex, vtimelvl, vdim, 2,
+		     gs, di, ds, len, 
+		     &dptr, dirsize, &totsize);
 
 }
+  /*$GetLocalSlab(cctkGH, vdim, 2,
+	       gs, di, ds, len,
+	       locstart, locend, dirsize, &totsize);$*/
diff --git a/src/make.code.defn b/src/make.code.defn
index a86191c..b448166 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -2,7 +2,7 @@
 # $Header$
 
 # Source files in this directory
-SRCS = GetLocalLine.c TestSlab.c TestSlab2D.c CollectData1D.c CollectData2D.c Misc.c
+SRCS = GetLocalLine.c GetLocalSlab.c TestSlab.c TestSlab2D.c CollectData1D.c Misc.c CollectDataND.c
 
 # Subdirectories containing source files
 SUBDIRS =