initial checkin of new LocalInterp thorn

git-svn-id: http://svn.cactuscode.org/arrangements/CactusBase/LocalInterp/trunk@2 df1f8a13-aa1d-4dd4-9681-27ded5b42416

6 files changed, 1157 insertions, 0 deletions

diff --git a/src/UniformCartesian/Interpolate.c b/src/UniformCartesian/Interpolate.c
new file mode 100644
index 0000000..95b44ce
--- /dev/null
+++ b/src/UniformCartesian/Interpolate.c
@@ -0,0 +1,586 @@
+/*@@
+  @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.  Jonathan Thornburg added some
+             additional comments in October 2001.
+  @enddesc
+
+  @history
+  @date      Wed 17 Jan 2001
+  @author    Thomas Radke
+  @hdesc     Translation from Fortran to C
+  @date      Thu 18 Oct 2001
+  @author    Jonathan Thornburg
+  @hdesc     Add lots of comments, LOCALINTERP_VERBOSE_DEBUG debugging code
+  @date      22 Jan 2002
+  @author    Jonathan Thornburg
+  @hdesc     Move all local-interpolation code from LocalInterp to here
+  @endhistory
+  
+  @version   $Header$
+  @@*/
+
+#include <math.h>
+#include <stdlib.h>
+#include <string.h>
+#include <stdio.h>
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+#include "Interpolate.h"
+
+/* the rcs ID and its dummy function to use it */
+static const char *rcsid = "$Header$";
+CCTK_FILEVERSION(CactusBase_LocalInterp_UniformCartesian_Interpolate_c)
+
+#ifdef LOCALINTERP_VERBOSE_DEBUG
+  /* if this is >= 0, we print verbose debugging information at this point */
+  int LocalInterp_verbose_debug_n = -1;
+#endif
+
+/* 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
+
+/******************************************************************************/
+
+/*@@
+  @routine INTERPOLATE (macro)
+  @date    18 Oct 2001
+  @author  code by ???, these comments by Jonathan Thornburg
+  @desc
+       This macro does the interpolation of in_array[] to compute
+       a single value  out_array[n]  (actually  out_array[n]subpart ;
+       see the comments below for details).  The data to be interpolated
+       must be real numbers of some type, i.e. if the arrays are
+       complex this macro must be called separately for the real
+       and imaginary parts.
+  @enddesc
+
+  @var    cctk_type
+  @vdesc  C type of input and output array elements (might be complex)
+  @endvar
+
+  @var    cctk_subtype
+  @vdesc  C type of actual numbers being interpolated (must be real)
+  @endvar
+
+  @var    subpart
+  @vdesc  string to be suffixed to input/output array element to get
+          to get real number, i.e. empty string if cctk_type is real,
+          .Re or .Im as appropriate if cctk_type is complex
+  @endvar
+
+  @var    in_array
+  @vdesc  A pointer to array to be interpolated (strictly speaking, to
+          the array's [0][0]...[0] element); this is typically passed
+          as a  void *  pointer so we typecast it as necessary.
+  @endvar
+
+  @var    out_array
+  @vdesc  A 1-dimensional array where the interpolation result should be
+          stored; this is typically passed as a  void *  pointer so we
+          typecast it as necessary.
+  @endvar
+
+  @var    order
+  @vdesc  The order of the interpolation (1=linear, 2=quadratic, 3=cubic, ...)
+  @endvar
+
+  @var    point
+  @vdesc  [MAXDIM] array of integers giving the integer grid coordinates
+          of the closest grid point to the interpolation point; the
+          interpolation molecule/stencil is centered at this point.
+  @endvar
+
+  @var    dims
+  @vdesc  [MAXDIM] array of integers giving the dimensions of  in_array .
+  @endvar
+
+  @var    n
+  @vdesc  Position in  out_array  where we should store the interpolation
+          result.
+  @endvar
+
+  @var    coeff
+  @vdesc  [MAXDIM][MAX_ORDER+1] array of (floating-point) interpolation
+          coefficients; detailed semantics are that coeff[axis][m] is the
+          coefficient of y[m] when the 1-dimensional Lagrange interpolation
+          polynomial passing through the  order+1  points
+             {(0,y[0]), (1,y[1]), ..., (order,y[order])}
+          is evaluated at the position x=offset[axis].
+  @endvar
+  @@*/
+/*
+ * The basic idea here is that conceptually we first interpolate the
+ * (say) 3D gridfn in the x direction at each y and z grid point,
+ * then interpolate that 2D plane of values in the y direction at
+ * each z grid point, and finally interpolate that 1D line of values
+ * in the z direction.  The implementation actually interleaves the
+ * different directions' interpolations so that only 3 scalar temporaries
+ * are needed.
+ */
+#define INTERPOLATE(cctk_type, cctk_subtype, subpart, in_array, out_array,    \
+                    order, point, dims, n, coeff)                             \
+    {                                                                         \
+      int ii, jj, kk;                                                         \
+      const 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-MAXDIM: support >3D arrays by adding more loops */              \
+      for (kk = 0; kk <= order; kk++)                                         \
+      {                                                                       \
+        fk[0] = 0;                                                            \
+        for (jj = 0; jj <= order; jj++)                                       \
+        {                                                                     \
+          /* NOTE-MAXDIM: for >3D arrays adapt the index calculation here */  \
+          fi = (const 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];                            \
+          }                                                                   \
+          /* at this point we have just computed */                           \
+          /* fj[0] = in_array[*][jj][kk] interpolated to x=offset[0] */       \
+                                                                              \
+          fk[0] += fj[0] * coeff[1][jj];                                      \
+        }                                                                     \
+        /* at this point we have just computed */                             \
+        /* fk[0] = fj[0][*][kk] interpolated to y=offset[1] */                \
+        /*       = in_array[*][*][kk] interpolated to */                      \
+        /*                            x=offset[0], y=offset[1] */             \
+                                                                              \
+        interp_result += fk[0] * coeff[2][kk];                                \
+      }                                                                       \
+      /* at this point we have just computed */                               \
+      /* interp_result = fk[0][*] interpolated to z=offset[2] */              \
+      /*               = in_array[*][*][*] interpolated to */                 \
+      /*                 x=offset[0], y=offset[1], z=offset[2] */             \
+                                                                              \
+      /* assign the result */                                                 \
+      ((cctk_type *) out_array)[n]subpart = interp_result;                    \
+    }                                                      /* end of macro */
+
+/* this is needed by some preprocessors to pass into INTERPOLATE
+   as a dummy macro */
+#define NOTHING
+
+/******************************************************************************/
+
+/*@@
+  @routine    LocalInterp_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
+              generically in that it can easily be extended to support
+              higher-dimensional arrays or more interpolation orders.
+              Places where the code would need to be changed to do this,
+              are marked with NOTE-MAXDIM and/or NOTE-MAXORDER comments
+              as appropriate.
+  @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 LocalInterp_Interpolate (int order,
+                            int num_points,
+                            int num_dims,
+                            int num_arrays,
+                            const int dims[],
+                            const CCTK_REAL coord[],
+                            const CCTK_REAL origin[],
+                            const CCTK_REAL delta[],
+                            const int in_types[],
+                            const void *const in_arrays[],
+                            const int out_types[],
+                            void *const out_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++)
+    {
+      /* closest grid point for stencil */
+      point[i] = floor ((coord[num_dims*n + i] - origin[i]) / delta[i]
+                        - 0.5 * (order - 1));
+
+      /* if beyond lower bound shift the grid point to the right */
+      shift = point[i];
+      if (shift < 0)
+      {
+        point[i] -= shift;
+      }
+
+      /* 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];
+    }
+
+#ifdef LOCALINTERP_VERBOSE_DEBUG
+if (n == LocalInterp_verbose_debug_n)
+	{
+	int ii;
+	printf("out_of_bounds = %d\n", out_of_bounds);
+		for (ii = 0 ; ii < num_dims ; ++ii)
+		{
+		printf("offset[%d] = %g\n", ii, (double) offset[ii]);
+		}
+	}
+#endif /* LOCALINTERP_VERBOSE_DEBUG */
+
+    /* 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.)
+     *
+     * These formulas are "just" the coefficients of the classical
+     * Lagrange interpolation polynomials along each dimension.
+     * For example, in 1 dimension the unique quadratic passing
+     * through the 3 points {(x0,y0), (x1,y1), (x2,y2)} is:
+     *    ( x-x1)( x-x2)        ( x-x0)( x-x2)        ( x-x0)( x-x1)
+     *    -------------- y0  +  -------------- y1  +  -------------- y2
+     *    (x0-x1)(x0-x2)        (x1-x0)(x1-x2)        (x2-x0)(x2-x1)
+     * (It's easy to see this: each of the terms is yi if x=xi, or
+     * zero if x=any other xj.)  To get the formulas below, just negate
+     * each (x-x) factor, and substitute the values xi=i.
+     */
+    /*
+     * NOTE-MAXORDER: 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");
+    }
+
+#ifdef LOCALINTERP_VERBOSE_DEBUG
+if (n == LocalInterp_verbose_debug_n)
+	{
+	int ii,mm;
+		for (ii = 0 ; ii < num_dims ; ++ii)
+		{
+			for (mm = 0 ; mm <= order ; ++mm)
+			{
+			printf("coeff[%d][%d] = %g\n",
+			       ii, mm, (double) coeff[ii][mm]);
+			}
+		}
+	}
+#endif /* LOCALINTERP_VERBOSE_DEBUG */
+
+    /* 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/UniformCartesian/Interpolate.h b/src/UniformCartesian/Interpolate.h
new file mode 100644
index 0000000..35dac21
--- /dev/null
+++ b/src/UniformCartesian/Interpolate.h
@@ -0,0 +1,45 @@
+/*@@
+  @file      Interpolate.h
+  @date      22 Jan 2002
+  @author    Jonathan Thornburg
+  @desc      function prototypes
+  @enddesc
+  @version   $Header$
+  @history
+  @date      22 Jan 2002
+  @author    Jonathan Thornburg
+  @hdesc     created by editing down LocalInterp::pughInterpGH.h
+             (originally dated 4 July 1999, by Thomas Radke)
+  @endhistory
+  @@*/
+
+/* prototypes of interpolation operator to be registered */
+int LocalInterp_InterpLocal(cGH *GH,
+                            int order,
+                            int num_points,
+                            int num_dims,
+                            int num_in_arrays,
+                            int num_out_arrays,
+                            const int coord_dims[],
+                            const void *const coord_arrays[],
+                            const int coord_array_types[],
+                            const void *const interp_coord_arrays[],
+                            const int interp_coord_array_types[],
+                            const void *const in_arrays[],
+                            const int in_array_types[],
+                            void *const out_arrays[],
+                            const int out_array_types[]);
+
+/* prototypes of routines used internally */
+int LocalInterp_Interpolate(int order,
+                            int num_points,
+                            int num_dims,
+                            int num_arrays,
+                            const int dims[ /* num_dims */ ],
+                            const CCTK_REAL coord[ /* num_dims*num_points */ ],
+                            const CCTK_REAL origin[ /* num_dims */ ],
+                            const CCTK_REAL delta[ /* num_dims */ ],
+                            const int in_types[ /* num_arrays */ ],
+                            const void *const in_arrays[ /* num_arrays */ ],
+                            const int out_types[ /* num_arrays */ ],
+                            void *const out_arrays[ /* num_arrays */ ]);
diff --git a/src/UniformCartesian/Operator.c b/src/UniformCartesian/Operator.c
new file mode 100644
index 0000000..28de748
--- /dev/null
+++ b/src/UniformCartesian/Operator.c
@@ -0,0 +1,323 @@
+/*@@
+  @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
+  @date      22 Jan 2002
+  @author    Jonathan Thornburg
+  @hdesc     Move all local-interpolation code from LocalInterp to here
+  @endhistory
+
+  @version   $Header$
+  @@*/
+
+#include <stdlib.h>
+#include <math.h>
+#include <string.h>
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+#include "CactusPUGH/PUGH/src/include/pugh.h"
+#include "Interpolate.h"
+
+/* the rcs ID and its dummy function to use it */
+static const char *rcsid = "$Header$";
+CCTK_FILEVERSION(CactusBase_LocalInterp_UniformCartesian_Operator_c)
+
+/* uncomment this to get some debugging output */
+/* #define PUGHINTERP_DEBUG 1 */
+
+
+
+/* prototypes of routines defined in this source file */
+static int LocalInterp_CheckArguments (cGH *GH,
+                                      int num_dims,
+                                      int num_points,
+                                      int num_in_arrays,
+                                      int num_out_arrays,
+                                      const int interp_coord_array_types[]);
+
+/******************************************************************************/
+
+/*@@
+  @routine    LocalInterp_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 *const [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 *const [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 *const [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 *const [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 LocalInterp_InterpLocal (cGH *GH,
+                            int order,
+                            int num_points,
+                            int num_dims,
+                            int num_in_arrays,
+                            int num_out_arrays,
+                            const int coord_dims[],
+                            const void *const coord_arrays[],
+                            const int coord_array_types[],
+                            const void *const interp_coord_arrays[],
+                            const int interp_coord_array_types[],
+                            const void *const in_arrays[],
+                            const int in_array_types[],
+                            void *const out_arrays[],
+                            const int out_array_types[])
+{
+  int dim, point, retval;
+  CCTK_REAL *coords, *origin, *delta;
+  const CCTK_REAL *const *data;
+
+  /* check arguments */
+  retval = LocalInterp_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   = (const CCTK_REAL *const *) 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   = (const CCTK_REAL *const *) 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 = LocalInterp_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                              */
+/**************************************************************************/
+
+/*@@
+  @routine    LocalInterp_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 LocalInterp'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 LocalInterp_CheckArguments (cGH *GH,
+                                      int num_dims,
+                                      int num_points,
+                                      int num_in_arrays,
+                                      int num_out_arrays,
+                                      const 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 || 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);
+    }
+  }
+
+  /* LocalInterp'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/UniformCartesian/README b/src/UniformCartesian/README
new file mode 100644
index 0000000..2342312
--- /dev/null
+++ b/src/UniformCartesian/README
@@ -0,0 +1,14 @@
+This directory contains the interpolation functions registered for
+   CCTK_InterpLocal()
+under the names
+   "first-order uniform cartesian"
+   "second-order uniform cartesian"
+   "third-order uniform cartesian"
+
+Implementation Notes
+====================
+The interpolation operators registered for different orders are mapped
+via wrappers (in "Startup.c") onto a single routine (in "Operator.c"),
+just passing the order as an additional argument.
+
+The actual core interpolation routine is located in "Interpolate.c".
diff --git a/src/UniformCartesian/Startup.c b/src/UniformCartesian/Startup.c
new file mode 100644
index 0000000..b25dc9f
--- /dev/null
+++ b/src/UniformCartesian/Startup.c
@@ -0,0 +1,185 @@
+/*@@
+  @file      Startup.c
+  @date      Sun Jul 04 1999
+  @author    Thomas Radke
+  @desc
+             Startup routines for LocalInterp/UniformCartesian
+  @enddesc
+
+  @history
+  @date      22 Jan 2002
+  @author    Jonathan Thornburg
+  @hdesc     Move all local-interpolation code from LocalInterp to here
+  @endhistory
+
+  @version   $Header$
+  @@*/
+
+#include <stdlib.h>
+
+#include "cctk.h"
+#include "cctk_Interp.h"
+#include "Interpolate.h"
+
+/* the rcs ID and its dummy function to use it */
+static const char *rcsid = "$Header$";
+CCTK_FILEVERSION(CactusBase_LocalInterp_UniformCartesian_Startup_c)
+
+
+/* prototypes of externally-visible routines defined in this source file */
+void LocalInterp_UC_Startup(void);
+
+/* prototypes of static routines defined in this source file */
+static int LocalInterp_InterpLocal_1stOrder(cGH *GH,
+                                            int num_points,
+                                            int num_dims,
+                                            int num_in_arrays,
+                                            int num_out_arrays,
+                                            const int coord_dims[],
+                                            const void *const coord_arrays[],
+                                            const int coord_array_types[],
+                                            const void *const interp_coord_arrays[],
+                                            const int interp_coord_array_types[],
+                                            const void *const in_arrays[],
+                                            const int in_array_types[],
+                                            void *const out_arrays[],
+                                            const int out_array_types[]);
+static int LocalInterp_InterpLocal_2ndOrder(cGH *GH,
+                                            int num_points,
+                                            int num_dims,
+                                            int num_in_arrays,
+                                            int num_out_arrays,
+                                            const int coord_dims[],
+                                            const void *const coord_arrays[],
+                                            const int coord_array_types[],
+                                            const void *const interp_coord_arrays[],
+                                            const int interp_coord_array_types[],
+                                            const void *const in_arrays[],
+                                            const int in_array_types[],
+                                            void *const out_arrays[],
+                                            const int out_array_types[]);
+static int LocalInterp_InterpLocal_3rdOrder(cGH *GH,
+                                            int num_points,
+                                            int num_dims,
+                                            int num_in_arrays,
+                                            int num_out_arrays,
+                                            const int coord_dims[],
+                                            const void *const coord_arrays[],
+                                            const int coord_array_types[],
+                                            const void *const interp_coord_arrays[],
+                                            const int interp_coord_array_types[],
+                                            const void *const in_arrays[],
+                                            const int in_array_types[],
+                                            void *const out_arrays[],
+                                            const int out_array_types[]);
+
+/******************************************************************************/
+
+/*@@
+  @routine   LocalInterp_InterpLocal_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 LocalInterp_InterpLocal_1stOrder(cGH *GH,
+                                            int num_points,
+                                            int num_dims,
+                                            int num_in_arrays,
+                                            int num_out_arrays,
+                                            const int coord_dims[],
+                                            const void *const coord_arrays[],
+                                            const int coord_array_types[],
+                                            const void *const interp_coord_arrays[],
+                                            const int interp_coord_array_types[],
+                                            const void *const in_arrays[],
+                                            const int in_array_types[],
+                                            void *const out_arrays[],
+                                            const int out_array_types[])
+{
+  return (LocalInterp_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 LocalInterp_InterpLocal_2ndOrder(cGH *GH,
+                                            int num_points,
+                                            int num_dims,
+                                            int num_in_arrays,
+                                            int num_out_arrays,
+                                            const int coord_dims[],
+                                            const void *const coord_arrays[],
+                                            const int coord_array_types[],
+                                            const void *const interp_coord_arrays[],
+                                            const int interp_coord_array_types[],
+                                            const void *const in_arrays[],
+                                            const int in_array_types[],
+                                            void *const out_arrays[],
+                                            const int out_array_types[])
+{
+  return (LocalInterp_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 LocalInterp_InterpLocal_3rdOrder(cGH *GH,
+                                            int num_points,
+                                            int num_dims,
+                                            int num_in_arrays,
+                                            int num_out_arrays,
+                                            const int coord_dims[],
+                                            const void *const coord_arrays[],
+                                            const int coord_array_types[],
+                                            const void *const interp_coord_arrays[],
+                                            const int interp_coord_array_types[],
+                                            const void *const in_arrays[],
+                                            const int in_array_types[],
+                                            void *const out_arrays[],
+                                            const int out_array_types[])
+{
+  return (LocalInterp_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   LocalInterp_Startup
+  @date      Sun Jul 04 1999
+  @author    Thomas Radke
+  @desc
+             The startup registration routine for LocalInterp.
+             Registers the interpolation operators with the flesh.
+  @enddesc
+  @calls     CCTK_InterpRegisterOperatorLocal
+  @@*/
+void LocalInterp_UC_Startup(void)
+{
+  CCTK_InterpRegisterOperatorLocal (LocalInterp_InterpLocal_1stOrder,
+                                    "first-order uniform cartesian");
+  CCTK_InterpRegisterOperatorLocal (LocalInterp_InterpLocal_2ndOrder,
+                                    "second-order uniform cartesian");
+  CCTK_InterpRegisterOperatorLocal (LocalInterp_InterpLocal_3rdOrder,
+                                    "third-order uniform cartesian");
+}
diff --git a/src/UniformCartesian/make.code.defn b/src/UniformCartesian/make.code.defn
new file mode 100644
index 0000000..4a506e8
--- /dev/null
+++ b/src/UniformCartesian/make.code.defn
@@ -0,0 +1,4 @@
+# $Header$
+
+# Source files in this directory
+SRCS = Startup.c Operator.c Interpolate.c