Remove Jonathan Thornburg's interpolator

	(which is GPL and thus not allowed to be in CactusBase;
	it now lives in AEIThorns/AEILocalInterp/)
so this thorn contains only the interpolator written by Thomas Radke
in early 2001.

The files for this interpolator are now in  src/ ; before this commit
they were in  src/UniformCartesian/ .


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

1 files changed, 323 insertions, 0 deletions

diff --git a/src/Operator.c b/src/Operator.c
new file mode 100644
index 0000000..5d17ab8
--- /dev/null
+++ b/src/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_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);
+}