path: root/src/Operator.c

/*@@
  @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);
}