Changes to support both the new and the old interpolation API.

A local interpolator named "uniform cartesian" is now provided
under the new API.


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

7 files changed, 429 insertions, 463 deletions

diff --git a/schedule.ccl b/schedule.ccl
index 6a7bcc9..969f1b9 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -4,4 +4,4 @@
 Schedule LocalInterp_Startup at STARTUP after Driver_Startup
 {
 LANG:C
-} "register CCTK_InterpLocal() interpolation operators"
+} "register LocalInterp's interpolation operators"
diff --git a/src/Interpolate.c b/src/Interpolate.c
index 6e9af89..4dd6067 100644
--- a/src/Interpolate.c
+++ b/src/Interpolate.c
@@ -22,8 +22,8 @@
   @author    Jonathan Thornburg
   @hdesc     Move all local-interpolation code from LocalInterp to here
   @endhistory
-  
-  @version   $Header$
+
+  @version   $Id$
   @@*/
 
 #include <math.h>
@@ -224,12 +224,12 @@ CCTK_FILEVERSION(CactusBase_LocalInterp_Interpolate_c)
   @endvar
   @var        dims
   @vdesc      dimensions of the input arrays
-  @vtype      int[ num_dims ]
+  @vtype      CCTK_INT[ num_dims ]
   @vio        in
   @endvar
   @var        coord
-  @vdesc      coordinates to interpolate at
-  @vtype      CCTK_REAL coord[ num_dims * num_points ]
+  @vdesc      list of coordinates to interpolate at
+  @vtype      CCTK_REAL coord[ num_dims ][ num_points ]
   @vio        in
   @endvar
   @var        origin
@@ -244,7 +244,7 @@ CCTK_FILEVERSION(CactusBase_LocalInterp_Interpolate_c)
   @endvar
   @var        in_types
   @vdesc      CCTK variable types of input arrays
-  @vtype      int in_types[ num_arrays ]
+  @vtype      CCTK_INT in_types[ num_arrays ]
   @vio        in
   @endvar
   @var        in_arrays
@@ -254,7 +254,7 @@ CCTK_FILEVERSION(CactusBase_LocalInterp_Interpolate_c)
   @endvar
   @var        out_types
   @vdesc      CCTK variable types of output arrays
-  @vtype      int out_types[ num_arrays ]
+  @vtype      CCTK_INT out_types[ num_arrays ]
   @vio        in
   @endvar
   @var        out_arrays
@@ -273,13 +273,13 @@ int LocalInterp_Interpolate (int order,
                             int num_points,
                             int num_dims,
                             int num_arrays,
-                            const int dims[],
-                            const CCTK_REAL coord[],
+                            const CCTK_INT dims[],
+                            const CCTK_REAL *const coord[],
                             const CCTK_REAL origin[],
                             const CCTK_REAL delta[],
-                            const int in_types[],
+                            const CCTK_INT in_types[],
                             const void *const in_arrays[],
-                            const int out_types[],
+                            const CCTK_INT out_types[],
                             void *const out_arrays[])
 {
   int retval;
@@ -287,12 +287,13 @@ int LocalInterp_Interpolate (int order,
 #if 0
   int out_of_bounds;
 #endif
-  int max_dims[MAXDIM], point[MAXDIM];
+  CCTK_INT max_dims[MAXDIM], point[MAXDIM];
   CCTK_REAL delta_inv[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)
@@ -345,7 +346,7 @@ int LocalInterp_Interpolate (int order,
      (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));
+  memcpy (max_dims, dims, (num_dims - 1) * sizeof (*max_dims));
 
   /* 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 */
@@ -370,7 +371,7 @@ int LocalInterp_Interpolate (int order,
     for (i = 0; i < num_dims; i++)
     {
       /* closest grid point for stencil/molecule */
-      point[i] = floor ((coord[num_dims*n + i] - origin[i]) * delta_inv[i]
+      point[i] = floor ((coord[i][n] - origin[i]) * delta_inv[i]
                         - 0.5 * (order - 1));
 
 #if 0
@@ -402,19 +403,19 @@ int LocalInterp_Interpolate (int order,
       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_inv[i];
+      offset[i] = (coord[i][n] - below[i]) * delta_inv[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]);
-		}
-	}
+        {
+        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 */
 
 #if 0
@@ -425,7 +426,7 @@ if (n == LocalInterp_verbose_debug_n)
 
 
       /* put all information into a single message string for output */
-      msg = (char *) malloc (100 + num_dims*(10 + 4*20));
+      msg = malloc (100 + num_dims*(10 + 4*20));
       sprintf (msg,
                "Interpolation stencil/molecule out of bounds at grid point [%d",
                point[0]);
@@ -516,17 +517,17 @@ if (n == LocalInterp_verbose_debug_n)
 
 #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]);
-			}
-		}
-	}
+        {
+        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 */
@@ -539,6 +540,12 @@ if (n == LocalInterp_verbose_debug_n)
         continue;
       }
 
+      /* skip this array if it's a query call only */
+      if (! out_arrays[a])
+      {
+        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)
diff --git a/src/Interpolate.h b/src/Interpolate.h
index 35dac21..c493f92 100644
--- a/src/Interpolate.h
+++ b/src/Interpolate.h
@@ -13,33 +13,50 @@
   @endhistory
   @@*/
 
+#ifndef _LOCALINTERP_INTERPOLATE_H_
+#define _LOCALINTERP_INTERPOLATE_H_  1
+
+#ifdef __cplusplus
+extern "C"
+{
+#endif
+
 /* 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[]);
+int LocalInterp_InterpLocalUniform (int num_dims,
+                                    int param_table_handle,
+                                    /***** coordinate system *****/
+                                    const CCTK_REAL coord_origin[],
+                                    const CCTK_REAL coord_delta[],
+                                    /***** interpolation points *****/
+                                    int num_interp_points,
+                                    int interp_coords_type_code,
+                                    const void *const interp_coords[],
+                                    /***** input arrays *****/
+                                    int num_input_arrays,
+                                    const CCTK_INT input_array_dims[],
+                                    const CCTK_INT input_array_type_codes[],
+                                    const void *const input_arrays[],
+                                    /***** output arrays *****/
+                                    int num_output_arrays,
+                                    const CCTK_INT output_array_type_codes[],
+                                    void *const output_arrays[]);
 
 /* 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 */ ]);
+int LocalInterp_Interpolate (int order,
+                             int num_points,
+                             int num_dims,
+                             int num_arrays,
+                             const CCTK_INT dims[ /* num_dims */ ],
+                             const CCTK_REAL *const coords[],
+                             const CCTK_REAL origin[ /* num_dims */ ],
+                             const CCTK_REAL delta[ /* num_dims */ ],
+                             const CCTK_INT in_types[ /* num_arrays */ ],
+                             const void *const in_arrays[ /* num_arrays */ ],
+                             const CCTK_INT out_types[ /* num_arrays */ ],
+                             void *const out_arrays[ /* num_arrays */ ]);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif  /* _LOCALINTERP_INTERPOLATE_H_ */
diff --git a/src/Operator.c b/src/Operator.c
index 5d17ab8..08f7027 100644
--- a/src/Operator.c
+++ b/src/Operator.c
@@ -2,9 +2,9 @@
   @file      Operator.c
   @date      Tue Apr 15 18:22:45 1997
   @author    Paul Walker
-  @desc 
+  @desc
              Definition of interpolation operators for regular uniform grids.
-  @enddesc 
+  @enddesc
 
   @history
   @date      Sun Jul 04 1999
@@ -18,7 +18,7 @@
   @hdesc     Move all local-interpolation code from LocalInterp to here
   @endhistory
 
-  @version   $Header$
+  @version   $Id$
   @@*/
 
 #include <stdlib.h>
@@ -27,297 +27,122 @@
 
 #include "cctk.h"
 #include "cctk_Parameters.h"
-#include "CactusPUGH/PUGH/src/include/pugh.h"
+#include "util_ErrorCodes.h"
+#include "util_Table.h"
 #include "Interpolate.h"
 
 /* the rcs ID and its dummy function to use it */
 static const char *rcsid = "$Header$";
 CCTK_FILEVERSION(CactusBase_LocalInterp_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[])
+/********************************************************************
+ ********************    External Routines   ************************
+ ********************************************************************/
+int LocalInterp_InterpLocalUniform (int num_dims,
+                                    int table,
+                                    /***** coordinate system *****/
+                                    const CCTK_REAL coord_origin[],
+                                    const CCTK_REAL coord_delta[],
+                                    /***** interpolation points *****/
+                                    int num_interp_points,
+                                    int interp_coords_type_code,
+                                    const void *const interp_coords[],
+                                    /***** input arrays *****/
+                                    int num_input_arrays,
+                                    const CCTK_INT input_array_dims[],
+                                    const CCTK_INT input_array_type_codes[],
+                                    const void *const input_arrays[],
+                                    /***** output arrays *****/
+                                    int num_output_arrays,
+                                    const CCTK_INT output_array_type_codes[],
+                                    void *const output_arrays[])
 {
-  int dim, point, retval;
-  CCTK_REAL *coords, *origin, *delta;
-  const CCTK_REAL *const *data;
+  int iterator, retval;
+  char key[128];
+  CCTK_INT order, type, nelems;
 
-  /* 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++)
+  /* check for invalid arguments */
+  if (num_dims < 0 || num_interp_points < 0 ||
+      num_input_arrays < 0 || num_output_arrays < 0)
   {
-    origin[dim] = data[dim][0];
-    delta[dim]  = data[dim][1] - data[dim][0];
+    return (UTIL_ERROR_BAD_INPUT);
   }
 
-  /* 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++)
+  /* this interpolation operator computes one output array per input array */
+  if (num_input_arrays != num_output_arrays)
   {
-    for (dim = 0; dim < num_dims; dim++)
-    {
-      *coords++ = data[dim][point];
-    }
+    CCTK_WARN (1, "Number of input arrays must match number of output arrays");
+    return (UTIL_ERROR_BAD_INPUT);
   }
-  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)
+  if (interp_coords_type_code != CCTK_VARIABLE_REAL)
   {
-    return (-1);
+    CCTK_WARN (1, "Interpolation coordinates must be of type CCTK_REAL");
+    return (UTIL_ERROR_BAD_INPUT);
   }
 
   /* 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)
+  if (num_dims == 0 || num_input_arrays == 0)
   {
     return (0);
   }
 
-  /* for now we can only deal with coordinates of type CCTK_REAL */
-  for (i = 0; i < num_dims; i++)
+
+  /* get the interpolation order from the user-supplied parameter table */
+  order = 1;
+  if (table >= 0)
   {
-    if (interp_coord_array_types[i] != CCTK_VARIABLE_REAL)
+    /* loop through all table options */
+    iterator = Util_TableItCreate (table);
+    for (iterator = Util_TableItCreate (table);
+         Util_TableItQueryIsNonNull (iterator) > 0 &&
+         Util_TableItQueryKeyValueInfo (iterator, sizeof (key), key, &type,
+                                        &nelems) > 0;
+         Util_TableItAdvance (iterator))
     {
-      CCTK_WARN (1, "Interpolation coordinates must be of type CCTK_REAL");
-      return (-1);
+      if (CCTK_Equals (key, "order"))
+      {
+        if (type == CCTK_VARIABLE_INT && nelems == 1)
+        {
+          Util_TableGetInt (table, &order, "order");
+        }
+        else
+        {
+          CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+                      "Invalid value for option 'order' in interpolation "
+                      "parameter options table (must be CCTK_INT scalar value)",
+                      key);
+        }
+      }
+      else if (CCTK_Equals (key, "N_boundary_points_to_omit") ||
+               CCTK_Equals (key, "boundary_off_centering_tolerance") ||
+               CCTK_Equals (key, "boundary_extrapolation_tolerance"))
+      {
+        /* warn about unsupported options */
+        CCTK_VWarn (4, __LINE__, __FILE__, CCTK_THORNSTRING,
+                    "Option with key '%s' in interpolation parameter options "
+                    "table is not not supported (will be ignored)", key);
+      }
+      else
+      {
+        /* warn about other options */
+        CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+                    "Found option with unrecognized key '%s' in interpolation "
+                    "parameter options table (will be ignored)", key);
+      }
     }
+    Util_TableItDestroy (iterator);
   }
 
-  /* 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);
-  }
+  /* call the interpolator function */
+  retval = LocalInterp_Interpolate (order, num_interp_points, num_dims,
+                                    num_output_arrays, input_array_dims,
+                                    (const CCTK_REAL *const *) interp_coords,
+                                    coord_origin, coord_delta,
+                                    input_array_type_codes, input_arrays,
+                                    output_array_type_codes, output_arrays);
 
-  return (1);
+  return (retval);
 }
diff --git a/src/README b/src/README
index e82e759..46e71bf 100644
--- a/src/README
+++ b/src/README
@@ -4,6 +4,11 @@ under the names
    "first-order uniform cartesian"
    "second-order uniform cartesian"
    "third-order uniform cartesian"
+and for
+   CCTK_InterpLocalUniform()
+under the name
+   "uniform cartesian"
+
 
 Implementation Notes
 ====================
diff --git a/src/Startup.c b/src/Startup.c
index c51b34a..c8277f6 100644
--- a/src/Startup.c
+++ b/src/Startup.c
@@ -19,6 +19,8 @@
 
 #include "cctk.h"
 #include "cctk_Interp.h"
+#include "util_ErrorCodes.h"
+#include "util_Table.h"
 #include "Interpolate.h"
 
 /* the rcs ID and its dummy function to use it */
@@ -26,160 +28,273 @@ static const char *rcsid = "$Header$";
 CCTK_FILEVERSION(CactusBase_LocalInterp_Startup_c)
 
 
-/* prototypes of externally-visible routines defined in this source file */
+/********************************************************************
+ ***************    External Routine Prototypes   *******************
+ ********************************************************************/
 void LocalInterp_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[]);
-
-/******************************************************************************/
+/********************************************************************
+ ***************    Internal Routine Prototypes   *******************
+ ********************************************************************/
+static int 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 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 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[]);
+static int InterpLocal_NthOrder (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[]);
 
+
+/********************************************************************
+ ********************    External Routines   ************************
+ ********************************************************************/
+
+/*@@
+  @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_Startup (void)
+{
+  CCTK_InterpRegisterOpLocalUniform (LocalInterp_InterpLocalUniform,
+                                     "uniform cartesian", CCTK_THORNSTRING);
+
+  CCTK_InterpRegisterOperatorLocal (InterpLocal_1stOrder,
+                                    "first-order uniform cartesian");
+  CCTK_InterpRegisterOperatorLocal (InterpLocal_2ndOrder,
+                                    "second-order uniform cartesian");
+  CCTK_InterpRegisterOperatorLocal (InterpLocal_3rdOrder,
+                                    "third-order uniform cartesian");
+}
+
+
+/********************************************************************
+ ********************    Internal Routines   ************************
+ ********************************************************************/
 /*@@
-  @routine   LocalInterp_InterpLocal_NthOrder
+  @routine   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.
+             LocalInterp's first/second/third order interpolation operators
+             which are registered with the flesh's old interpolation API.
+             These routines are just wrappers which translate their arguments
+             to call LocalInterp's interpolation operator for the new API.
   @enddesc
 
   @returntype int
   @returndesc
-              the return code of the common interpolation routine
+              the return code of LocalInterp's interpolation operator
   @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[])
+ @@*/
+static int 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));
+  return (InterpLocal_NthOrder (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[])
+static int 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));
+  return (InterpLocal_NthOrder (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[])
+static int 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));
+  return (InterpLocal_NthOrder (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_Startup(void)
+static int InterpLocal_NthOrder (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[])
 {
-  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");
+  int i, table, retval;
+  CCTK_INT *_coord_dims, *_in_array_types, *_out_array_types;
+  CCTK_REAL *origin, *delta;
+  char string[32];
+
+
+  /* no information needed from the GH */
+  (void *) (GH + 0);
+
+  /* need to turn int's into CCTK_INT's */
+  _coord_dims = malloc ((num_dims + num_in_arrays + num_out_arrays) *
+                        sizeof (CCTK_INT));
+  _in_array_types = _coord_dims + num_dims;
+  _out_array_types = _in_array_types + num_in_arrays;
+  for (i = 0; i < num_in_arrays; i++)
+  {
+    _in_array_types[i] = in_array_types[i];
+  }
+  for (i = 0; i < num_out_arrays; i++)
+  {
+    _out_array_types[i] = out_array_types[i];
+  }
+
+  /* for now we can only deal with coordinates of type CCTK_REAL */
+  for (i = 0; i < num_dims; i++)
+  {
+    if (coord_array_types[i] != CCTK_VARIABLE_REAL)
+    {
+      CCTK_WARN (1, "Coordinates must be of type CCTK_REAL");
+      return (UTIL_ERROR_BAD_INPUT);
+    }
+    if (interp_coord_array_types[i] != CCTK_VARIABLE_REAL)
+    {
+      CCTK_WARN (1, "Interpolation coordinates must be of type CCTK_REAL");
+      return (UTIL_ERROR_BAD_INPUT);
+    }
+
+    _coord_dims[i] = coord_dims[i];
+  }
+
+  /* get the grid spacings - this assumes a uniform cartesian grid */
+  origin = malloc (2 * num_dims * sizeof (CCTK_REAL));
+  delta  = origin + num_dims;
+  for (i = 0; i < num_dims; i++)
+  {
+    origin[i] = ((const CCTK_REAL *const *) coord_arrays)[i][0];
+    delta[i]  = origin[i+1] - origin[i];
+  }
+
+  /* create a table with the interpolation order information */
+  sprintf (string, "order = %d", order);
+  table = Util_TableCreateFromString (string);
+
+  /* call the interpolator function */
+  retval = LocalInterp_InterpLocalUniform (num_dims, table, origin, delta,
+                                           num_points, CCTK_VARIABLE_REAL,
+                                           interp_coord_arrays, num_in_arrays,
+                                           _coord_dims, _in_array_types,
+                                           in_arrays, num_out_arrays,
+                                           _out_array_types, out_arrays);
+
+  /* free allocated resources */
+  Util_TableDestroy (table);
+  free (origin);
+  free (_coord_dims);
+
+  return (retval);
 }
diff --git a/src/make.code.defn b/src/make.code.defn
index 0e61525..a937952 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -3,6 +3,3 @@
 
 # Source files in this directory
 SRCS = Startup.c Operator.c Interpolate.c
-
-# Subdirectories containing source files
-SUBDIRS =