Removed old hyperslabbing code for diagonals which was working only for

cubic grids. The new code (which is already available via the new
hyperslab API) can do diagonals from non-cubic variables also (but still
only for non-staggered 3D variables).


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

6 files changed, 255 insertions, 669 deletions

diff --git a/src/CollectData1D.c b/src/CollectData1D.c
deleted file mode 100644
index 38b5f60..0000000
--- a/src/CollectData1D.c
+++ /dev/null
@@ -1,240 +0,0 @@
-
-#include <stdio.h>
-#include <stdlib.h>
-
-#include "cctk.h"
-#include "CactusPUGH/PUGH/src/include/pugh.h"
-#include "PUGHSlab.h"
-
-static const char *rcsid = "$Header$";
-
-CCTK_FILEVERSION(CactusPUGH_PUGHSlab_CollectData1D_c)
-
-/* #define DEBUGME  1 */
-
-
-#define PICKUP_LINE_DATA(GH, cctk_type, vsize, start, directions, vdata, hdata)\
-          {                                                                   \
-            int i, dim, point [3];                                            \
-            cctk_type *typed_data = (cctk_type *) hdata;                      \
-                                                                              \
-                                                                              \
-            for (i = 0; i < vsize; i++)                                       \
-            {                                                                 \
-              for (dim = 0; dim < 3; dim++)                                   \
-                point [dim] = start [dim] + i*directions [dim];               \
-              *typed_data++ = ((cctk_type *) vdata)                           \
-                [CCTK_GFINDEX3D (GH, point [0], point [1], point [2])];       \
-            }                                                                 \
-          }
-
-/* local function prototypes */
-int Hyperslab_CollectData1D(const cGH *GH, int vindex, int vtimelvl,
-                            const int *origin, const int *directions,
-                            int downsample, int length,
-                            void **hdata, int *hsize, int proc);
-int GetLocalLine(const cGH *GH, pGExtras *extras,
-                 int vardim, const int *S_l, const int *u_l, int ds, int len_l,
-                 int *locstart, int *locend, int *nlocpoints);
-
-
-static int CollectLocalData1D(const cGH *GH, int vindex, int vtimelvl,
-                              const int *origin, const int *directions,
-                              int downsample, int length,
-                              void **hdata, int *hsize)
-{
-  int idim, ierr;
-  cGroup groupinfo;
-  int locstart [3], locend [3];
-  void *vdata = CCTK_VarDataPtrI (GH, vtimelvl, vindex);
-  pGExtras *extras = ((pGA *) PUGH_pGH (GH)->variables [vindex][vtimelvl])->extras;
-
-
-
-  CCTK_GroupData (CCTK_GroupIndexFromVarI (vindex), &groupinfo);
-
-  if (groupinfo.dim != 3)
-  {
-    CCTK_WARN (1, "CollectLocalData1D() processes 3D variables only");
-    return (-1);
-  }
-
-  ierr = GetLocalLine(GH, extras, 3, origin, directions, downsample, length, 
-                          locstart, locend, hsize);
-  if (ierr < 0)
-  {
-    CCTK_WARN (1, "GetLocalLine() failed");
-    return (ierr);
-  }
-
-  /* Transform from global to local index coordinate */
-  for(idim=0;idim<3;idim++)
-  {  
-    locstart[idim]-=GH->cctk_lbnd[idim];
-    locend[idim]  -=GH->cctk_lbnd[idim];
-  }
-
-#ifdef DEBUGME
-  printf("CollectLocal1D: %d npoints   lb: %d %d %d\n",
-         *hsize,GH->cctk_lbnd[0],GH->cctk_lbnd[1],GH->cctk_lbnd[2]);
-#endif
-
-  *hdata= malloc (*hsize * CCTK_VarTypeSize (groupinfo.vartype));
-
-  /* CCTK_VARIABLE_INT */
-  switch (groupinfo.vartype)
-  {
-    case CCTK_VARIABLE_CHAR:
-      PICKUP_LINE_DATA(GH, CCTK_BYTE, *hsize, locstart, directions, vdata, *hdata);
-      break;
-
-    case CCTK_VARIABLE_INT:
-      PICKUP_LINE_DATA(GH, CCTK_INT, *hsize, locstart, directions, vdata, *hdata);
-      break;
-
-    case CCTK_VARIABLE_REAL:
-      PICKUP_LINE_DATA(GH, CCTK_REAL, *hsize, locstart, directions, vdata, *hdata);
-      break;
-
-    case CCTK_VARIABLE_COMPLEX:
-      PICKUP_LINE_DATA(GH, CCTK_COMPLEX, *hsize, locstart, directions, vdata, *hdata);
-      break;
-
-    default:
-      CCTK_WARN (1, "Unsupported variable type");
-      ierr = -1;
-      break;
-  }
-
-  return (ierr);
-}
-
-
-int Hyperslab_CollectData1D(const cGH *GH, int vindex, int vtimelvl,
-                            const int *origin, const int *directions,
-                            int downsample, int length,
-                            void **hdata, int *hsize, int proc)
-{
-
-  void *loc_dptr;
-  int   loc_dcount;
-  void *all_dptr;
-  int   all_dcount=0;
-  int vtype, vtypesize;
-  int myproc,nprocs;
-  int ierr;
-
-#ifdef CCTK_MPI
-  int iproc;
-  const pGH *pughGH;   
-  CCTK_INT  tmp, *rem_dcount;
-  MPI_Datatype vmpitype;
-  int *recvcnts, *displs;
-#endif
-
-
-  vtype    = CCTK_VarTypeI(vindex);
-  vtypesize= CCTK_VarTypeSize(vtype);
-
-  myproc   = CCTK_MyProc(GH); 
-  nprocs   = CCTK_nProcs(GH);
-
-
-  ierr = CollectLocalData1D(GH, vindex, vtimelvl,
-                            origin, directions, downsample, length, 
-                            &loc_dptr, &loc_dcount);
-
-  if (ierr)
-  {
-    CCTK_WARN (1, "CollectLocalData1D() failed");
-    return (ierr);
-  }
-
-#ifdef CCTK_MPI
-
-  pughGH   = PUGH_pGH(GH);
-
-  rem_dcount = NULL;
-  recvcnts = NULL;
-  displs = NULL;
-  vmpitype = PUGH_MPIDataType (pughGH, vtype);
-
-  /* Send local number of data elements */
-  tmp = (CCTK_INT) loc_dcount;
-
-  if (proc < 0 || proc == myproc)
-  {
-    rem_dcount= (CCTK_INT *) malloc (nprocs * sizeof (CCTK_INT));
-    recvcnts  = (int *) malloc (nprocs * sizeof (int));
-    displs    = (int *) malloc (nprocs * sizeof (int));
-  }
-
-  if (proc < 0)
-  {
-    CACTUS_MPI_ERROR (MPI_Allgather (&tmp, 1, PUGH_MPI_INT,
-                                     rem_dcount, 1, PUGH_MPI_INT,
-                                     pughGH->PUGH_COMM_WORLD));
-  } else {
-    CACTUS_MPI_ERROR (MPI_Gather (&tmp, 1, PUGH_MPI_INT,
-                                  rem_dcount, 1, PUGH_MPI_INT,
-                                  proc, pughGH->PUGH_COMM_WORLD));
-  }
-
-  if (proc < 0 || proc == myproc)
-  {
-    displs[0]  = 0;
-    all_dcount = recvcnts[0] = (int) rem_dcount[0];
-    for (iproc=1;iproc<nprocs;iproc++)
-    {
-      recvcnts[iproc] = (int) rem_dcount[iproc];
-      displs[iproc] = displs[iproc-1] + recvcnts[iproc-1];
-      all_dcount   += recvcnts[iproc];
-    }
-
-#ifdef DEBUGME
-    printf ("  local hsize: %d\n", loc_dcount);
-    for (iproc=0;iproc<nprocs;iproc++)
-      printf(" DISPLS[%d] rem: %d   displs: %d   all %d \n",
-             iproc, recvcnts[iproc],displs[iproc],all_dcount);
-    fflush(stdout);
-#endif
-
-    all_dptr = malloc (all_dcount * vtypesize);
-  }
-  else
-  {
-    all_dptr = NULL;
-  }
-
-  if (proc < 0)
-  {
-    CACTUS_MPI_ERROR (MPI_Allgatherv(loc_dptr, loc_dcount, vmpitype, 
-                                     all_dptr, recvcnts, displs, vmpitype,
-                                     pughGH->PUGH_COMM_WORLD));
-  } else {
-    CACTUS_MPI_ERROR (MPI_Gatherv(loc_dptr, loc_dcount, vmpitype, 
-                                  all_dptr, recvcnts, displs, vmpitype,
-                                  proc, pughGH->PUGH_COMM_WORLD));
-  }
-
-  if (proc < 0 || proc == myproc)
-  {
-    free (rem_dcount);
-    free (recvcnts);
-    free (displs);
-  }
-
-#else
-  all_dptr = loc_dptr;
-  all_dcount = loc_dcount;
-#endif
-
-  /* assign return values */
-  if (proc < 0 || proc == myproc)
-  {
-    *hdata = all_dptr;
-    *hsize = all_dcount;
-  }
-
-  return (0);
-}
diff --git a/src/GetHyperslab.c b/src/GetHyperslab.c
index b6e0af0..b1206e2 100644
--- a/src/GetHyperslab.c
+++ b/src/GetHyperslab.c
@@ -39,19 +39,12 @@ static int GetLocalHyperslab (const cGH *GH,
                               int timelevel,
                               int hdatatype,
                               void *hdata);
-/********************************************************************
- ********************    External Routines   ************************
- ********************************************************************/
-/* Gerd's routine to get 1D lines from a 3D array */
-int Hyperslab_CollectData1D (const cGH *GH, int vindex, int vtimelvl,
-                             const int *origin,
-                             const int *directions,
-                             int downsampling,
-                             int length,
-                             void **hdata,
-                             int *hsize,
-                             int proc);
-
+static int GetDiagonalFromFrom3D (const cGH *GH,
+                                  const hslab_mapping_t *mapping,
+                                  int vindex,
+                                  int timelevel,
+                                  int hdatatype,
+                                  void *hdata);
 
 
 CCTK_INT Hyperslab_Get (const cGH *GH,
@@ -274,27 +267,11 @@ static int GetLocalHyperslab (const cGH *GH,
     return (0);
   }
 
-  /* FIXME: hack for getting diagonals from 3D variables
-            This is calling Gerd's CollectData1D() routine which can
-            extract non-axis-parallel lines too but is fixed to 3D data. */
+  /* diagonals from 3D variables are treated special */
   if (mapping->is_diagonal_in_3D)
   {
-    const int origin[] = {0, 0, 0};
-    const int directions[] = {1, 1, 1};
-    int dummy_hsize;
-    void *dummy_hdata;
-
-
-    retval = Hyperslab_CollectData1D (GH, vindex, timelevel, origin, directions,
-                                      1, -1, &dummy_hdata, &dummy_hsize,
-                                      mapping->target_proc);
-    if (! retval)
-    {
-      memcpy (hdata, dummy_hdata,
-              mapping->totals * CCTK_VarTypeSize (vinfo.vartype));
-      free (dummy_hdata);
-    }
-
+    retval = GetDiagonalFromFrom3D (GH, mapping, vindex, timelevel, hdatatype,
+                                    hdata);
     return (retval);
   }
 
@@ -427,7 +404,7 @@ static int GetLocalHyperslab (const cGH *GH,
 
       /* get the byte pointer into the source array */
       typed_vdata = (const char *) vdata + point[0];
-#if 1
+#if 0
 fprintf (stderr, "***** base vdata %p offset %d '%s'\n", vdata, point[0], CCTK_FullName (vindex));
 #endif
       for (i = 1; i < vinfo.dim; i++)
@@ -445,7 +422,7 @@ fprintf (stderr, "***** base vdata %p offset %d '%s'\n", vdata, point[0], CCTK_F
         }
         else
         {
-#if 1
+#if 0
 fprintf (stderr, "***** copying %d bytes from %p tp %p\n", dim0_hsize, typed_vdata, typed_hdata);
 #endif
           memcpy (typed_hdata, typed_vdata, dim0_hsize);
@@ -489,3 +466,156 @@ fprintf (stderr, "***** copying %d bytes from %p tp %p\n", dim0_hsize, typed_vda
 
   return (0);
 }
+
+
+static int GetDiagonalFromFrom3D (const cGH *GH,
+                                  const hslab_mapping_t *mapping,
+                                  int vindex,
+                                  int timelevel,
+                                  int hdatatype,
+                                  void *hdata)
+{
+  int i, j, k, myproc, nprocs, linear_idx;
+  CCTK_INT local_npoints, npoints;
+  int vdatatype, htypesize, vtypesize;
+  const char *vdata;
+  char *local_hdata;
+  const pGH *pughGH;
+  const pGExtras *extras;
+#ifdef CCTK_MPI
+  CCTK_INT *global_npoints;
+  int *recvcnts, *displs;
+  MPI_Datatype mpidatatype;
+#endif
+
+
+
+  /* get the pGH pointer and the variable's GA structure */
+  pughGH = PUGH_pGH (GH);
+  extras = ((const pGA *) pughGH->variables[vindex][timelevel])->extras;
+
+  vdatatype = CCTK_VarTypeI (vindex);
+  if (hdatatype < 0)
+  {
+    hdatatype = vdatatype;
+  }
+  htypesize = CCTK_VarTypeSize (hdatatype);
+  vtypesize = CCTK_VarTypeSize (vdatatype);
+  vdata = (const char *) CCTK_VarDataPtrI (GH, timelevel, vindex);
+
+  myproc = CCTK_MyProc (GH);
+  nprocs = CCTK_nProcs (GH);
+  if (nprocs == 1)
+  {
+    local_hdata = (char *) hdata;
+  }
+  else
+  {
+    local_hdata = (char *) malloc (mapping->global_hsize[0] * htypesize);
+  }
+
+  i = mapping->global_startpoint[0] - extras->lb[myproc][0];
+  j = mapping->global_startpoint[1] - extras->lb[myproc][1];
+  k = mapping->global_startpoint[2] - extras->lb[myproc][2];
+  local_npoints = 0;
+  for (npoints = 0; npoints < mapping->global_hsize[0]; npoints++)
+  {
+    if (i >= extras->ownership[0][0][0] && i < extras->ownership[0][1][0] &&
+        j >= extras->ownership[0][0][1] && j < extras->ownership[0][1][1] &&
+        k >= extras->ownership[0][0][2] && k < extras->ownership[0][1][2])
+    {
+      linear_idx = i + j*extras->hyper_volume[1] + k*extras->hyper_volume[2];
+      if (vdatatype != hdatatype)
+      {
+        mapping->conversion_fn (vdata + linear_idx*vtypesize, local_hdata,
+                                1, 1, 1);
+      }
+      else
+      {
+        memcpy (local_hdata, vdata + linear_idx*vtypesize, htypesize);
+      }
+      local_hdata += htypesize;
+      local_npoints++;
+    }
+    i += mapping->downsample[0];
+    j += mapping->downsample[1];
+    k += mapping->downsample[2];
+  }
+  local_hdata -= local_npoints * htypesize;
+
+#ifdef CCTK_MPI
+  if (nprocs > 1)
+  {
+    /* allocate communication buffers */
+    myproc = CCTK_MyProc (GH);
+    if (mapping->target_proc < 0 || mapping->target_proc == myproc)
+    {
+      global_npoints = (CCTK_INT *) malloc (nprocs * sizeof (CCTK_INT));
+      recvcnts = (int *) malloc (2 * nprocs * sizeof (int));
+      displs   = recvcnts + nprocs;
+    }
+    else
+    {
+      global_npoints = NULL; recvcnts = displs = NULL; hdata = NULL;
+    }
+
+    /* gather the number of local points on each processor */
+    if (mapping->target_proc < 0)
+    {
+      CACTUS_MPI_ERROR (MPI_Allgather (&local_npoints, 1, PUGH_MPI_INT,
+                                       global_npoints, 1, PUGH_MPI_INT,
+                                       pughGH->PUGH_COMM_WORLD));
+    }
+    else
+    {
+      CACTUS_MPI_ERROR (MPI_Gather (&local_npoints, 1, PUGH_MPI_INT,
+                                    global_npoints, 1, PUGH_MPI_INT,
+                                    mapping->target_proc,
+                                    pughGH->PUGH_COMM_WORLD));
+    }
+
+    /* compute the receive count and displacement vectors */
+    if (mapping->target_proc < 0 || mapping->target_proc == myproc)
+    {
+      for (i = 0; i < nprocs; i++)
+      {
+        recvcnts[i] = (int) global_npoints[i];
+        displs[i] = i == 0 ? 0 : displs[i-1] + recvcnts[i-1];
+      }
+    }
+
+    /* get the MPI datatype for the hyperslab data */
+    mpidatatype = PUGH_MPIDataType (pughGH, hdatatype);
+
+    /* gather the local hyperslab data from each processor */
+    if (mapping->target_proc < 0)
+    {
+      CACTUS_MPI_ERROR (MPI_Allgatherv (local_hdata, local_npoints,
+                                        mpidatatype, hdata, recvcnts, displs,
+                                        mpidatatype, pughGH->PUGH_COMM_WORLD));
+    }
+    else
+    {
+      CACTUS_MPI_ERROR (MPI_Gatherv (local_hdata, local_npoints,
+                                     mpidatatype, hdata, recvcnts, displs,
+                                     mpidatatype, mapping->target_proc,
+                                     pughGH->PUGH_COMM_WORLD));
+    }
+
+    /* free communication buffers */
+    if (mapping->target_proc < 0 || mapping->target_proc == myproc)
+    {
+      free (global_npoints);
+      free (recvcnts);
+    }
+  }
+#endif  /* CCTK_MPI */
+
+  /* free allocated resources */
+  if (nprocs > 1)
+  {
+    free (local_hdata);
+  }
+
+  return (0);
+}
diff --git a/src/GetLocalLine.c b/src/GetLocalLine.c
deleted file mode 100644
index 78a9b57..0000000
--- a/src/GetLocalLine.c
+++ /dev/null
@@ -1,334 +0,0 @@
-#include <stdio.h>
-#include <stdlib.h>
-#include <math.h>  
-
-#include "cctk.h"
-#include "CactusPUGH/PUGH/src/include/pugh.h"
-
-static const char *rcsid = "$Header$";
-
-CCTK_FILEVERSION(CactusPUGH_PUGHSlab_GetLocalLine_c)
-
-#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(extras,P) ((P[0])+ (extras)->nsize[0]*((P[1]) + (extras)->nsize[1]*(P[2])))
-
-#define MIN(a,b) ((a)<(b) ? (a) : (b))
-#define MAX(a,b) ((a)>(b) ? (a) : (b))
-
-/*$#define LINE_DEBUG$*/
-
-int IsPointOnProc_INT (int dim, int *gownership, const int *pnt);
-int IsPointOnProc_REAL(int dim, int *gownership, CCTK_REAL *pnt);
-
-int IsPointOnProc_INT(int dim, int *gownership, const int *pnt)
-{
-  int retval=1; /*assume YES */
-  int idim;
-  for (idim=0;idim<dim;idim++)
-  {
-    /* out of range -> set retval to zero;*/
-    if ((gownership[2*idim  ]>pnt[idim]) || 
-	(gownership[2*idim+1]<pnt[idim]))
-      retval&=0;  /* incremental AND on the return value */
-  }
-  return(retval);
-}
-
-int IsPointOnProc_REAL(int dim, int *gownership, CCTK_REAL *pnt)
-{
-  int retval=1; /*assume YES */
-  int idim;
-  for (idim=0;idim<dim;idim++)
-  {
-    /* out of range -> set retval to zero;*/
-    if ((gownership[2*idim  ]>pnt[idim]) || 
-	(gownership[2*idim+1]<pnt[idim]))
-      retval&=0;  /* incremental AND on the return value */
-  }
-  return(retval);
-}
-
-int GetLocalLine(cGH *GH, pGExtras *extras,
-		 int vardim, const int *S_l, const 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 gownership[MAX_FACE];
-
-  int notcoll[MAX_FACE];
-  CCTK_REAL t_pnt1[MAX_DIM];
-
-  int prochasline=0;
-
-  int tmp1,tmp2;
-  CCTK_REAL rtmp1;
-
-  int idim, iface, linsec;
-  int myproc;
-
-
-
-#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
-
-  myproc = CCTK_MyProc (GH);
-
-  for (idim=0;idim<vardim;idim++)
-  {
-    /* GLOBAL OWNERSHIP: lower/upperbnd are INCLUSIVE */
-    gownership[2*idim]  = extras->lb[myproc][idim];
-    if (GH->cctk_bbox[2*idim]==0) 
-      gownership[2*idim]  += extras->nghostzones[idim];
-
-    gownership[2*idim+1]= extras->ub[myproc][idim];
-    if (GH->cctk_bbox[2*idim+1]==0) 
-      gownership[2*idim+1]-= (extras->nghostzones[idim]);    
-  }
-
-  
-  for (idim=0;idim<vardim;idim++)
-  {
-    P_s[idim][0] = gownership[2*idim];    /* Aufpunkt, minus faces */
-    P_s[idim][2] = gownership[2*idim];    /* Aufpunkt, minus faces */
-    P_s[idim][4] = gownership[2*idim];    /* Aufpunkt, minus faces */
-    P_s[idim][1] = gownership[2*idim+1];  /* Aufpunkt, plus faces */
-    P_s[idim][3] = gownership[2*idim+1];  /* Aufpunkt, plus faces */
-    P_s[idim][5] = gownership[2*idim+1];  /* Aufpunkt, plus faces */
-  }
-  
-#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++)
-  {
-    tmp1=0.0; tmp2=0.0;
-
-    for (idim=0;idim<vardim;idim++) 
-    {
-      tmp1+= (S_l[idim]-P_s[idim][iface])*n_s[idim][iface];
-      tmp2+= (ds*u_l[idim]*n_s[idim][iface]);
-    }
-    if (tmp2==0) 
-    {
-      lmd_s[iface]=BAD;
-      notcoll[iface]=0;
-    }
-    else
-    {
-      notcoll[iface]=1;
-      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
-    }
-  }
-
-   /* Find 2 valid intersection with the proc surfaces */
-  /* tricky: duplicate points corners and the startpoint */
-  /* We have to restrict the line to positive lambda, otherwise
-     it is getting extremly ugly in finding out the local
-     points */
-
-  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;
-    linsec      = LIN_INDEX3D(extras,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];
-
-#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
-
-	/* CASE1: first valid face: this point is good */
-	if (linsec<0)
-	{
-#ifdef LINE_DEBUG
-	  printf(" valid first \n");
-#endif
-	  lmd[0]      = lmd_s[iface];
-	  linsec      = LIN_INDEX3D(extras,t_pnt1);
-	}            
-	/* 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(extras,t_pnt1) && lmd_s[iface]>EPS)
-	{
-#ifdef LINE_DEBUG
-	  printf(" but duplicate, skip \n");
-#endif
-	}
-	else
-	{
-#ifdef LINE_DEBUG
-	  printf(" valid second \n");
-#endif
-	  lmd[1]      =lmd_s[iface];
-	}
-      }
-      else
-      {
-#ifdef LINE_DEBUG
-	printf("     ... NO \n");
-#endif
-      }
-    }
-    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 
-     if we have any  */
-  if (lmd[0]!=BAD && lmd[1]!=BAD)
-  {   
-    prochasline=1;
-
-    /* switch: lmd[0] has to be the closer one to S_l */
-    if (lmd[0]>lmd[1])
-    {
-      t_lmd  = lmd[0];
-      lmd[0] = lmd[1];
-      lmd[1] = t_lmd;
-    }
-    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: you are screwed. lmd 0/1  %4.2f %4.2f\n",
-	   lmd[0],lmd[1]);
-  }
-  
-  /* patch if lambdas are greater specified line length */
-  if (len_l>=0) 
-  {
-    if (len_l>lmd[0] && len_l<lmd[1])
-      lmd[1]=floor(len_l);
-    else if (len_l<lmd[0])
-    {
-      lmd[0]=BAD;
-      lmd[1]=BAD;
-      prochasline=0;
-    }
-  }
-
-  /* Calculate the intersections points in global index space*/
-  if (prochasline) 
-  {
-    nlocpoints[0]= ABS(lmd[0]-lmd[1])+1;
-    for (idim=0;idim<vardim;idim++)
-    {
-      locstart[idim] = (S_l[idim]+lmd[0]*ds*u_l[idim]);
-      locend[idim]   = (S_l[idim]+lmd[1]*ds*u_l[idim]);
-    }
-  } 
-  else
-  {
-    nlocpoints[0]  = 0;
-    for (idim=0;idim<vardim;idim++)
-    {
-      locstart[idim] = BAD;
-      locend[idim]   = BAD;  
-    }
-  }
-
-#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
-
-  return(0);
-  
-}
diff --git a/src/Hyperslab.c b/src/Hyperslab.c
index 9c482e6..dd5f840 100644
--- a/src/Hyperslab.c
+++ b/src/Hyperslab.c
@@ -131,16 +131,6 @@ CCTK_FILEVERSION(CactusPUGH_PUGHSlab_Hyperslab_c)
     }
 
 
-/*** local function prototypes ***/
-/* Gerd's routine to get 1D lines from a 3D array */
-int Hyperslab_CollectData1D (const cGH *GH, int vindex, int vtimelvl,
-                             const int *origin,
-                             const int *directions,
-                             int downsampling,
-                             int length,
-                             void **hdata,
-                             int *hsize,
-                             int proc);
 /* routine to check parameters passed to the Hyperslab routines */
 static const char *checkParameters (const cGH *GH, int vindex, int vtimelvl,
                                     int hdim,
@@ -746,46 +736,6 @@ int Hyperslab_GetHyperslab (const cGH *GH,
 #endif
 
 
-  /* do some plausibility checks */
-  errormsg = checkParameters (GH, vindex, vtimelvl, hdim,
-                              global_origin, directions, extents,
-                              downsample_, hdata, hsize);
-
-  /* FIXME: hack for getting diagonals from 3D variables
-            This is calling Gerd's CollectData1D() routine which can
-            extract non-axis-parallel lines too but is fixed to 3D data. */
-  if (errormsg && ! strcmp (errormsg, "Given normal vector isn't axis-parallel"))
-  {
-    int length;
-
-
-    /* get the info on the variable to extract a hyperslab from */
-    CCTK_GroupData (CCTK_GroupIndexFromVarI (vindex), &vinfo);
-
-    if (hdim != 1 || vinfo.dim != 3)
-    {
-      CCTK_WARN (1, "Non-axis-parallel hyperslabs are supported as 1D lines "
-                    "from 3D variables only");
-      return (-1);
-    }
-
-    length = extents[0];
-    if (length > 0)
-    {
-      length--;
-    }
-    return (Hyperslab_CollectData1D (GH, vindex, vtimelvl, global_origin,
-                                     directions, downsample_[0], length,
-                                     hdata, hsize, target_proc));
-  }
-
-  /* immediately return in case of errors */
-  if (errormsg)
-  {
-    CCTK_WARN (1, errormsg);
-    return (-1);
-  }
-
   /* get processor information */
   myproc = CCTK_MyProc (GH);
   nprocs = CCTK_nProcs (GH);
@@ -812,6 +762,55 @@ int Hyperslab_GetHyperslab (const cGH *GH,
   /* get the info on the variable to extract a hyperslab from */
   CCTK_GroupData (CCTK_GroupIndexFromVarI (vindex), &vinfo);
 
+  /* do some plausibility checks */
+  errormsg = checkParameters (GH, vindex, vtimelvl, hdim,
+                              global_origin, directions, extents,
+                              downsample_, hdata, hsize);
+
+  /* FIXME: hack for getting diagonals from 3D variables */
+  if (errormsg && ! strcmp (errormsg, "Given normal vector isn't axis-parallel"))
+  {
+    int mapping;
+
+
+    if (hdim != 1 || vinfo.dim != 3)
+    {
+      CCTK_WARN (1, "Non-axis-parallel hyperslabs are supported as diagonals "
+                    "from non-staggered 3D variables only");
+      return (-1);
+    }
+
+    mapping = Hyperslab_DefineGlobalMappingByIndex (GH, vindex, 1, directions,
+                                                    global_origin, extents,
+                                                    downsample_, -1,
+                                                    target_proc, NULL, hsize);
+    if (mapping < 0)
+    {
+      CCTK_WARN (1, "Failed to define hyperslab mapping for 3D diagonal");
+      return (-1);
+    }
+    if (hsize[0] > 0)
+    {
+      *hdata = malloc (hsize[0] * CCTK_VarTypeSize (vinfo.vartype));
+      retval = Hyperslab_Get (GH, mapping, vindex, 0, vinfo.vartype, *hdata);
+    }
+    else
+    {
+      *hdata = NULL;
+      retval = 0;
+    }
+    Hyperslab_FreeMapping (mapping);
+
+    return (retval);
+  }
+
+  /* immediately return in case of errors */
+  if (errormsg)
+  {
+    CCTK_WARN (1, errormsg);
+    return (-1);
+  }
+
   /* set the pointers to pass to Hyperslab_GetLocalHyperslab() */
   if (nprocs == 1)
   {
diff --git a/src/Mapping.c b/src/Mapping.c
index c74be93..8e80893 100644
--- a/src/Mapping.c
+++ b/src/Mapping.c
@@ -218,6 +218,17 @@ CCTK_INT Hyperslab_DefineGlobalMappingByIndex (
       return (-7);
     }
   }
+
+  /* diagonals can be extracted from non-staggered 3D variables only */ 
+  if (mapping->is_diagonal_in_3D && vinfo.stagtype != 0)
+  {
+    free (mapping->vectors); free (mapping);
+    CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+                "Hyperslab_DefineGlobalMappingByIndex: diagonals can be "
+                "extracted from non-staggered 3D variables only");
+    return (-7);
+  }
+
   for (vdim = 0; vdim < (unsigned int) vinfo.dim; vdim++)
   {
     mapping->do_dir[vdim] = 0;
@@ -256,6 +267,15 @@ CCTK_INT Hyperslab_DefineGlobalMappingByIndex (
       }
       hdim++;
     }
+    else if (mapping->is_diagonal_in_3D &&
+             origin[vdim] + extent[0] > GA->extras->nsize[vdim])
+    {
+      free (mapping->vectors); free (mapping);
+      CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+                  "Hyperslab_DefineGlobalMappingByIndex: extent in "
+                  "%d-direction exceeds grid size", vdim);
+      return (-8);
+    }
   }
 
   /* now fill out the hyperslab mapping structure */
@@ -281,6 +301,23 @@ CCTK_INT Hyperslab_DefineGlobalMappingByIndex (
       }
       hdim++;
     }
+    else if (mapping->is_diagonal_in_3D)
+    {
+      mapping->totals = extent[0] / mapping->downsample[0];
+      if (extent[0] % mapping->downsample[0])
+      {
+        mapping->totals++;
+      }
+      /* subtract ghostzones for periodic BC */
+      if (GA->connectivity->perme[vdim])
+      {
+        mapping->totals -= 2 * GA->extras->nghostzones[vdim];
+      }
+      if ((unsigned int) mapping->global_hsize[0] > mapping->totals)
+      {
+        mapping->global_hsize[0] = mapping->totals;
+      }
+    }
   }
 
   /* check whether the full local data patch was requested as hyperslab */
@@ -293,17 +330,12 @@ CCTK_INT Hyperslab_DefineGlobalMappingByIndex (
   }
   else if (mapping->is_diagonal_in_3D)
   {
-    /* diagonals in 3D are special: global_hsize is minimum of grid extents */
-    mapping->global_hsize[0] = GH->cctk_gsh[0];
-    if (mapping->global_hsize[0] > GH->cctk_gsh[1])
-    {
-      mapping->global_hsize[0] = GH->cctk_gsh[1];
-    }
-    if (mapping->global_hsize[0] > GH->cctk_gsh[2])
+    /* just initialize the downsample and global_startpoint vectors */
+    for (vdim = 0; vdim < (unsigned int) vinfo.dim; vdim++)
     {
-      mapping->global_hsize[0] = GH->cctk_gsh[2];
+      mapping->downsample[vdim] = mapping->downsample[0];
+      mapping->global_startpoint[vdim] = origin[vdim];
     }
-    mapping->totals = mapping->global_hsize[0];
   }
   else
   {
@@ -371,7 +403,6 @@ CCTK_INT Hyperslab_DefineGlobalMappingByIndex (
         mapping->local_endpoint[vdim] = -1;
       }
 
-#define DEBUG 1
 #ifdef DEBUG
       printf ("direction %d: global ranges [%d, %d), local ranges[%d, %d)\n",
               vdim,
diff --git a/src/make.code.defn b/src/make.code.defn
index 7d322e8..d3c05c7 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -2,4 +2,4 @@
 # $Header$
 
 # Source files in this directory
-SRCS = Mapping.c GetHyperslab.c Hyperslab.c NewHyperslab.c DatatypeConversion.c CollectData1D.c GetLocalLine.c
+SRCS = Mapping.c GetHyperslab.c Hyperslab.c NewHyperslab.c DatatypeConversion.c