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