From 2a2204af7b85195002e67e3ac5218da8926e2ff8 Mon Sep 17 00:00:00 2001 From: tradke Date: Tue, 8 Jul 2003 10:15:40 +0000 Subject: 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 --- src/Interpolate.c | 77 ++++++----- src/Interpolate.h | 71 ++++++---- src/Operator.c | 357 +++++++++++++------------------------------------- src/README | 5 + src/Startup.c | 377 ++++++++++++++++++++++++++++++++++------------------- src/make.code.defn | 3 - 6 files changed, 428 insertions(+), 462 deletions(-) (limited to 'src') 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 @@ -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 @@ -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 = -- cgit v1.2.3