This commit was generated by cvs2svn to compensate for changes in r2, which

included commits to RCS files with non-trunk default branches.

git-svn-id: http://svn.cactuscode.org/arrangements/CactusPUGH/PUGHInterp/trunk@3 1c20744c-e24a-42ec-9533-f5004cb800e5

5 files changed, 2049 insertions, 0 deletions

diff --git a/src/Interpolate.c b/src/Interpolate.c
new file mode 100644
index 0000000..392cb18
--- /dev/null
+++ b/src/Interpolate.c
@@ -0,0 +1,427 @@
+ /*@@
+   @file      Interpolate.c
+   @date      Wed 17 Jan 2001
+   @author    Thomas Radke
+   @desc
+              Interpolation of arrays to arbitrary points
+
+              This interpolator is based on the Cactus 3.x Fortran version
+              written by Paul Walker. It also contains some nice optimization
+              features from Erik Schnetter.
+   @enddesc
+   @history
+   @date      Wed 17 Jan 2001
+   @author    Thomas Radke
+   @hdesc     Translation from Fortran to C
+   @endhistory
+   @version   $Id$
+ @@*/
+
+#include <stdlib.h>
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+#include "pughInterpGH.h"
+
+/* the rcs ID and its dummy function to use it */
+static char *rcsid = "$Id$";
+CCTK_FILEVERSION(CactusPUGH_PUGHInterp_Interpolate_c)
+
+
+/* the highest order of interpolation we support so far */
+#define MAXORDER  3
+
+/* the highest dimension for variables we can deal with (so far) */
+#define MAXDIM    3
+
+/* macro to do the interpolation of in_array[n] to out_array[n] at point[] */
+#define INTERPOLATE(cctk_type, cctk_subtype, subpart, in_array, out_array,    \
+                    order, point, dims, n, coeff)                             \
+    {                                                                         \
+      int ii, jj, kk;                                                         \
+      cctk_type *fi;                                                          \
+      /* for some reason (probably loop unrolling & pipelining) the compiler  \
+         produces faster code if arrays are used instead of scalars */        \
+      cctk_subtype interp_result, fj[1], fk[1];                               \
+                                                                              \
+                                                                              \
+      interp_result = 0;                                                      \
+                                                                              \
+      /* NOTE: support >3D arrays by adding more loops */                     \
+      for (kk = 0; kk <= order; kk++)                                         \
+      {                                                                       \
+        fk[0] = 0;                                                            \
+        for (jj = 0; jj <= order; jj++)                                       \
+        {                                                                     \
+          /* NOTE: for >3D arrays adapt the index calculation here */         \
+          fi = (cctk_type *) in_array +                                       \
+               point[0] + dims[0]*(point[1]+jj + dims[1]*(point[2]+kk));      \
+          fj[0] = 0;                                                          \
+          for (ii = 0; ii <= order; ii++)                                     \
+          {                                                                   \
+            fj[0] += fi[ii]subpart * coeff[0][ii];                            \
+          }                                                                   \
+          fk[0] += fj[0] * coeff[1][jj];                                      \
+        }                                                                     \
+        interp_result += fk[0] * coeff[2][kk];                                \
+      }                                                                       \
+                                                                              \
+      /* assign the result */                                                 \
+      ((cctk_type *) out_array)[n]subpart = interp_result;                    \
+    }
+
+/* this is needed by some preprocessors to pass into INTERPOLATE
+   as a dummy macro */
+#define NOTHING
+
+
+/*@@
+   @routine    PUGHInterp_Interpolate
+   @date       Wed 17 Jan 2001
+   @author     Thomas Radke
+   @desc
+               This routine interpolates a set of input arrays
+               to a set of output arrays (one-to-one) at arbitrary points\                     which are given by their coordinates and the underlying
+               regular, uniform grid.
+
+               Current limitations of this implementation are:
+
+                - arrays up to three (MAXDIM) dimensions only can be handled
+                - interpolation orders up to three (MAXORDER) only are supported
+                - coordinates must be given as CCTK_REAL types
+                - input and output array types must be the same
+                  (no type casting of interpolation results supported)
+
+               Despite of these limitations, the code was programmed almost
+               generic in that it can easily be extended to support higher-
+               dimensional arrays or more interpolation orders.
+               Please see the NOTES in this source file for details.
+   @enddesc
+
+   @var        num_points
+   @vdesc      number of points to interpolate at
+   @vtype      int
+   @vio        in
+   @endvar
+   @var        num_dims
+   @vdesc      dimensionality of the input arrays
+   @vtype      int
+   @vio        in
+   @endvar
+   @var        num_arrays
+   @vdesc      number of input/output arrays
+   @vtype      int
+   @vio        in
+   @endvar
+   @var        dims
+   @vdesc      dimensions of the input arrays
+   @vtype      int[ num_dims ]
+   @vio        in
+   @endvar
+   @var        coord
+   @vdesc      coordinates to interpolate at
+   @vtype      CCTK_REAL coord[ num_dims * num_points ]
+   @vio        in
+   @endvar
+   @var        origin
+   @vdesc      origin of the underlying grid
+   @vtype      CCTK_REAL origin[ num_dims ]
+   @vio        in
+   @endvar
+   @var        delta
+   @vdesc      deltas of the underlying grid
+   @vtype      CCTK_REAL delta[ num_dims ]
+   @vio        in
+   @endvar
+   @var        in_types
+   @vdesc      CCTK variable types of input arrays
+   @vtype      int in_types[ num_arrays ]
+   @vio        in
+   @endvar
+   @var        in_arrays
+   @vdesc      list of input arrays
+   @vtype      void *in_arrays[ num_arrays ]
+   @vio        in
+   @endvar
+   @var        out_types
+   @vdesc      CCTK variable types of output arrays
+   @vtype      int out_types[ num_arrays ]
+   @vio        in
+   @endvar
+   @var        out_arrays
+   @vdesc      list of output arrays
+   @vtype      void *out_arrays[ num_arrays ]
+   @vio        out
+   @endvar
+
+   @returntype int
+   @returndesc
+               0  - successful interpolation
+               -1 - in case of any errors
+   @endreturndesc
+@@*/
+int PUGHInterp_Interpolate (int order,
+                            int num_points,
+                            int num_dims,
+                            int num_arrays,
+                            int dims[ /* num_dims */ ],
+                            CCTK_REAL coord[ /* num_dims * num_points */ ],
+                            CCTK_REAL origin[ /* num_dims */ ],
+                            CCTK_REAL delta[ /* num_dims */ ],
+                            int in_types[ /* num_arrays */ ],
+                            void *in_arrays[ /* num_arrays */ ],
+                            int out_types[ /* num_arrays */ ],
+                            void *out_arrays[ /* num_arrays */ ])
+{
+  int retval;
+  int i, a, n, out_of_bounds, shift;
+  int max_dims[MAXDIM], point[MAXDIM];
+  CCTK_REAL below[MAXDIM];
+  CCTK_REAL offset[MAXDIM];
+  CCTK_REAL coeff[MAXDIM][MAXORDER + 1];
+
+
+  /* verify parameters and check against our restrictions */
+  retval = -1;
+  if (num_dims < 1)
+  {
+    CCTK_WARN (1, "Number of dimensions must be positive");
+  }
+  else if (num_dims > MAXDIM)
+  {
+    CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+                "Interpolation of %d-dimensional arrays not implemented",
+                num_dims);
+  }
+  else if (order < 1)
+  {
+    CCTK_WARN (1, "Inperpolation order must be positive");
+  }
+  else if (order > MAXORDER)
+  {
+    CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+                "Interpolation order %d not implemented", order);
+  }
+  else if (num_points < 0)
+  {
+    CCTK_WARN (1, "Negative number of points given");
+  }
+  else
+  {
+    retval = 0;
+  }
+
+  /* immediately return in case of errors */
+  if (retval)
+  {
+    return (retval);
+  }
+
+  /* also immediately return if there's nothing to do */
+  if (num_points == 0)
+  {
+    return (retval);
+  }
+
+  /* duplicate the dims[] vector into one with MAXDIM-1 elements
+     (with the remaining elements zeroed out)
+     so that we can use nested loops over MAXDIM dimensions later on */
+  memset (max_dims, 0, sizeof (max_dims));
+  memcpy (max_dims, dims, (num_dims - 1) * sizeof (int));
+
+  /* zero out the coefficients and set the elements with index 'order' to one
+     so that we can use nested loops over MAXDIM dimensions later on */
+  memset (coeff, 0, sizeof (coeff));
+  for (i = num_dims; i < MAXDIM; i++)
+  {
+    coeff[i][0] = 1;
+  }
+
+  /* zero out the iterator */
+  memset (point, 0, sizeof (point));
+
+  /* loop over all points to interpolate at */
+  for (n = 0; n < num_points; n++)
+  {
+    /* reset the out-of-bounds flag */
+    out_of_bounds = 0;
+
+    /* loop over all dimensions */
+    for (i = 0; i < num_dims; i++)
+    {
+      /* nearest grid point for stencil */
+      point[i] = (coord[num_dims*n + i] - origin[i]) / delta[i];
+
+      /* if beyond upper bound shift the grid point to the left */
+      shift = point[i] + order - (dims[i] - 1);
+      if (shift > 0)
+      {
+        point[i] -= shift;
+      }
+
+      /* physical coordinate of that grid point */
+      below[i] = origin[i] + point[i] * delta[i];
+
+      /* offset from that grid point, in fractions of grid points */
+      offset[i] = (coord[num_dims*n + i] - below[i]) / delta[i];
+
+      /* test bounds */
+      out_of_bounds |= point[i] < 0 || point[i]+order >= dims[i];
+    }
+
+    /* check bounds */
+    if (out_of_bounds)
+    {
+      char *msg;
+
+
+      /* put all information into a single message string for output */
+      msg = (char *) malloc (100 + num_dims*(10 + 4*20));
+      sprintf (msg, "Interpolation stencil out of bounds at grid point [%d",
+               point[0]);
+      for (i = 1; i < num_dims; i++)
+      {
+        sprintf (msg, "%s, %d", msg, point[i]);
+      }
+      sprintf (msg, "%s]\nrange would be min/max [%f / %f", msg,
+               (double) below[0], (double) (below[0] + offset[0]));
+      for (i = 1; i < num_dims; i++)
+      {
+        sprintf (msg, "%s, %f / %f", msg,
+                 (double) below[i], (double) (below[i] + offset[i]));
+      }
+      sprintf (msg, "%s]\ngrid is min/max        [%f / %f", msg,
+               (double) origin[0], (double) (origin[0] + dims[0]*delta[0]));
+      for (i = 1; i < num_dims; i++)
+      {
+        sprintf (msg, "%s, %f / %f", msg,
+                 (double) origin[i], (double) (origin[i] + dims[i]*delta[i]));
+      }
+      sprintf (msg, "%s]", msg);
+      CCTK_WARN (1, msg);
+      free (msg);
+
+      continue;
+    }
+
+    /* compute the interpolation coefficients according to the order
+       Thanks to Erik for formulating these so nicely. */
+    /* NOTE: support higher interpolation orders by adding the
+             appropriate coefficients in another else branch */
+    if (order == 1)
+    {
+      /* first order (linear) 1D interpolation */
+      for (i = 0; i < num_dims; i++)
+      {
+        coeff[i][0] = 1 - offset[i];
+        coeff[i][1] =     offset[i];
+      }
+    }
+    else if (order == 2)
+    {
+      /* second order (quadratic) 1D interpolation */
+      for (i = 0; i < num_dims; i++)
+      {
+        coeff[i][0] = (1-offset[i]) * (2-offset[i]) / (  2  *   1 );
+        coeff[i][1] = ( -offset[i]) * (2-offset[i]) / (  1  * (-1));
+        coeff[i][2] = ( -offset[i]) * (1-offset[i]) / ((-1) * (-2));
+      }
+    }
+    else if (order == 3)
+    {
+      /* third order (cubic) 1D interpolation */
+      for (i = 0; i < num_dims; i++)
+      {
+        coeff[i][0] = (1-offset[i]) * (2-offset[i]) * (3-offset[i]) /
+                      (  3  *   2  *   1 );
+        coeff[i][1] = ( -offset[i]) * (2-offset[i]) * (3-offset[i]) /
+                      (  2  *   1  * (-1));
+        coeff[i][2] = ( -offset[i]) * (1-offset[i]) * (3-offset[i]) /
+                      (  1  * (-1) * (-2));
+        coeff[i][3] = ( -offset[i]) * (1-offset[i]) * (2-offset[i]) /
+                      ((-1) * (-2) * (-3));
+      }
+    }
+    else
+    {
+      CCTK_WARN (0, "Implementation error");
+    }
+
+    /* now loop over all arrays to interpolate at the current point */
+    for (a = 0; a < num_arrays; a++)
+    {
+      /* sorry, for now input and output arrays must be of same type */
+      if (in_types[a] != out_types[a])
+      {
+        CCTK_WARN (1, "Type casting of interpolation results not implemented");
+        continue;
+      }
+
+      /* now do the interpolation according to the array type
+         we support all kinds of CCTK_REAL* and CCTK_COMPLEX* types here */
+      if (in_types[a] == CCTK_VARIABLE_REAL)
+      {
+        INTERPOLATE (CCTK_REAL, CCTK_REAL, NOTHING, in_arrays[a],
+                     out_arrays[a], order, point, max_dims, n, coeff);
+      }
+      else if (in_types[a] == CCTK_VARIABLE_COMPLEX)
+      {
+        INTERPOLATE (CCTK_COMPLEX, CCTK_REAL, .Re, in_arrays[a],
+                     out_arrays[a], order, point, max_dims, n, coeff);
+        INTERPOLATE (CCTK_COMPLEX, CCTK_REAL, .Im, in_arrays[a],
+                     out_arrays[a], order, point, max_dims, n, coeff);
+      }
+#ifdef CCTK_REAL4
+      else if (in_types[a] == CCTK_VARIABLE_REAL4)
+      {
+        INTERPOLATE (CCTK_REAL4, CCTK_REAL4, NOTHING, in_arrays[a],
+                     out_arrays[a], order, point, max_dims, n, coeff);
+      }
+      else if (in_types[a] == CCTK_VARIABLE_COMPLEX8)
+      {
+        INTERPOLATE (CCTK_COMPLEX8, CCTK_REAL4, .Re, in_arrays[a],
+                     out_arrays[a], order, point, max_dims, n, coeff);
+        INTERPOLATE (CCTK_COMPLEX8, CCTK_REAL4, .Im, in_arrays[a],
+                     out_arrays[a], order, point, max_dims, n, coeff);
+      }
+#endif
+#ifdef CCTK_REAL8
+      else if (in_types[a] == CCTK_VARIABLE_REAL8)
+      {
+        INTERPOLATE (CCTK_REAL8, CCTK_REAL8, NOTHING, in_arrays[a],
+                     out_arrays[a], order, point, max_dims, n, coeff);
+      }
+      else if (in_types[a] == CCTK_VARIABLE_COMPLEX16)
+      {
+        INTERPOLATE (CCTK_COMPLEX16, CCTK_REAL8, .Re, in_arrays[a],
+                     out_arrays[a], order, point, max_dims, n, coeff);
+        INTERPOLATE (CCTK_COMPLEX16, CCTK_REAL8, .Im, in_arrays[a],
+                     out_arrays[a], order, point, max_dims, n, coeff);
+      }
+#endif
+#ifdef CCTK_REAL16
+      else if (in_types[a] == CCTK_VARIABLE_REAL16)
+      {
+        INTERPOLATE (CCTK_REAL16, CCTK_REAL16, NOTHING, in_arrays[a],
+                     out_arrays[a], order, point, max_dims, n, coeff);
+      }
+      else if (in_types[a] == CCTK_VARIABLE_COMPLEX32)
+      {
+        INTERPOLATE (CCTK_COMPLEX32, CCTK_REAL16, .Re, in_arrays[a],
+                     out_arrays[a], order, point, max_dims, n, coeff);
+        INTERPOLATE (CCTK_COMPLEX32, CCTK_REAL16, .Im, in_arrays[a],
+                     out_arrays[a], order, point, max_dims, n, coeff);
+      }
+#endif
+      else
+      {
+        CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+                    "Unsupported variable type %d", (int) in_types[a]);
+      }
+    } /* end of loop over all arrays */
+
+  } /* end of loop over all points to interpolate at */
+
+  /* we're done */
+  return (retval);
+}
diff --git a/src/Operator.c b/src/Operator.c
new file mode 100644
index 0000000..76a71f5
--- /dev/null
+++ b/src/Operator.c
@@ -0,0 +1,1293 @@
+ /*@@
+   @file      Operator.c
+   @date      Tue Apr 15 18:22:45 1997
+   @author    Paul Walker
+   @desc 
+              Definition of interpolation operators for regular uniform grids.
+   @enddesc 
+
+   @history
+   @date      Sun Jul 04 1999
+   @author    Thomas Radke
+   @hdesc     conversion to Cactus 4.0 (copied from pughGetPoints.c)
+   @date      Wed 31 Jan 2001
+   @author    Thomas Radke
+   @hdesc     translation of fortran interpolators into C
+   @endhistory
+   @version   $Id$
+ @@*/
+
+#include <stdlib.h>
+#include <math.h>
+#include <string.h>
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+#include "CactusPUGH/PUGH/src/include/pugh.h"
+#include "pughInterpGH.h"
+
+/* the rcs ID and its dummy function to use it */
+static char *rcsid = "$Header$";
+CCTK_FILEVERSION(CactusPUGH_PUGHInterp_Operator_c)
+
+/* uncomment this to get some debugging output */
+/* #define PUGHINTERP_DEBUG 1 */
+
+/* macro do sort interpolation results from a single communication buffer
+   into their appropriate output arrays */
+#define SORT_TYPED_ARRAY(cctk_type, source, dest, points, index_array, offset)\
+          {                                                                   \
+            int i;                                                            \
+            cctk_type *src, *dst;                                             \
+                                                                              \
+                                                                              \
+            src = (cctk_type *) source;                                       \
+            dst = (cctk_type *) dest;                                         \
+            for (i = 0; i < points; i++)                                      \
+            {                                                                 \
+              dst[index_array[i + offset]] = *src++;                          \
+            }                                                                 \
+            source = (char *) src;                                            \
+          }
+
+
+
+#ifdef CCTK_MPI
+/* internal structure describing a handle for a single CCTK data type */
+typedef struct
+{
+  int vtypesize;              /* variable type's size in bytes */
+  MPI_Datatype mpi_type;      /* corresponding MPI datatype */
+  int num_arrays;             /* number of in/out arrays */
+  int array;                  /* iterator for num_arrays */
+  void *send_buffer;          /* communication send buffer for this type */
+  void *recv_buffer;          /* communication receive buffer for this type */
+  char *buffer;               /* work pointer for send_buffer */
+} type_desc_t;
+#endif
+
+
+/* prototypes of routines defined in this source file */
+static int PUGHInterp_CheckArguments (cGH *GH,
+                                      int num_dims,
+                                      int num_points,
+                                      int num_in_arrays,
+                                      int num_out_arrays,
+                                      int interp_coord_array_types[]);
+#ifdef CCTK_MPI
+static int GetLocalCoords (cGH *GH,
+                           int num_points,
+                           const char *coord_system,
+                           pGExtras *extras,
+                           CCTK_REAL *coords[],
+                           int *num_local_points,
+                           CCTK_REAL **local_coords);
+#endif
+
+
+/*@@
+  @routine    PUGHInterp_InterpGV
+  @date       Sun Jul 04 1999
+  @author     Thomas Radke
+  @desc
+              The interpolation operator registered with the CCTK
+              under the name "regular uniform cartesian".
+
+              Interpolates a list of CCTK variables (domain-decomposed
+              grid functions or arrays) to a list of output arrays
+              (one-to-one) at a given number of interpolation points
+              (indicated by their coordinates). The points are located
+              on a coordinate system which is assumed to be a uniform
+              cartesian.
+  @enddesc
+
+  @var        GH
+  @vdesc      Pointer to CCTK grid hierarchy
+  @vtype      cGH *
+  @vio        in
+  @endvar
+  @var        order
+  @vdesc      interpolation order
+  @vtype      int
+  @vio        in
+  @endvar
+  @var        coord_system
+  @vdesc      name of coordinate system to use for interpolation
+  @vtype      const char *
+  @vio        in
+  @endvar
+  @var        num_points
+  @vdesc      number of points to be interpolated on this processor
+  @vtype      int
+  @vio        in
+  @endvar
+  @var        num_in_array_indices
+  @vdesc      number of input arrays (given by their indices)
+              to interpolate from
+  @vtype      int
+  @vio        in
+  @endvar
+  @var        num_out_arrays
+  @vdesc      number of output arrays to interpolate to
+  @vtype      int
+  @vio        in
+  @endvar
+  @var        interp_coord_arrays
+  @vdesc      coordinates of points to interpolate at
+  @vtype      void *[num_dims]
+  @vio        in
+  @endvar
+  @var        interp_coord_array_types
+  @vdesc      CCTK data type of coordinate arrays
+  @vtype      int [num_dims]
+  @vio        in
+  @endvar
+  @var        in_arrays
+  @vdesc      list of input arrays to interpolate on
+  @vtype      void *[num_in_arrays]
+  @vio        in
+  @endvar
+  @var        in_array_types
+  @vdesc      CCTK data types of input arrays
+  @vtype      int [num_in_arrays]
+  @vio        in
+  @endvar
+  @var        out_arrays
+  @vdesc      list of output arrays to interpolate to
+  @vtype      void *[num_out_arrays]
+  @vio        out
+  @endvar
+  @var        out_array_types
+  @vdesc      CCTK data types of output arrays
+  @vtype      int [num_out_arrays]
+  @vio        in
+  @endvar
+
+  @returntype int
+  @returndesc
+              0  - successful interpolation
+              -1 - in case of any errors
+  @endreturndesc
+@@*/
+int PUGHInterp_InterpGV (cGH *GH,
+                         int order,
+                         const char *coord_system,
+                         int num_points,
+                         int num_in_array_indices,
+                         int num_out_arrays,
+                         void *interp_coord_arrays[],
+                         int interp_coord_array_types[],
+                         int in_array_indices[],
+                         void *out_arrays[],
+                         int out_array_types[])
+{
+  CCTK_REAL **data, *interp_local_coords;
+  CCTK_REAL *origin, *delta;
+  int nprocs;
+  int timelevel;
+  int dim, num_dims, *dims;
+  int array;
+  int point;
+  int retval;
+  void **in_arrays;
+  pGH *pughGH;
+  pGExtras *extras;
+  cGroupDynamicData group_data;
+#ifdef CCTK_MPI
+  int num_local_points;
+  int offset;
+  int proc;
+  int type;
+  void **local_out_arrays;
+  pughInterpGH *myGH;
+  int max_type;
+  type_desc_t *type_desc;
+#endif
+
+
+  /* get dimensionality of the coordinate system */
+  num_dims = CCTK_CoordSystemDim (coord_system);
+  if (num_dims <= 0)
+  {
+    CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+                "Cannot get dimensions of coordinate system '%s'",
+                coord_system);
+    return (-1);
+  }
+
+  /* check other arguments */
+  retval = PUGHInterp_CheckArguments (GH, num_dims, num_points,
+                                      num_in_array_indices, num_out_arrays,
+                                      interp_coord_array_types);
+  if (retval <= 0)
+  {
+    return (retval);
+  }
+
+  /* get dimensions, origin, and delta of the processor-local grid */
+  /* NOTE: getting the dimensions should be a flesh routine as well
+           for now we get the dimensions of every coordinate and take the
+           i'th element - this is inconsistent !! */
+  dims = (int *) malloc (2 * num_dims * sizeof (int));
+  origin = (CCTK_REAL *) malloc (2 * num_dims * sizeof (CCTK_REAL));
+  delta  = origin + num_dims;
+  in_arrays = malloc (num_in_array_indices * sizeof (void *));
+  for (dim = 0; dim < num_dims; dim++)
+  {
+    CCTK_GrouplshVI (GH, num_dims, dims + num_dims,
+                     CCTK_CoordIndex (dim + 1, NULL, coord_system));
+    dims[dim] = dims[dim + num_dims];
+
+    CCTK_CoordLocalRange (GH, &origin[dim], &delta[dim], dim + 1, NULL,
+                          coord_system);
+    delta[dim] = (delta[dim] - origin[dim]) / dims[dim];
+  }
+
+  /* get extension handle for PUGH */
+  pughGH = (pGH *) CCTK_GHExtension (GH, "PUGH");
+  extras = NULL;
+
+  /* check that the input arrays dimensions match the coordinate system
+     (but their dimensionality can be less) */
+  for (array = 0; array < num_in_array_indices; array++)
+  {
+    if (CCTK_GroupDynamicData (GH,
+                               CCTK_GroupIndexFromVarI(in_array_indices[array]),
+                               &group_data) < 0)
+    {
+      CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+                  "Invalid input array index %d",
+                  in_array_indices[array]);
+      retval = -1;
+      continue;
+    }
+
+    if (group_data.dim > num_dims)
+    {
+      CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+                  "Input array variable with index %d has more dimensions "
+                  "than coordinate system '%s'",
+                  in_array_indices[array], coord_system);
+      retval = -1;
+      continue;
+    }
+
+    if (memcmp (group_data.lsh, dims, group_data.dim * sizeof (int)))
+    {
+      CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+                  "Dimensions of input array variable with index %d "
+                  "doesn't match with coordinate system '%s'",
+                  in_array_indices[array], coord_system);
+      retval = -1;
+    }
+
+    /* the extras pointer will refer to the first input array with the
+       highest dimension (need this information later on) */
+    if (! extras || extras->dim < group_data.dim)
+    {
+      extras = ((pGA *) pughGH->variables[in_array_indices[array]][0])->extras;
+    }
+
+    /* get the data pointer to the input array (use current timelevel) */
+    timelevel = CCTK_NumTimeLevelsFromVarI (in_array_indices[array]) - 1;
+    if (timelevel > 0)
+    {
+      timelevel--;
+    }
+    in_arrays[array] = CCTK_VarDataPtrI(GH, timelevel, in_array_indices[array]);
+  }
+  if (retval < 0)
+  {
+    free (in_arrays);
+    free (origin);
+    free (dims);
+    return (retval);
+  }
+
+  /* single-processor case is easy: no communication or buffering, just direct
+     interpolation of interp_coord_arrays from in_arrays into out_arrays */
+  nprocs = CCTK_nProcs (GH);
+  if (nprocs == 1)
+  {
+    /* sort the individual interpolation coordinate arrays into a single one */
+    interp_local_coords = (CCTK_REAL *) malloc (num_dims * num_points *
+                                                sizeof (CCTK_REAL));
+    data = (CCTK_REAL **) interp_coord_arrays;
+    for (point = 0; point < num_points; point++)
+    {
+      for (dim = 0; dim < num_dims; dim++)
+      {
+        *interp_local_coords++ = data[dim][point];
+      }
+    }
+    interp_local_coords -= num_dims * num_points;
+
+    /* call the interpolator function */
+    retval = PUGHInterp_Interpolate (order,
+                                     num_points, num_dims, num_out_arrays,
+                                     dims, interp_local_coords,
+                                     origin, delta,
+                                     out_array_types, in_arrays,
+                                     out_array_types, out_arrays);
+
+    /* free allocated resources */
+    free (interp_local_coords);
+    free (in_arrays);
+    free (origin);
+    free (dims);
+
+    return (retval);
+  }
+
+#ifdef CCTK_MPI
+  /*** Here follows the multi-processor case:
+       All processors locate their points to interpolate at
+       and exchange the coordinates so that every processor gets
+       those points which it can process locally.
+       After interpolation the results have to be send back to the
+       requesting processors.
+       For both communications MPI_Alltoallv() is used.
+
+       In order to minimize the total number of MPI_Alltoallv() calls
+       (which are quite expensive) we collect the interpolation results
+       for all output arrays of the same CCTK data type into a single
+       communication buffer. That means, after communication the data
+       needs to be resorted from the buffer into the output arrays.
+   ***/
+
+  /* first of all, set up a structure with information of the
+     CCTK data types we have to deal with */
+
+  /* get the maximum value of the output array CCTK data types
+     NOTE: we assume that CCTK data types are defined as consecutive
+           positive constants starting from zero */
+  max_type = 0;
+  for (array = 0; array < num_out_arrays; array++)
+  {
+    if (max_type < out_array_types[array])
+    {
+      max_type = out_array_types[array];
+    }
+  }
+
+  /* now allocate an array of structures for all potential types */
+  type_desc = (type_desc_t *) calloc (max_type + 1, sizeof (type_desc_t));
+
+  /* count the number of arrays of same type
+     (the num_arrays element was already initialized to zero by calloc() */
+  for (array = 0; array < num_out_arrays; array++)
+  {
+    type_desc[out_array_types[array]].num_arrays++;
+  }
+
+  /* fill in the type description information */
+  retval = 0;
+  for (type = 0; type <= max_type; type++)
+  {
+    if (type_desc[type].num_arrays > 0)
+    {
+      /* get the variable type size in bytes */
+      type_desc[type].vtypesize = CCTK_VarTypeSize (type);
+      if (type_desc[type].vtypesize <= 0)
+      {
+        CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+                    "Invalid variable type %d passed, "
+                    "arrays of such type will be skipped during interpolation",
+                    type);
+        type_desc[type].num_arrays = 0;
+        continue;
+      }
+
+      /* get the MPI data type to use for communicating such a CCTK data type */
+      type_desc[type].mpi_type = PUGH_MPIDataType (pughGH, type);
+      if (type_desc[type].mpi_type == MPI_DATATYPE_NULL)
+      {
+        CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+                    "No MPI data type defined for variable type %d, "
+                    "arrays of such type will be skipped during interpolation",
+                    type);
+        type_desc[type].num_arrays = 0;
+        continue;
+      }
+
+      retval++;
+    }
+  }
+
+  /* check that there's at least one array with a valid CCTK data type */
+  if (retval <= 0)
+  {
+    free (in_arrays);
+    free (dims);
+    free (origin);
+    free (type_desc);
+    return (-1);
+  }
+
+  /* map the requested points to interpolate at onto the processors
+     they belong to and gather the coordinates of all the points
+     that this processor owns
+     the number of processor-local points is returned in num_local_points,
+     their coordinates in interp_local_coords */
+  retval = GetLocalCoords (GH, num_points, coord_system, extras,
+                           (CCTK_REAL **) interp_coord_arrays,
+                           &num_local_points, &interp_local_coords);
+  if (retval)
+  {
+    free (in_arrays);
+    free (origin);
+    free (dims);
+    free (type_desc);
+    return (retval);
+  }
+
+  /* allocate contiguous communication buffers for each CCTK data type
+     holding the local interpolation results from all input arrays
+     of that type
+     If there are no points to process on this processor
+     set the buffer pointer to an invalid but non-NULL value
+     otherwise we might get trouble with NULL pointers in MPI_Alltoallv () */
+  for (type = 0; type <= max_type; type++)
+  {
+    if (type_desc[type].num_arrays > 0 && num_local_points > 0)
+    {
+      type_desc[type].send_buffer = malloc (num_local_points *
+                                            type_desc[type].num_arrays *
+                                            type_desc[type].vtypesize);
+      type_desc[type].buffer = (char *) type_desc[type].send_buffer;
+      type_desc[type].array = 0;
+    }
+    else
+    {
+      /* dereferencing such an address should code crash on most systems */
+      type_desc[type].send_buffer = (void *) 1;
+    }
+  }
+
+  /* get extension handle for interp */
+  myGH = (pughInterpGH *) CCTK_GHExtension (GH, "PUGHInterp");
+
+  /* allocate new out_arrays array for local interpolation results
+     from this processor */
+  local_out_arrays = calloc (num_out_arrays, sizeof (void *));
+
+  /* now, in a loop over all processors, do the interpolation
+     and put the results in the communication buffer at the proper offset */
+  for (proc = 0; proc < nprocs; proc++)
+  {
+    for (type = 0; type <= max_type; type++)
+    {
+      if (type_desc[type].num_arrays > 0)
+      {
+        for (array = 0; array < num_out_arrays; array++)
+        {
+          if (out_array_types[array] == type)
+          {
+            local_out_arrays[array] = type_desc[type].buffer;
+            type_desc[type].buffer += myGH->num_points_to[proc] *
+                                      type_desc[type].vtypesize;
+          }
+        }
+      }
+    }
+
+    /* call the interpolation operator to process all points of all
+       output arrays for this processor */
+    PUGHInterp_Interpolate (order,
+                            myGH->num_points_to[proc], num_dims, num_out_arrays,
+                            dims, interp_local_coords, origin, delta,
+                            out_array_types, in_arrays,
+                            out_array_types, local_out_arrays);
+
+    /* have to add offset for this processor to coordinates array */
+    interp_local_coords += myGH->num_points_to[proc] * num_dims;
+
+  } /* end of loop over all processors */
+
+  /* don't need these anymore */
+  if (num_local_points > 0)
+  {
+    interp_local_coords -= num_local_points * num_dims;
+    free (interp_local_coords);
+  }
+  free (local_out_arrays);
+  free (in_arrays);
+  free (origin);
+  free (dims);
+
+  /* now send the interpolation results back to the processors they were
+     requested from, also receive my own results that were computed
+     by other processors
+     Since all the locally computed results are in a single contiguous buffer
+     we need to call MPI_Alltoall() only once for each CCTK data type. */
+  for (type = 0; type <= max_type; type++)
+  {
+    /* skip unused types */
+    if (type_desc[type].num_arrays <= 0)
+    {
+      continue;
+    }
+
+    /* set up the communication (this is type-independent) */
+    myGH->send_count[0] = type_desc[type].num_arrays * myGH->num_points_to[0];
+    myGH->recv_count[0] = type_desc[type].num_arrays * myGH->num_points_from[0];
+    myGH->send_displ[0] = myGH->recv_displ[0] = 0;
+    for (proc = 1; proc < pughGH->nprocs; proc++)
+    {
+      myGH->send_count[proc] = type_desc[type].num_arrays *
+                               myGH->num_points_to[proc];
+      myGH->recv_count[proc] = type_desc[type].num_arrays *
+                               myGH->num_points_from[proc];
+      myGH->send_displ[proc] = myGH->send_displ[proc-1] +
+                               myGH->send_count[proc-1];
+      myGH->recv_displ[proc] = myGH->recv_displ[proc-1] +
+                               myGH->recv_count[proc-1];
+    }
+
+    /* allocate buffer for receiving my own requested points */
+    /* avoid NULL pointers here because MPI_Alltoallv() doesn't like it */
+    if (num_points > 0)
+    {
+      type_desc[type].recv_buffer = malloc (num_points *
+                                            type_desc[type].num_arrays *
+                                            type_desc[type].vtypesize);
+    }
+    else
+    {
+      /* access to such a fake address should crash the code on most systems */
+      type_desc[type].recv_buffer = (void *) 1;
+    }
+
+    /* now exchange the data for this CCTK data type */
+    CACTUS_MPI_ERROR (MPI_Alltoallv (type_desc[type].send_buffer,
+                                     myGH->send_count, myGH->send_displ,
+                                     type_desc[type].mpi_type,
+                                     type_desc[type].recv_buffer,
+                                     myGH->recv_count, myGH->recv_displ,
+                                     type_desc[type].mpi_type,
+                                     pughGH->PUGH_COMM_WORLD));
+
+    /* now that the data is sent we don't need the buffer anymore */
+    if (num_local_points > 0)
+    {
+      free (type_desc[type].send_buffer);
+    }
+
+    /* no sort neccessary if there are no points */
+    if (num_points <= 0)
+    {
+      continue;
+    }
+
+    /* Go back from processor-sorted data to input-ordered data.
+       The creation of the indices array above makes this not so bad. */ 
+    offset = 0;
+    type_desc[type].buffer = (char *) type_desc[type].recv_buffer;
+    for (proc = 0; proc < nprocs; proc++)
+    {
+      for (array = 0; array < num_out_arrays; array++)
+      {
+        if (out_array_types[array] != type)
+        {
+          continue;
+        }
+
+        /* now do the sorting according to the CCTK data type */
+        if (out_array_types[array] == CCTK_VARIABLE_CHAR)
+        {
+          SORT_TYPED_ARRAY (CCTK_CHAR, type_desc[type].buffer,
+                            out_arrays[array], myGH->num_points_from[proc],
+                            myGH->indices, offset);
+        }
+        else if (out_array_types[array] == CCTK_VARIABLE_INT)
+        {
+          SORT_TYPED_ARRAY (CCTK_INT, type_desc[type].buffer,
+                            out_arrays[array], myGH->num_points_from[proc],
+                            myGH->indices, offset);
+        }
+        else if (out_array_types[array] == CCTK_VARIABLE_REAL)
+        {
+          SORT_TYPED_ARRAY (CCTK_REAL, type_desc[type].buffer,
+                            out_arrays[array], myGH->num_points_from[proc],
+                            myGH->indices, offset);
+        }
+        else if (out_array_types[array] == CCTK_VARIABLE_COMPLEX)
+        {
+          SORT_TYPED_ARRAY (CCTK_COMPLEX, type_desc[type].buffer,
+                            out_arrays[array], myGH->num_points_from[proc],
+                            myGH->indices, offset);
+        }
+#ifdef CCTK_REAL4
+        else if (out_array_types[array] == CCTK_VARIABLE_REAL4)
+        {
+          SORT_TYPED_ARRAY (CCTK_REAL4, type_desc[type].buffer,
+                            out_arrays[array], myGH->num_points_from[proc],
+                            myGH->indices, offset);
+        }
+        else if (out_array_types[array] == CCTK_VARIABLE_COMPLEX8)
+        {
+          SORT_TYPED_ARRAY (CCTK_COMPLEX8, type_desc[type].buffer,
+                            out_arrays[array], myGH->num_points_from[proc],
+                            myGH->indices, offset);
+        }
+#endif
+#ifdef CCTK_REAL8
+        else if (out_array_types[array] == CCTK_VARIABLE_REAL8)
+        {
+          SORT_TYPED_ARRAY (CCTK_REAL8, type_desc[type].buffer,
+                            out_arrays[array], myGH->num_points_from[proc],
+                            myGH->indices, offset);
+        }
+        else if (out_array_types[array] == CCTK_VARIABLE_COMPLEX16)
+        {
+          SORT_TYPED_ARRAY (CCTK_COMPLEX16, type_desc[type].buffer,
+                            out_arrays[array], myGH->num_points_from[proc],
+                            myGH->indices, offset);
+        }
+#endif
+#ifdef CCTK_REAL16
+        else if (out_array_types[array] == CCTK_VARIABLE_REAL16)
+        {
+          SORT_TYPED_ARRAY (CCTK_REAL16, type_desc[type].buffer,
+                            out_arrays[array], myGH->num_points_from[proc],
+                            myGH->indices, offset);
+        }
+        else if (out_array_types[array] == CCTK_VARIABLE_COMPLEX32)
+        {
+          SORT_TYPED_ARRAY (CCTK_COMPLEX32, type_desc[type].buffer,
+                            out_arrays[array], myGH->num_points_from[proc],
+                            myGH->indices, offset);
+        }
+#endif
+        else
+        {
+          CCTK_WARN (0, "Implementation error");
+        }
+
+      } /* end of loop over all output arrays */
+
+      /* advance the offset into the communication receive buffer */
+      offset += myGH->num_points_from[proc];
+
+    } /* end of loop over all processors */
+
+    /* this communication receive buffer isn't needed anymore */
+    free (type_desc[type].recv_buffer);
+
+  } /* end of loop over all types */
+
+  /* free remaining resources allocated within this run */
+  if (myGH->whichproc)
+  {
+    free (myGH->whichproc);
+    myGH->whichproc = NULL;
+  }
+  free (type_desc);
+#endif /* MPI */
+
+  return (0);
+}
+
+
+/*@@
+  @routine    PUGHInterp_InterpLocal
+  @date       Sun Jul 04 1999
+  @author     Thomas Radke
+  @desc
+              The interpolation operator registered with the CCTK
+              under the name "uniform cartesian".
+
+              Interpolates a list of non-distributed, processor-local
+              input arrays to a list of output arrays (one-to-one)
+              at a given number of interpolation points (indicated by
+              their coordinates). The points are located on a coordinate
+              system which is assumed to be a uniform cartesian.
+  @enddesc
+
+  @var        GH
+  @vdesc      Pointer to CCTK grid hierarchy
+  @vtype      cGH *
+  @vio        in
+  @endvar
+  @var        order
+  @vdesc      interpolation order
+  @vtype      int
+  @vio        in
+  @endvar
+  @var        num_points
+  @vdesc      number of points to be interpolated on this processor
+  @vtype      int
+  @vio        in
+  @endvar
+  @var        num_dims
+  @vdesc      dimensionality of the underlying grid
+  @vtype      int
+  @vio        in
+  @endvar
+  @var        num_in_arrays
+  @vdesc      number of input arrays to interpolate from
+  @vtype      int
+  @vio        in
+  @endvar
+  @var        num_out_arrays
+  @vdesc      number of output arrays to interpolate to
+  @vtype      int
+  @vio        in
+  @endvar
+  @var        coord_dims
+  @vdesc      dimensions of the underlying grid
+  @vtype      int [num_dims]
+  @vio        in
+  @endvar
+  @var        coord_arrays
+  @vdesc      list of grid coordinate arrays
+  @vtype      void *[num_dims]
+  @vio        in
+  @endvar
+  @var        coord_array_types
+  @vdesc      CCTK data type of coordinate arrays
+  @vtype      int [num_dims]
+  @vio        int
+  @endvar
+  @var        interp_coord_arrays
+  @vdesc      coordinates of points to interpolate at
+  @vtype      void *[num_dims]
+  @vio        in
+  @endvar
+  @var        interp_coord_array_types
+  @vdesc      CCTK data type of coordinate arrays to interpolate at
+  @vtype      int [num_dims]
+  @vio        out
+  @endvar
+  @var        in_arrays
+  @vdesc      list of input arrays to interpolate on
+  @vtype      void *[num_in_arrays]
+  @vio        in
+  @endvar
+  @var        in_array_types
+  @vdesc      CCTK data types of input arrays
+  @vtype      int [num_in_arrays]
+  @vio        in
+  @endvar
+  @var        out_arrays
+  @vdesc      list of output arrays to interpolate to
+  @vtype      void *[num_out_arrays]
+  @vio        out
+  @endvar
+  @var        out_array_types
+  @vdesc      CCTK data types of output arrays
+  @vtype      int [num_out_arrays]
+  @vio        in
+  @endvar
+
+  @returntype int
+  @returndesc
+              0  - successful interpolation
+              -1 - in case of any errors
+  @endreturndesc
+@@*/
+int PUGHInterp_InterpLocal (cGH *GH,
+                            int order,
+                            int num_points,
+                            int num_dims,
+                            int num_in_arrays,
+                            int num_out_arrays,
+                            int coord_dims[],
+                            void *coord_arrays[],
+                            int coord_array_types[],
+                            void *interp_coord_arrays[],
+                            int interp_coord_array_types[],
+                            void *in_arrays[],
+                            int in_array_types[],
+                            void *out_arrays[],
+                            int out_array_types[])
+{
+  int dim, point, retval;
+  CCTK_REAL *coords, *origin, *delta, **data;
+
+
+  /* check arguments */
+  retval = PUGHInterp_CheckArguments (GH, num_dims, num_points,
+                                      num_in_arrays, num_out_arrays,
+                                      interp_coord_array_types);
+  if (retval <= 0)
+  {
+    return (retval);
+  }
+  for (dim = 0; dim < num_dims; dim++)
+  {
+    if (coord_array_types[dim] != CCTK_VARIABLE_REAL)
+    {
+      CCTK_WARN (1, "Coordinates should be of type CCTK_REAL");
+      return (-1);
+    }
+    if (interp_coord_array_types[dim] != CCTK_VARIABLE_REAL)
+    {
+      CCTK_WARN (1, "Interpolation coordinates should be of type CCTK_REAL");
+      return (-1);
+    }
+  }
+
+  /* get the grid spacings - this assumes a cartesian grid */
+  origin = (CCTK_REAL *) malloc (2 * num_dims * sizeof (CCTK_REAL));
+  delta  = origin + num_dims;
+  data   = (CCTK_REAL **) coord_arrays;
+  for (dim = 0; dim < num_dims; dim++)
+  {
+    origin[dim] = data[dim][0];
+    delta[dim]  = data[dim][1] - data[dim][0];
+  }
+
+  /* sort the individual interpolation coordinate arrays into a single one */
+  coords = (CCTK_REAL *) malloc (num_dims * num_points * sizeof (CCTK_REAL));
+  data   = (CCTK_REAL **) interp_coord_arrays;
+  for (point = 0; point < num_points; point++)
+  {
+    for (dim = 0; dim < num_dims; dim++)
+    {
+      *coords++ = data[dim][point];
+    }
+  }
+  coords -= num_dims * num_points;
+
+  /* call the interpolator function */
+  retval = PUGHInterp_Interpolate (order,
+                                   num_points, num_dims, num_out_arrays,
+                                   coord_dims, coords, origin, delta,
+                                   in_array_types, in_arrays,
+                                   out_array_types, out_arrays);
+
+  /* free allocated resources */
+  free (coords);
+  free (origin);
+
+  return (retval);
+}
+
+
+/**************************************************************************/
+/*                            local routines                              */
+/**************************************************************************/
+
+#ifdef CCTK_MPI
+
+/*@@
+  @routine   GetLocalCoords
+  @date      Sun Jul 04 1999
+  @author    Thomas Radke
+  @desc 
+  Collect the coordinates of all points to be processed locally
+  into an array coords[num_dims][num_local_points].
+  <B>
+  This means for the single-processor case to sort
+
+   inCoords1[num_points], inCoords2[npoints], ..., inCoords<num_dims>[npoints]
+    into coords[num_dims][num_points]
+
+  where num_points == num_local_points.
+  <B>
+  In the multiprocessor case all processors map their points' coordinates
+  to the processor that owns this point and exchange this information
+  via MPI_Alltoall ().
+  num_local_points is then the number of all processors' points to be
+  interpolated locally on this processor.
+  <B>
+  This routine returns the number of points to be processed locally and
+  a pointer to the allocated array of their coordinates.
+  @enddesc 
+
+  @var        GH
+  @vdesc      Pointer to CCTK grid hierarchy
+  @vtype      cGH
+  @vio        in
+  @endvar
+  @var        num_points
+  @vdesc      number of points to be interpolated on this processor
+  @vtype      int
+  @vio        in
+  @endvar
+  @var        isGlobal
+  @vdesc      flag indicating that coordinates are global (and need to be
+              collected over all processors)
+  @vtype      int
+  @vio        in
+  @endvar
+  @var        num_dims
+  @vdesc      number of coordinate dimensions for each point
+  @vtype      int
+  @vio        in
+  @endvar
+  @var        coords
+  @vdesc      coordinates of each point to be interpolated on this processor
+  @vtype      CCTK_REAL array of size num_dims
+  @vio        in
+  @endvar
+  @var        num_local_points
+  @vdesc      number of points to be processed by this processor
+  @vtype      int *
+  @vio        out
+  @endvar
+  @var        local_coords
+  @vdesc      coordinates of each point to be processed by this processor
+  @vtype      pointer to CCTK_REAL array of dims[num_dims][num_local_points]
+  @vio        out
+  @endvar
+
+@@*/
+static int GetLocalCoords (cGH *GH,
+                           int num_points,
+                           const char *coord_system,
+                           pGExtras *extras,
+                           CCTK_REAL *coords[],
+                           int *num_local_points,
+                           CCTK_REAL **local_coords)
+{
+  DECLARE_CCTK_PARAMETERS
+  int dim, point;
+  int tmp, proc, nprocs;
+  pGH *pughGH;
+  pughInterpGH *myGH;
+  CCTK_REAL *range_min, *range_max;
+  CCTK_REAL *origin, *delta;
+  CCTK_REAL *proc_coords;
+#define FUDGE 0.00001
+
+
+  /* Do a fliperooni in octant mode */
+  /* This flips the sign around the origin only */
+  /* if the left side (cx0) is bigger than the position to interpolate at */
+  /* BoxInBox needs this for interpolating onto grids */ 
+  if (CCTK_Equals (domain, "octant"))
+  {
+    for (point = 0; point < num_points; point++)
+    {
+      tmp = 0;
+      for (dim = 0; dim < extras->dim; dim++)
+      {
+        if (coords[dim][point] <= GH->cctk_origin_space[dim])
+        {
+          tmp = 1;                 /* mark flipping */
+          break;
+        }
+      }
+/*** FIXME *** really change input arrays here ??? */
+      if (tmp)
+      {
+        for (dim = 0; dim < extras->dim; dim++)
+        {
+          coords[dim][point] = fabs (coords[dim][point]);
+        }
+      }
+    }
+  }
+
+  /* get GH extension handles for PUGHInterp and PUGH */
+  myGH = (pughInterpGH *) CCTK_GHExtension (GH, "PUGHInterp");
+  pughGH = (pGH *) CCTK_GHExtension (GH, "PUGH");
+
+  /* This holds the proccessor for *each* of num_points points */
+  if (num_points > 0)
+  {
+    myGH->whichproc = (int *) malloc (2 * num_points * sizeof (int));
+  }
+  else
+  {
+    myGH->whichproc = NULL;
+  }
+  /* indices[] is used to make the sorting easier
+     when receiving the output data */
+  myGH->indices = myGH->whichproc + num_points;
+
+  /* initialize whichproc with invalid processor number -1 */
+  for (point = 0; point < num_points; point++)
+  {
+    myGH->whichproc[point] = -1;
+  }
+  
+  /* initialize num_points_from to 0 for counting it up in the following loop */
+  nprocs = CCTK_nProcs (GH);
+  memset (myGH->num_points_from, 0, nprocs * sizeof (CCTK_INT));
+
+  /* allocate the ranges for my local coordinates */
+  range_min = (CCTK_REAL *) malloc (4 * extras->dim * sizeof (CCTK_REAL));
+  range_max = range_min + 1*extras->dim;
+  origin    = range_min + 2*extras->dim;
+  delta     = range_min + 3*extras->dim;
+
+  /* get the global origin and delta of the coordinate system */
+  for (dim = 0; dim < extras->dim; dim++)
+  {
+    CCTK_CoordRange (GH, &origin[dim], &delta[dim], dim+1, NULL, coord_system);
+    delta[dim] = (delta[dim] - origin[dim]) / extras->nsize[dim];
+  }
+
+  /* locate the points to interpolate at */
+  for (proc = 0; proc < nprocs; proc++)
+  {
+    for (dim = 0; dim < extras->dim; dim++)
+    {
+      /* compute the coordinate ranges */
+      range_min[dim] = origin[dim] + (extras->lb[proc][dim] - FUDGE)*delta[dim];
+      range_max[dim] = origin[dim] + (extras->ub[proc][dim] + FUDGE)*delta[dim];
+    }
+
+    /* and now which point will be processed by what processor */
+    for (point = 0; point < num_points; point++)
+    {
+      /* skip points which have already been located */
+      if (myGH->whichproc[point] >= 0)
+      {
+        continue;
+      }
+
+      /* check whether the point belongs to this processor
+         (must be within min/max in all dimensions) */
+      tmp = 0;
+      for (dim = 0; dim < extras->dim; dim++)
+      {
+        if (coords[dim][point] >= range_min[dim] &&
+            coords[dim][point] <= range_max[dim])
+        {
+          tmp++;
+        }
+      }
+      if (tmp == extras->dim)
+      {
+        myGH->whichproc[point] = proc;
+        myGH->num_points_from[proc]++;
+      }
+    }
+  }
+  /* don't need this anymore */
+  free (range_min);
+
+  /* make sure that all points could be mapped onto a processor */
+  tmp = 0;
+  for (point = 0; point < num_points; point++)
+  {
+    if (myGH->whichproc[point] < 0)
+    {
+      int i;
+      char *msg = (char *) malloc (80 + extras->dim*20);
+
+
+      sprintf (msg, "Unable to locate point %d [%f",
+               point, coords[0][point]);
+      for (i = 1; i < extras->dim; i++)
+      {
+        sprintf (msg, "%s %f", msg, coords[i][point]);
+      }
+      sprintf (msg, "%s]", msg);
+      CCTK_WARN (1, msg);
+      free (msg);
+      tmp = 1;                       /* mark as error */
+    }
+  }
+  if (tmp)
+  {
+    if (myGH->whichproc)
+    {
+      free (myGH->whichproc);
+      myGH->whichproc = NULL;
+    }
+    return (-1);
+  }
+
+  /* Now we want to resolve the num_points_from[]. Currently this is
+     the form of ( in 2 proc mode )
+     P1:  Num from P1  NFP2
+     P2:  NFP1         NFP2
+
+     and this needs to become
+     P1:  Num to P1    NTP2
+     P2:  NTP1         NTP1
+
+     Since NTP1 = NFP2 (and this works in more proc mode too)
+     this is an all-to-all communication.
+     */
+  CACTUS_MPI_ERROR (MPI_Alltoall (myGH->num_points_from, 1, PUGH_MPI_INT,
+                                  myGH->num_points_to, 1, PUGH_MPI_INT,
+                                  pughGH->PUGH_COMM_WORLD));
+
+#ifdef PUGHINTERP_DEBUG
+  for (proc = 0; proc < nprocs; proc++)
+  {
+    printf ("processor %d <-> %d  From: %d  To: %d\n",
+            CCTK_MyProc (GH), proc, myGH->num_points_from[proc],
+            myGH->num_points_to[proc]);
+  }
+#endif
+
+  /* Great. Now we know how many to expect from each processor,
+     and how many to send to each processor. So first we have
+     to send the locations to the processors which hold our data.
+     This means I send coords[dim][point] to whichproc[point].
+     I have num_points_from[proc] to send to each processor.
+     */
+  
+  /* This is backwards in the broadcast location; the number of points
+     we are getting is how many everyone else is sending to us,
+     eg, num_points_to, not how many we get back from everyone else,
+     eg, num_points_from. The number we are sending, of course, is
+     all of our locations, eg, num_points */
+  *num_local_points = 0;
+  for (proc = 0; proc < nprocs; proc++)
+  {
+    *num_local_points += myGH->num_points_to[proc];
+  }
+
+#ifdef PUGHINTERP_DEBUG
+  printf ("processor %d gets %d points in total\n",
+           CCTK_MyProc (GH), *num_local_points);
+#endif
+
+  /* allocate the local coordinates array (sorted in processor order)
+     and the resulting coordinates array that I have to process */
+  proc_coords = (CCTK_REAL *) malloc (extras->dim * num_points *
+                                      sizeof (CCTK_REAL));
+  *local_coords = (CCTK_REAL *) malloc (extras->dim * *num_local_points *
+                                        sizeof (CCTK_REAL));
+
+  /* now sort my own coordinates as tupels of [extras->dim] */
+  tmp = 0;
+  for (proc = 0; proc < nprocs; proc++)
+  {
+    for (point = 0; point < num_points; point++)
+    {
+      if (myGH->whichproc[point] == proc)
+      {
+        for (dim = 0; dim < extras->dim; dim++)
+        {
+          *proc_coords++ = coords[dim][point];
+        }
+        myGH->indices[tmp++] = point;
+      }
+    }
+  }
+  proc_coords -= tmp * extras->dim;
+
+  /* So load up the send and recv stuff */
+  /* Send extras->dim elements per data point */
+  myGH->send_count[0] = extras->dim * myGH->num_points_from[0];
+  myGH->recv_count[0] = extras->dim * myGH->num_points_to[0];
+  myGH->send_displ[0] = myGH->recv_displ[0] = 0;
+  for (proc = 1; proc < nprocs; proc++)
+  {
+    myGH->send_count[proc] = extras->dim * myGH->num_points_from[proc];
+    myGH->recv_count[proc] = extras->dim * myGH->num_points_to[proc];
+    myGH->send_displ[proc] = myGH->send_displ[proc-1] +
+                              myGH->send_count[proc-1];
+    myGH->recv_displ[proc] = myGH->recv_displ[proc-1] +
+                              myGH->recv_count[proc-1];
+  }
+ 
+  /* Great, and now exchange the coordinates and collect the ones
+     that I have to process in *local_coords[] */
+  CACTUS_MPI_ERROR (MPI_Alltoallv (proc_coords, myGH->send_count,
+                                   myGH->send_displ, PUGH_MPI_REAL,
+                                   *local_coords, myGH->recv_count,
+                                   myGH->recv_displ, PUGH_MPI_REAL,
+                                   pughGH->PUGH_COMM_WORLD));
+
+  /* don't need this anymore */
+  free (proc_coords);
+
+  return (0);
+}
+#endif /* CCTK_MPI */
+
+
+ /*@@
+   @routine    PUGHInterp_CheckArguments
+   @date       Thu 25 Jan 2001
+   @author     Thomas Radke
+   @desc
+               Checks the interpolation arguments passed in via
+               the flesh's general interpolation calling interface
+
+               This routine also verifies that the parameters meet
+               the limitations of PUGHInterp's interpolation operators.
+   @enddesc
+
+   @var        GH
+   @vdesc      Pointer to CCTK grid hierarchy
+   @vtype      cGH *
+   @vio        in
+   @endvar
+   @var        num_dims
+   @vdesc      dimensionality of the underlying grid
+   @vtype      int
+   @vio        in
+   @endvar
+   @var        num_points
+   @vdesc      number of points to interpolate at
+   @vtype      int
+   @vio        in
+   @endvar
+   @var        num_in_arrays
+   @vdesc      number of passed input arrays
+   @vtype      int
+   @vio        in
+   @endvar
+   @var        num_out_arrays
+   @vdesc      number of passed input arrays
+   @vtype      int
+   @vio        in
+   @endvar
+   @var        interp_coord_array_types
+   @vdesc      types of passed coordinates to interpolate at
+   @vtype      int [num_dims]
+   @vio        in
+   @endvar
+
+   @returntype int
+   @returndesc
+               +1 for success
+                0 for success but nothing to do
+               -1 for failure (wrong parameters passed or limitations not met)
+   @endreturndesc
+@@*/
+static int PUGHInterp_CheckArguments (cGH *GH,
+                                      int num_dims,
+                                      int num_points,
+                                      int num_in_arrays,
+                                      int num_out_arrays,
+                                      int interp_coord_array_types[])
+{
+  int i;
+
+
+  /* check for invalid arguments */
+  if (num_dims < 0 || num_points < 0 || num_in_arrays < 0 || num_out_arrays < 0)
+  {
+    return (-1);
+  }
+
+  /* check if there's anything to do at all */
+  /* NOTE: num_points can be 0 in a collective call */
+  if (num_dims == 0 || (CCTK_nProcs (GH) == 1 && num_points == 0) ||
+      num_in_arrays == 0 || num_out_arrays == 0)
+  {
+    return (0);
+  }
+
+  /* for now we can only deal with coordinates of type CCTK_REAL */
+  for (i = 0; i < num_dims; i++)
+  {
+    if (interp_coord_array_types[i] != CCTK_VARIABLE_REAL)
+    {
+      CCTK_WARN (1, "Interpolation coordinates must be of type CCTK_REAL");
+      return (-1);
+    }
+  }
+
+  /* PUGHInterp's interpolation operators compute one output array
+     per input array */
+  if (num_in_arrays != num_out_arrays)
+  {
+    CCTK_WARN (1, "Number of input arrays must match number of output arrays");
+    return (-1);
+  }
+
+  return (1);
+}
diff --git a/src/Startup.c b/src/Startup.c
new file mode 100644
index 0000000..729e712
--- /dev/null
+++ b/src/Startup.c
@@ -0,0 +1,253 @@
+ /*@@
+   @file      Startup.c
+   @date      Sun Jul 04 1999
+   @author    Thomas Radke
+   @desc
+              Startup routines for PUGHInterp
+   @enddesc
+   @history
+   @endhistory
+   @version   $Id$
+ @@*/
+
+#include <stdlib.h>
+
+#include "cctk.h"
+#include "cctk_Interp.h"
+#include "pughInterpGH.h"
+
+/* the rcs ID and its dummy function to use it */
+static char *rcsid = "$Header$";
+CCTK_FILEVERSION(CactusPUGH_PUGHInterp_Startup_c)
+
+
+/* prototypes of routines defined in this source file */
+void PUGHInterp_Startup (void);
+#ifdef CCTK_MPI
+static void *PUGHInterp_SetupGH (tFleshConfig *config,
+                                 int convergence_level,
+                                 cGH *GH);
+#endif
+
+
+#ifdef CCTK_MPI
+ /*@@
+   @routine   PUGHInterp_SetupGH
+   @date      Sun Jul 04 1999
+   @author    Thomas Radke
+   @desc
+              Allocates the GH extension structure and - for the MPI case -
+              the count/displacement buffers used in MPI_Alltoall()
+              and MPI_Alltoallv()
+   @enddesc
+   @calls     CCTK_nProcs
+
+   @returntype void *
+   @returndesc
+               pointer to the allocated GH extension structure
+   @endreturndesc
+@@*/
+static void *PUGHInterp_SetupGH (tFleshConfig *config,
+                                 int convergence_level,
+                                 cGH *GH)
+{
+  int nprocs;
+  pughInterpGH *myGH;
+
+
+  /* suppress compiler warnings about unused variables */
+  config = config;
+  convergence_level = convergence_level;
+
+  myGH = (pughInterpGH *) malloc (sizeof (pughInterpGH));
+
+  /* allocate once for all fields of same type */
+  nprocs = CCTK_nProcs (GH);
+  myGH->send_count = (int *) malloc (4 * nprocs * sizeof (int));
+  myGH->send_displ = myGH->send_count + nprocs;
+  myGH->recv_count = myGH->send_displ + nprocs;
+  myGH->recv_displ = myGH->recv_count + nprocs;
+
+  myGH->num_points_from = (CCTK_INT *) malloc (2 * nprocs * sizeof (CCTK_INT));
+  myGH->num_points_to   = myGH->num_points_from + nprocs;
+
+  return (myGH);
+}
+#endif /* CCTK_MPI */
+
+
+ /*@@
+   @routine   PUGHInterp_InterpGV_NthOrder
+   @date      Wed 14 Feb 2001
+   @author    Thomas Radke
+   @desc
+              Wrappers for the different interpolation operators 
+              registered for first/second/third order interpolation.
+              These wrappers just call the common interpolation routine
+              passing all arguments plus the interpolation order.
+   @enddesc
+
+   @returntype int
+   @returndesc
+               the return code of the common interpolation routine
+   @endreturndesc
+@@*/
+static int PUGHInterp_InterpGV_1stOrder (cGH *GH,
+                                         const char *coord_system,
+                                         int num_points,
+                                         int num_in_array_indices,
+                                         int num_out_arrays,
+                                         void *interp_coord_arrays[],
+                                         int interp_coord_array_types[],
+                                         int in_array_indices[],
+                                         void *out_arrays[],
+                                         int out_array_types[])
+{
+  return (PUGHInterp_InterpGV (GH, 1, coord_system, num_points,
+                               num_in_array_indices, num_out_arrays,
+                               interp_coord_arrays, interp_coord_array_types,
+                               in_array_indices, out_arrays, out_array_types));
+}
+                       
+
+static int PUGHInterp_InterpGV_2ndOrder (cGH *GH,
+                                         const char *coord_system,
+                                         int num_points,
+                                         int num_in_array_indices,
+                                         int num_out_arrays,
+                                         void *interp_coord_arrays[],
+                                         int interp_coord_array_types[],
+                                         int in_array_indices[],
+                                         void *out_arrays[],
+                                         int out_array_types[])
+{
+  return (PUGHInterp_InterpGV (GH, 2, coord_system, num_points,
+                               num_in_array_indices, num_out_arrays,
+                               interp_coord_arrays, interp_coord_array_types,
+                               in_array_indices, out_arrays, out_array_types));
+}
+                       
+
+static int PUGHInterp_InterpGV_3rdOrder (cGH *GH,
+                                         const char *coord_system,
+                                         int num_points,
+                                         int num_in_array_indices,
+                                         int num_out_arrays,
+                                         void *interp_coord_arrays[],
+                                         int interp_coord_array_types[],
+                                         int in_array_indices[],
+                                         void *out_arrays[],
+                                         int out_array_types[])
+{
+  return (PUGHInterp_InterpGV (GH, 3, coord_system, num_points,
+                               num_in_array_indices, num_out_arrays,
+                               interp_coord_arrays, interp_coord_array_types,
+                               in_array_indices, out_arrays, out_array_types));
+}
+                       
+
+static int PUGHInterp_InterpLocal_1stOrder (cGH *GH,
+                                            int num_points,
+                                            int num_dims,
+                                            int num_in_arrays,
+                                            int num_out_arrays,
+                                            int coord_dims[],
+                                            void *coord_arrays[],
+                                                   int coord_array_types[],
+                                            void *interp_coord_arrays[],
+                                            int interp_coord_array_types[],
+                                            void *in_arrays[],
+                                            int in_array_types[],
+                                            void *out_arrays[],
+                                            int out_array_types[])
+{
+  return (PUGHInterp_InterpLocal (GH, 1, num_points, num_dims,
+                                  num_in_arrays, num_out_arrays,
+                                  coord_dims, coord_arrays, coord_array_types,
+                                  interp_coord_arrays, interp_coord_array_types,
+                                  in_arrays, in_array_types,
+                                  out_arrays, out_array_types));
+}
+
+
+static int PUGHInterp_InterpLocal_2ndOrder (cGH *GH,
+                                            int num_points,
+                                            int num_dims,
+                                            int num_in_arrays,
+                                            int num_out_arrays,
+                                            int coord_dims[],
+                                            void *coord_arrays[],
+                                            int coord_array_types[],
+                                            void *interp_coord_arrays[],
+                                            int interp_coord_array_types[],
+                                            void *in_arrays[],
+                                            int in_array_types[],
+                                            void *out_arrays[],
+                                            int out_array_types[])
+{
+  return (PUGHInterp_InterpLocal (GH, 2, num_points, num_dims,
+                                  num_in_arrays, num_out_arrays,
+                                  coord_dims, coord_arrays, coord_array_types,
+                                  interp_coord_arrays, interp_coord_array_types,
+                                  in_arrays, in_array_types,
+                                  out_arrays, out_array_types));
+}
+
+
+static int PUGHInterp_InterpLocal_3rdOrder (cGH *GH,
+                                            int num_points,
+                                            int num_dims,
+                                            int num_in_arrays,
+                                            int num_out_arrays,
+                                            int coord_dims[],
+                                            void *coord_arrays[],
+                                            int coord_array_types[],
+                                            void *interp_coord_arrays[],
+                                            int interp_coord_array_types[],
+                                            void *in_arrays[],
+                                            int in_array_types[],
+                                            void *out_arrays[],
+                                            int out_array_types[])
+{
+  return (PUGHInterp_InterpLocal (GH, 3, num_points, num_dims,
+                                  num_in_arrays, num_out_arrays,
+                                  coord_dims, coord_arrays, coord_array_types,
+                                  interp_coord_arrays, interp_coord_array_types,
+                                  in_arrays, in_array_types,
+                                  out_arrays, out_array_types));
+}
+
+
+ /*@@
+   @routine   PUGHInterp_Startup
+   @date      Sun Jul 04 1999
+   @author    Thomas Radke
+   @desc
+              The startup registration routine for PUGHInterp.
+              Registers the PUGHInterp GH extensions setup routine
+              and the interpolation operators with the flesh.
+   @enddesc
+   @calls     CCTK_RegisterGHExtensionSetupGH
+              CCTK_InterpRegisterOperatorGV
+              CCTK_InterpRegisterOperatorLocal
+@@*/
+void PUGHInterp_Startup (void)
+{
+#ifdef CCTK_MPI
+  CCTK_RegisterGHExtensionSetupGH (CCTK_RegisterGHExtension ("PUGHInterp"),
+                                   PUGHInterp_SetupGH);
+#endif
+
+  CCTK_InterpRegisterOperatorGV (PUGHInterp_InterpGV_1stOrder,
+                                 "first-order uniform cartesian");
+  CCTK_InterpRegisterOperatorLocal (PUGHInterp_InterpLocal_1stOrder,
+                                    "first-order uniform cartesian");
+  CCTK_InterpRegisterOperatorGV (PUGHInterp_InterpGV_2ndOrder,
+                                 "second-order uniform cartesian");
+  CCTK_InterpRegisterOperatorLocal (PUGHInterp_InterpLocal_2ndOrder,
+                                    "second-order uniform cartesian");
+  CCTK_InterpRegisterOperatorGV (PUGHInterp_InterpGV_3rdOrder,
+                                 "third-order uniform cartesian");
+  CCTK_InterpRegisterOperatorLocal (PUGHInterp_InterpLocal_3rdOrder,
+                                    "third-order uniform cartesian");
+}
diff --git a/src/make.code.defn b/src/make.code.defn
new file mode 100644
index 0000000..b8d41c1
--- /dev/null
+++ b/src/make.code.defn
@@ -0,0 +1,5 @@
+# Main make.code.defn file for thorn PUGHInterp
+# $Header$
+
+# Source files in this directory
+SRCS = Startup.c Operator.c Interpolate.c
diff --git a/src/pughInterpGH.h b/src/pughInterpGH.h
new file mode 100644
index 0000000..e6eef76
--- /dev/null
+++ b/src/pughInterpGH.h
@@ -0,0 +1,71 @@
+ /*@@
+   @file      pughInterpGH.h
+   @date      Sun Jul 04 1999
+   @author    Thomas Radke
+   @desc
+              Definition of GH extensions of thorn PUGHInterp and
+              declaration of function prototypes
+   @enddesc
+   @history
+   @endhistory
+@@*/
+
+
+typedef struct
+{
+  int *send_count;             /* number of elements to be sent */
+  int *send_displ;             /* offset of each element to be sent */
+  int *recv_count;             /* number of elements to be received */
+  int *recv_displ;             /* offset of each element to be received */
+
+  CCTK_INT *num_points_to;     /* number of coordinate points to be sent */
+  CCTK_INT *num_points_from;   /* number of coordinate points to be received */
+
+  int *whichproc;              /* which processor does point i belong to */
+  int *indices;                /* indices to sort from processor-sorted into
+                                  point-sorted order */
+} pughInterpGH;
+
+
+/* prototypes of interpolation operators to be registered */
+int PUGHInterp_InterpGV (cGH *GH,
+                         int order,
+                         const char *coord_system,
+                         int num_points,
+                         int num_in_array_indices,
+                         int num_out_arrays,
+                         void *interp_coord_arrays[],
+                         int interp_coord_array_types[],
+                         int in_array_indices[],
+                         void *out_arrays[],
+                         int out_array_types[]);
+
+int PUGHInterp_InterpLocal (cGH *GH,
+                            int order,
+                            int num_points,
+                            int num_dims,
+                            int num_in_arrays,
+                            int num_out_arrays,
+                            int coord_dims[],
+                            void *coord_arrays[],
+                            int coord_array_types[],
+                            void *interp_coord_arrays[],
+                            int interp_coord_array_types[],
+                            void *in_arrays[],
+                            int in_array_types[],
+                            void *out_arrays[],
+                            int out_array_types[]);
+
+/* prototypes of routines used internally */
+int PUGHInterp_Interpolate (int order,
+                            int num_points,
+                            int num_dims,
+                            int num_arrays,
+                            int dims[ /* num_dims */ ],
+                            CCTK_REAL coord[ /* num_dims * num_points */ ],
+                            CCTK_REAL origin[ /* num_dims */ ],
+                            CCTK_REAL delta[ /* num_dims */ ],
+                            int in_types[ /* num_arrays */ ],
+                            void *in_arrays[ /* num_arrays */ ],
+                            int out_types[ /* num_arrays */ ],
+                            void *out_arrays[ /* num_arrays */ ]);