/*@@ @file Interp.c @date July 07 1999 @author Thomas Radke @desc Registration and invocation routines for interpolation operators. @enddesc @history @date July 07 1999 @author Thomas Radke @hdesc Just copied from Reduction.c @date Wed Feb 20 17:52:55 CET 2002 @author Jonathan Thornburg @hdesc * revise registration routines to pull out common code into macro and new static function * permute arguments to the CCTKi_* registration functions * add CCTK_InterpRegisterOpLocalUniform() to register new-API interpolators @endhistory @version $Id$ @@*/ /* #define DEBUG_INTERP 1 */ #include #include #include #include #include "cctk_Constants.h" #include "cctk_Flesh.h" #include "cctk_FortranString.h" #include "cctk_Groups.h" #include "cctk_GroupsOnGH.h" #include "StoreHandledData.h" #include "cctk_Interp.h" #include "cctk_Comm.h" #include "cctk_WarnLevel.h" #include "cctk_Coord.h" #include "cctk_ActiveThorns.h" #include "util_ErrorCodes.h" static const char *rcsid = "$Header$"; CCTK_FILEVERSION(comm_Interp_c); /****************************************************************************** ************************* External Routines ***************************** ******************************************************************************/ /* prototypes for external C routines are declared in header cctk_Interp.h here only follow the fortran wrapper prototypes */ void CCTK_FCALL CCTK_FNAME (CCTK_InterpHandle) (int *handle, ONE_FORTSTRING_ARG); void CCTK_FCALL CCTK_FNAME (CCTK_InterpGV) (int *fortranreturn, cGH *GH, const int *operator_handle, const int *coord_system_handle, const int *num_points, const int *num_in_array_indices, const int *num_out_arrays, ...); void CCTK_FCALL CCTK_FNAME (CCTK_InterpLocal) (int *fortranreturn, cGH *GH, const int *operator_handle, const int *num_points, const int *num_dims, const int *num_in_arrays, const int *num_out_arrays, ...); void CCTK_FCALL CCTK_FNAME (CCTK_InterpGridArrays) (int *ierror, const cGH *GH, const int *N_dims, const int *param_table_handle, const int *local_interp_handle, const int *coord_system_handle, const int *N_interp_points, const int *interp_coords_type, const void *const interp_coords[], const int *N_input_arrays, const CCTK_INT input_array_indices[], const int *N_output_arrays, const CCTK_INT output_array_types[], void *const output_arrays[]); /****************************************************************************** ************************* Internal Data Structures ************************ ******************************************************************************/ /* structure holding the routines for a registered interpolation operator */ struct interp_info { const char *thorn_name; const char *implementation_name; const char *operator_name; cInterpOperatorGV interp_operator_GV; cInterpOperatorLocal interp_operator_local; cInterpOpLocalUniform interp_op_local_uniform; }; /****************************************************************************** ************************* Static Variables ****************************** ******************************************************************************/ /* static data: interpolation operator database and counter for registered operators */ static cHandledData *interp_operators = NULL; static int num_interp_operators = 0; /****************************************************************************** ****************** Prototypes for Functions Local to this File *************** ******************************************************************************/ static int GetOrSetupInterpInfo(const char *thorn_name, const char *operator_name, struct interp_info **pp_interp_info); /****************************************************************************** ************************* Macro Definitions ***************************** ******************************************************************************/ /* macro to read a list of items from a variable argument list into an array */ #define VARARGS_TO_ARRAY(array, type, modifier, count, varargs_list) \ { \ int i; \ \ \ for (i = 0; i < (count); i++) \ { \ (array)[i] = modifier va_arg (varargs_list, type); \ } \ } /* empty define to pass into the VARARGS_TO_ARRAY macro (some preprocessors need that) */ #define NOTHING /****************************************************************************** ************************* Registration Functions ***************************** ******************************************************************************/ /*@@ @routine CCTK_NumInterpOperators @date Mon Oct 22 2001 @author Gabrielle Allen @desc The number of interp operators registered @enddesc @returntype int @returndesc number of interpolation operators @endreturndesc @@*/ int CCTK_NumInterpOperators() { return num_interp_operators; } /******************************************************************************/ /*@@ @routine CCTK_InterpOperatorImplementation @date Mon Oct 22 2001 @author Gabrielle Allen @desc Provide the implementation which provides an interpolation operator @enddesc @returntype int @returndesc Implementation which supplied the interpolation operator @endreturndesc @@*/ const char *CCTK_InterpOperatorImplementation(int handle) { struct interp_info *operator; const char *imp=NULL; operator = (struct interp_info *) Util_GetHandledData (interp_operators, handle); if (operator) { imp = operator->implementation_name; } return imp; } /******************************************************************************/ /*@@ @routine CCTK_INTERP_REGISTER_FN_BODY (macro) @date Fri Feb 22 15:45:42 CET 2002 @author Jonathan Thornburg @desc This macro expands into the *body* (not including the opening/closing { }) of a function which registers a user-specified function as an interpolation operator. @enddesc @var operator_ptr_fnarg @vdesc the "function pointer pointing to the interpolation operator" argument to the registration function @vtype cInterpOperatorGV or cInterpOperatorLocal or cInterpOpLocalUniform or ... as appropriate @vio in @endvar @var operator_name_fnarg @vdesc the "character-string name identifying the interpolation operator" argument to the registration function @vtype const char * @vio in @endvar @var thorn_name_fnarg @vdesc the "character-string name identifying which thorn provides the operator being registered" argument to the registration function @vtype const char * @vio in @endvar @var function_name_string @vdesc the name of the registration function, as a literal C string (encased in "double quotes") @vtype const char * @vio in @endvar @var struct_field_name @vdesc the name of the field in struct interp_info which holds the function pointer pointing to the interpolation operator @vtype const char * @vio in @endvar @returntype int @returndesc the registration function returns... the handle for the newly registered operator, or

-1 NULL pointer was passed as interpolation operator routine

-2 failed to allocate memory

-3 interpolation operator by given name already exists @endreturndesc @history @date Mon 12 Feb 2001 @author Thomas Radke @hdesc Original version @date Thu Feb 21 16:03:25 CET 2002 @author Jonathan Thornburg @hdesc * move common logic in all interpolator-registration functions into new GetOrSetupInterpInfo() function * convert remaining boilerplate which differs from one registration function to another, into this macro @endhistory @@*/ #define CCTK_INTERP_REGISTER_FN_BODY(operator_ptr_fnarg, \ operator_name_fnarg, \ thorn_name_fnarg, \ function_name_string, \ struct_field_name) \ /* begin function body */ \ int handle; \ struct interp_info *p_interp_info; \ \ if (operator_ptr_fnarg == NULL) \ { \ CCTK_VWarn(0, __LINE__, __FILE__, "Cactus", \ function_name_string ":\n" \ " (called from thorn %s)\n" \ " NULL function pointer passed\n" \ " for interpolation operator \"%s\"!" \ , \ thorn_name_fnarg, \ operator_name_fnarg); \ return -1; /*** ERROR RETURN ***/ \ } \ \ handle = GetOrSetupInterpInfo(thorn_name_fnarg, \ operator_name_fnarg, \ &p_interp_info); \ if (handle < 0) \ { \ return handle; /*** ERROR RETURN ***/ \ } \ \ /* check that the operator isn't already registered */ \ if (p_interp_info->struct_field_name != NULL) \ { \ CCTK_VWarn(0, __LINE__, __FILE__, "Cactus", \ function_name_string ":\n" \ " (called from thorn %s)\n" \ " Ignoring attempt to register operator \"%s\"\n" \ " because it has already been registered by thorn %s" \ , \ thorn_name_fnarg, \ operator_name_fnarg, \ p_interp_info->thorn_name); \ return -3; /*** ERROR RETURN ***/ \ } \ \ /* now, (finally!) the actual registration */ \ p_interp_info->struct_field_name = operator_ptr_fnarg; \ \ return handle; /*** NORMAL RETURN ***/ \ /* end function body */ /* end macro */ /******************************************************************************/ /*@@ @routine CCTKi_InterpRegisterOperatorGV @routine CCTKi_InterpRegisterOperatorLocal @routine CCTK_InterpRegisterOpLocalUniform @date Thu Feb 21 16:02:05 CET 2002 @author Jonathan Thornburg @desc Each of these functions registers a user-specified function as the corresponding type of interpolation operator. The function body is produced by the CCTK_INTERP_REGISTER_FN_BODY() macro defined above. @enddesc @var operator_ptr @vdesc function pointer pointing to the interpolation operator @vtype cInterpOperatorLocal or cInterpOpLocalUniform as appropriate @vio in @endvar @var operator_name @vdesc character-string name identifying the interpolation operator @vtype const char * @vio in @endvar @var thorn_name @vdesc character-string name identifying which thorn provides the operator being registered @vtype const char * @vio in @endvar @returntype int @returndesc the handle for the newly registered operator, or

-1 NULL pointer was passed as interpolation operator routine

-2 failed to allocate memory

-3 interpolation operator by given name already exists @endreturndesc @history @date Mon 12 Feb 2001 @author Thomas Radke @hdesc Original version @date Thu Feb 21 16:03:25 CET 2002 @author Jonathan Thornburg @hdesc * move common logic in all interpolator-registration functions into new GetOrSetupInterpInfo() function * convert remaining boilerplate which differs from one registration function to another, into this macro @endhistory @@*/ /**************************************/ int CCTKi_InterpRegisterOperatorGV(cInterpOperatorGV operator_ptr, const char *operator_name, const char *thorn_name) { CCTK_INTERP_REGISTER_FN_BODY(operator_ptr, operator_name, thorn_name, "CCTKi_InterpRegisterOperatorGV", interp_operator_GV) /* no semicolon here! */ } /**************************************/ int CCTKi_InterpRegisterOperatorLocal(cInterpOperatorLocal operator_ptr, const char *operator_name, const char *thorn_name) { CCTK_INTERP_REGISTER_FN_BODY(operator_ptr, operator_name, thorn_name, "CCTKi_InterpRegisterOperatorLocal", interp_operator_local) /* no semicolon here! */ } /**************************************/ int CCTK_InterpRegisterOpLocalUniform(cInterpOpLocalUniform operator_ptr, const char *operator_name, const char *thorn_name) { CCTK_INTERP_REGISTER_FN_BODY(operator_ptr, operator_name, thorn_name, "CCTK_InterpRegisterOpLocalUniform", interp_op_local_uniform) /* no semicolon here! */ } /******************************************************************************/ /*@@ @routine GetOrSetupInterpInfo @date Thu Feb 21 14:41:35 CET 2002 @author Jonathan Thornburg @desc This is an internal worker routine used as part of the process of registering an interpolation operator. It gets the interpolator handle and the struct interp_info in which the operator information will be stored. If some interpolation operator is already registered under the specified operator name, we use the existing struct interp_info . Otherwise, we allocate a new struct interp_info and set it up (thorn, implementation, and operator names assigned from this function's arguments, all operator pointers set to NULL), then use it. @enddesc @var thorn_name @vdesc the character-string name of the thorn doing the registration @vtype const char *thorn_name @vio in @endvar @var operator_name @vdesc name identifying the interpolation operator @vtype const char *operator_name @vio in @endvar @var pp_interp_info @vdesc pointer to a struct interp_info * pointer which this function will set to point to the appropriate struct interp_info @vtype struct interp_info **pp_interp_info @vio out @endvar @returntype int @returndesc the handle for the interpolation operator, or

-2 failed to allocate memory (in this case *pp_interp_info will be set to NULL)

@endreturndesc @@*/ static int GetOrSetupInterpInfo(const char *thorn_name, const char *operator_name, struct interp_info **pp_interp_info) { /* has some operator already been registered under this operator name? */ int handle = Util_GetHandle(interp_operators, operator_name, (void **) pp_interp_info); /* yes ==> use the existing handle */ if (handle >= 0) { return handle; /*** NORMAL RETURN ***/ } /* no ==> set up a new struct interp_info for the registration */ *pp_interp_info = (struct interp_info *) malloc(sizeof(struct interp_info)); if (*pp_interp_info == NULL) { return -2; /*** ERROR RETURN ***/ } (*pp_interp_info)->thorn_name = thorn_name; (*pp_interp_info)->implementation_name = CCTK_ThornImplementation(thorn_name); (*pp_interp_info)->operator_name = operator_name; (*pp_interp_info)->interp_operator_GV = NULL; (*pp_interp_info)->interp_operator_local = NULL; (*pp_interp_info)->interp_op_local_uniform = NULL; handle = Util_NewHandle(&interp_operators, operator_name, (void *) *pp_interp_info); num_interp_operators++; return handle; /*** NORMAL RETURN ***/ } /****************************************************************************** ************* User Functions to Get Interpolator Handle/Name/etc ************* ******************************************************************************/ /*@@ @routine CCTK_InterpHandle @date July 07 1999 @author Thomas Radke @desc Returns the handle of a given interpolation operator @enddesc @var name @vdesc String containing name of interpolation operator @vtype const char * @vio in @vcomment @endvar @returntype int @returndesc the handle for the operator registered by this name or negative otherwise @endreturndesc @@*/ int CCTK_InterpHandle (const char *name) { int handle; handle = Util_GetHandle (interp_operators, name, NULL); #ifdef DEBUG_INTERP printf("In CCTK_InterpHandle\n"); printf("--------------------------\n"); printf(" Got handle %d for %s\n",handle,name); #endif if (handle < 0) { CCTK_VWarn (1, __LINE__, __FILE__, "Cactus", "No handle found for interpolation operator '%s'", name); } return (handle); } /******************************************************************************/ void CCTK_FCALL CCTK_FNAME (CCTK_InterpHandle) (int *handle, ONE_FORTSTRING_ARG) { ONE_FORTSTRING_CREATE (name) *handle = CCTK_InterpHandle (name); free (name); } /******************************************************************************/ /*@@ @routine CCTK_InterpOperator @date December 27 2001 @author Gabrielle Allen @desc Returns the name of a interpolation operator @enddesc @var handle @vdesc Handle for interpolation operator @vtype int @vio in @vcomment @endvar @returntype const char * @returndesc The name of the interpolation operator, or NULL if the handle is invalid @endreturndesc @@*/ const char *CCTK_InterpOperator (int handle) { const char *name=NULL; struct interp_info *operator; if (handle < 0) { CCTK_VWarn (6, __LINE__, __FILE__, "Cactus", "CCTK_InterpHandle: Handle %d invalid", handle); } else { operator = (struct interp_info *) Util_GetHandledData (interp_operators, handle); if (operator) { name = operator->operator_name; } else { CCTK_VWarn (6, __LINE__, __FILE__, "Cactus", "CCTK_InterpHandle: Handle %d invalid", handle); } } return name; } /****************************************************************************** ****************** User Functions to Do Interpolation ************************ ******************************************************************************/ /*@@ @routine CCTK_InterpGV @date Sun 28 Jan 2001 @author Thomas Radke @desc The general CCTK interpolation routine for grid variables Just puts the arguments from the variable argument list into arrays and calls the appropriate interpolation operator. @enddesc @var GH @vdesc pointer to CCTK grid hierarchy @vtype cGH * @vio in @endvar @var operator_handle @vdesc handle for the interpolation operator @vtype int @vio in @endvar @var coord_system_handle @vdesc handle for the coordinate system @vtype int @vio in @endvar @var num_points @vdesc number of points to interpolate at @vtype int @vio in @endvar @var num_in_array_indices @vdesc number of passed input array indices @vtype int @vio in @endvar @var num_out_arrays @vdesc number of passed output arrays @vtype int @vio in @vcomment end of fixed argument list @endvar @var interp_coord_arrays @vdesc list of coordinate arrays for points to interpolate at @vtype void [dimensions of coordinate system][num_points] @vio in @endvar @var interp_coord_array_types @vdesc types of passed coordinate arrays @vtype int [dimensions of coordinate system] @vio in @endvar @var in_array_indices @vdesc list of grid variables (given by their indices) to interpolate @vtype int [num_in_array_indices] @vio in @endvar @var out_arrays @vdesc list of output arrays which receive the interpolation results @vtype void * [num_out_arrays] @vio out @endvar @var out_array_types @vdesc types of output arrays which receive the interpolation results @vtype int [num_out_arrays] @vio in @endvar @returntype int @returndesc return value of the interpolation operator, -1 if invalid interpolation operator handle was passed in -2 if invalid coordinate system handle was passed in @endreturndesc @@*/ int CCTK_InterpGV (cGH *GH, int operator_handle, int coord_system_handle, int num_points, int num_in_array_indices, int num_out_arrays, ...) { va_list indices; int num_dims, retcode; const char *coord_system; int *in_array_indices, *interp_coord_array_types, *out_array_types; const void **interp_coord_arrays; void **out_arrays; struct interp_info *operator; /* Get the interpolation operator routine and the coordinate system name */ operator = (struct interp_info *) Util_GetHandledData (interp_operators, operator_handle); coord_system = CCTK_CoordSystemName (coord_system_handle); if ((operator == NULL) || (operator->interp_operator_GV == NULL)) { CCTK_VWarn(0, __LINE__, __FILE__, "Cactus", "\n" " CCTK_InterpGV(): no interpolation operator is registered\n" " under the handle %d\n" " (did you activate PUGHInterp or some other thorn\n" " providing this interpolation operator?)" , operator_handle); retcode = -1; } else if (coord_system == NULL) { CCTK_Warn (3, __LINE__, __FILE__, "Cactus", "CCTK_InterpGV: Invalid coordinate system handle passed to CCTK_InterpGV"); retcode = -2; } else { num_dims = CCTK_CoordSystemDim (coord_system); interp_coord_arrays = (const void **) malloc (num_dims * sizeof (void *)); interp_coord_array_types = (int *) malloc (num_dims * sizeof (int)); in_array_indices = (int *) malloc (num_in_array_indices * sizeof (int)); out_arrays = (void **) malloc (num_out_arrays * sizeof (void *)); out_array_types = (int *) malloc (num_out_arrays * sizeof (int)); /* Fill in the arrays from the variable argument list */ va_start (indices, num_out_arrays); VARARGS_TO_ARRAY (interp_coord_arrays, void *, NOTHING, num_dims, indices); VARARGS_TO_ARRAY (interp_coord_array_types, int, NOTHING, num_dims,indices); VARARGS_TO_ARRAY (in_array_indices, int, NOTHING, num_in_array_indices, indices); VARARGS_TO_ARRAY (out_arrays, void *, NOTHING, num_out_arrays, indices); VARARGS_TO_ARRAY (out_array_types, int, NOTHING, num_out_arrays, indices); va_end (indices); retcode = operator->interp_operator_GV (GH, coord_system, num_points, num_in_array_indices,num_out_arrays, interp_coord_arrays, interp_coord_array_types, in_array_indices, out_arrays, out_array_types); free (out_array_types); free (out_arrays); free (in_array_indices); free (interp_coord_array_types); free (interp_coord_arrays); } return (retcode); } /******************************************************************************/ void CCTK_FCALL CCTK_FNAME (CCTK_InterpGV) (int *fortranreturn, cGH *GH, const int *operator_handle, const int *coord_system_handle, const int *num_points, const int *num_in_array_indices, const int *num_out_arrays, ...) { va_list indices; int num_dims, retcode; const char *coord_system; int *in_array_indices, *interp_coord_array_types, *out_array_types; const void **interp_coord_arrays; void **out_arrays; struct interp_info *operator; /* Get the interpolation operator and the coordinate system name */ operator = (struct interp_info *) Util_GetHandledData (interp_operators, *operator_handle); coord_system = CCTK_CoordSystemName (*coord_system_handle); if ((operator == NULL) || (operator->interp_operator_GV == NULL)) { CCTK_VWarn(0, __LINE__, __FILE__, "Cactus", "\n" " CCTK_InterpGV(): no interpolation operator is registered\n" " under the handle %d\n" " (did you activate PUGHInterp or some other thorn\n" " providing this interpolation operator?)" , *operator_handle); retcode = -1; } else if (coord_system == NULL) { CCTK_Warn (0, __LINE__, __FILE__, "Cactus", "CCTK_InterpGV: Invalid coordinate system handle passed to CCTK_InterpGV"); retcode = -1; } else { num_dims = CCTK_CoordSystemDim (coord_system); interp_coord_arrays = (const void **) malloc (num_dims * sizeof (void *)); interp_coord_array_types = (int *) malloc (num_dims * sizeof (int)); in_array_indices = (int *) malloc (*num_in_array_indices * sizeof (int)); out_arrays = (void **) malloc (*num_out_arrays * sizeof (void *)); out_array_types = (int *) malloc (*num_out_arrays * sizeof (int)); /* Fill in the arrays from the variable argument list */ va_start (indices, num_out_arrays); VARARGS_TO_ARRAY (interp_coord_arrays, void *, NOTHING, num_dims, indices); VARARGS_TO_ARRAY (interp_coord_array_types, int *, *, num_dims, indices); VARARGS_TO_ARRAY (in_array_indices, int *, *,*num_in_array_indices,indices); VARARGS_TO_ARRAY (out_arrays, void *, NOTHING, *num_out_arrays, indices); VARARGS_TO_ARRAY (out_array_types, int *, *, *num_out_arrays, indices); va_end (indices); retcode = operator->interp_operator_GV (GH, coord_system, *num_points, *num_in_array_indices, *num_out_arrays, interp_coord_arrays, interp_coord_array_types, in_array_indices, out_arrays, out_array_types); free (out_array_types); free (out_arrays); free (in_array_indices); free (interp_coord_array_types); free (interp_coord_arrays); } *fortranreturn = retcode; } /******************************************************************************/ /*@@ @routine CCTK_InterpLocal @date Thu 01 Feb 2001 @author Thomas Radke @desc The general CCTK interpolation routine for a list of local arrays. Just puts the arguments from the variable argument list into arrays and calls the appropriate interpolation operator. @enddesc @var GH @vdesc pointer to CCTK grid hierarchy @vtype cGH * @vio in @endvar @var operator_handle @vdesc handle for the interpolation operator @vtype int @vio in @endvar @var num_points @vdesc number of points to interpolate at @vtype int @vio in @endvar @var num_dims @vdesc number of dimensions of the underlying grid @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 output arrays @vtype int @vio in @vcomment end of fixed argument list @endvar @var dims @vdesc coordinate dimensions of the underlying grid @vtype int [num_dims] @vio in @endvar @var coord_arrays @vdesc coordinates of the underlying grid @vtype void * [num_dims] @vio in @endvar @var coord_array_types @vdesc CCTK data types of passed coordinate arrays @vtype int [num_dims] @vio in @endvar @var interp_coord_arrays @vdesc coordinates of points to interpolate at @vtype void * [num_dims] @vio in @endvar @var interp_coord_array_types @vdesc CCTK data types of passed interpolation coordinate arrays @vtype int [num_dims] @vio in @endvar @var in_arrays @vdesc list of input arrays to interpolate from @vtype void * [num_in_arrays] @vio in @endvar @var in_array_types @vdesc CCTK data types of input arrays to interpolate from @vtype int [num_in_arrays] @vio in @endvar @var out_arrays @vdesc list of output arrays to hold the interpolation results @vtype void * [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 return value of the interpolation operator, -1 if invalid interpolation operator handle was passed in @endreturndesc @@*/ int CCTK_InterpLocal (cGH *GH, int operator_handle, int num_points, int num_dims, int num_in_arrays, int num_out_arrays, ...) { va_list indices; int retcode; int *coord_dims; int *coord_array_types, *interp_coord_array_types; int *in_array_types, *out_array_types; const void **coord_arrays, **interp_coord_arrays, **in_arrays; void **out_arrays; struct interp_info *operator; /* Get the interpolation operator */ operator = (struct interp_info *) Util_GetHandledData (interp_operators, operator_handle); if ((operator == NULL) || (operator->interp_operator_local == NULL)) { CCTK_VWarn(0, __LINE__, __FILE__, "Cactus", "\n" " CCTK_InterpLocal(): no interpolation operator is registered\n" " under the handle %d\n" " (did you activate LocalInterp or some other thorn\n" " providing this interpolation operator?)" , operator_handle); retcode = -1; } else { coord_dims = (int *) malloc (num_dims * sizeof (int)); coord_arrays = (const void **) malloc (num_dims * sizeof (void *)); coord_array_types = (int *) malloc (num_dims * sizeof (int)); interp_coord_arrays = (const void **) malloc (num_dims * sizeof (void *)); interp_coord_array_types = (int *) malloc (num_dims * sizeof (int)); in_arrays = (const void **) malloc (num_in_arrays * sizeof (void *)); in_array_types = (int *) malloc (num_in_arrays * sizeof (int)); out_arrays = (void **) malloc (num_out_arrays * sizeof (void *)); out_array_types = (int *) malloc (num_out_arrays * sizeof (int)); /* Fill in the arrays from the variable argument list */ va_start (indices, num_out_arrays); VARARGS_TO_ARRAY (coord_dims, int, NOTHING, num_dims, indices); VARARGS_TO_ARRAY (coord_arrays, void *, NOTHING, num_dims, indices); VARARGS_TO_ARRAY (coord_array_types, int, NOTHING, num_dims, indices); VARARGS_TO_ARRAY (interp_coord_arrays, void *, NOTHING, num_dims, indices); VARARGS_TO_ARRAY (interp_coord_array_types, int, NOTHING, num_dims,indices); VARARGS_TO_ARRAY (in_arrays, void *, NOTHING, num_in_arrays, indices); VARARGS_TO_ARRAY (in_array_types, int, NOTHING, num_in_arrays, indices); VARARGS_TO_ARRAY (out_arrays, void *, NOTHING, num_out_arrays, indices); VARARGS_TO_ARRAY (out_array_types, int, NOTHING, num_out_arrays, indices); va_end (indices); retcode = operator->interp_operator_local (GH, 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); free (out_array_types); free (out_arrays); free (in_array_types); free (in_arrays); free (interp_coord_array_types); free (interp_coord_arrays); free (coord_array_types); free (coord_arrays); free (coord_dims); } return (retcode); } /******************************************************************************/ void CCTK_FCALL CCTK_FNAME (CCTK_InterpLocal) (int *fortranreturn, cGH *GH, const int *operator_handle, const int *num_points, const int *num_dims, const int *num_in_arrays, const int *num_out_arrays, ...) { va_list indices; int retcode; int *coord_dims; int *coord_array_types, *interp_coord_array_types; int *in_array_types, *out_array_types; const void **coord_arrays, **interp_coord_arrays, **in_arrays; void **out_arrays; struct interp_info *operator; /* Get the interpolation operator */ operator = (struct interp_info *) Util_GetHandledData (interp_operators, *operator_handle); if ((operator == NULL) || (operator->interp_operator_local == NULL)) { CCTK_VWarn(0, __LINE__, __FILE__, "Cactus", "\n" " CCTK_InterpLocal(): no interpolation operator is registered\n" " under the handle %d\n" " (did you activate LocalInterp or some other thorn\n" " providing this interpolation operator?)" , *operator_handle); retcode = -1; } else { coord_dims = (int *) malloc (*num_dims * sizeof (int)); coord_arrays = (const void **) malloc (*num_dims * sizeof (void *)); coord_array_types = (int *) malloc (*num_dims * sizeof (int)); interp_coord_arrays = (const void **) malloc (*num_dims * sizeof (void *)); interp_coord_array_types = (int *) malloc (*num_dims * sizeof (int)); in_arrays = (const void **) malloc (*num_in_arrays * sizeof (void *)); in_array_types = (int *) malloc (*num_in_arrays * sizeof (int)); out_arrays = (void **) malloc (*num_out_arrays * sizeof (void *)); out_array_types = (int *) malloc (*num_out_arrays * sizeof (int)); /* Fill in the arrays from the variable argument list */ va_start (indices, num_out_arrays); VARARGS_TO_ARRAY (coord_dims, int *, *, *num_dims, indices); VARARGS_TO_ARRAY (coord_arrays, void *, NOTHING, *num_dims, indices); VARARGS_TO_ARRAY (coord_array_types, int *, *, *num_dims, indices); VARARGS_TO_ARRAY (interp_coord_arrays, void *, NOTHING, *num_dims, indices); VARARGS_TO_ARRAY (interp_coord_array_types, int *, *, *num_dims, indices); VARARGS_TO_ARRAY (in_arrays, void *, NOTHING, *num_in_arrays, indices); VARARGS_TO_ARRAY (in_array_types, int *, *, *num_in_arrays, indices); VARARGS_TO_ARRAY (out_arrays, void *, NOTHING, *num_out_arrays, indices); VARARGS_TO_ARRAY (out_array_types, int *, *, *num_out_arrays, indices); va_end (indices); retcode = operator->interp_operator_local (GH, *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); free (out_array_types); free (out_arrays); free (in_array_types); free (in_arrays); free (interp_coord_array_types); free (interp_coord_arrays); free (coord_array_types); free (coord_arrays); free (coord_dims); } *fortranreturn = retcode; } /*@@ @routine CCTK_InterpGridArrays @date Mon 16 Dec 2002 @author Thomas Radke @desc The general CCTK interpolation routine for grid variables Here only the fortran wrapper is defined which calls the the C routine CCTK_InterpGridArrays(). This is an overloadable routine defined in src/comm/OverloadComm.c. @enddesc @var GH @vdesc pointer to CCTK grid hierarchy @vtype const cGH * @vio in @endvar @var N_dims @vdesc (reference to) number of dimensions for the interpolation @vtype const int * @vio in @endvar @endvar @var param_table_handle @vdesc (reference to) the parameter table handle for passing optional parameters to the interpolator routine @vtype const int * @vio in @endvar @var coord_system_handle @vdesc (reference to) the handle for the underlying coordinate system @vtype const int * @vio in @endvar @var N_interp_points @vdesc (reference to) the number of points to interpolate at @vtype const int * @vio in @endvar @var interp_coords_type @vdesc (reference to) the CCTK datatype of the coordinate arrays as passed via (common datatype for all arrays) @vtype const int * @vio in @endvar @var interp_coords @vdesc list of arrays with coordinate for points to interpolate at @vtype const void *const [] @vio in @endvar @var N_input_arrays @vdesc (reference to) the number of input arrays @vtype const int * @vio in @endvar @var input_array_indices @vdesc list of grid variables (given by their indices) to interpolate @vtype const CCTK_INT [] @vio in @endvar @var N_output_arrays @vdesc (reference to) the number of output arrays @vtype const int * @vio in @endvar @var out_array_types @vdesc list of requested CCTK datatypes for the output arrays @vtype const CCTK_INT [] @vio in @endvar @var output_arrays @vdesc list of output arrays (given by their pointers) which receive the interpolation results @vtype void *const [] @vio out @endvar @returntype int @returndesc return code from routine which overloades CCTK_InterpGridArrays() The return code is passed back in . @endreturndesc @@*/ void CCTK_FCALL CCTK_FNAME (CCTK_InterpGridArrays) (int *ierror, const cGH *GH, const int *N_dims, const int *param_table_handle, const int *local_interp_handle, const int *coord_system_handle, const int *N_interp_points, const int *interp_coords_type, const void *const interp_coords[], const int *N_input_arrays, const CCTK_INT input_array_indices[], const int *N_output_arrays, const CCTK_INT output_array_types[], void *const output_arrays[]) { *ierror = CCTK_InterpGridArrays (GH, *N_dims, *param_table_handle, *local_interp_handle, *coord_system_handle, *N_interp_points, *interp_coords_type, interp_coords, *N_input_arrays,input_array_indices, *N_output_arrays, output_array_types, output_arrays); } /******************************************************************************/ /*@@ @routine CCTK_InterpLocalUniform @date 22 Oct 2001 @author Jonathan Thornburg @desc This API interpolates a set of data arrays defined on a uniform N-dimensional tensor-product grid, to a set of interpolation points. A key-value table is used to pass options to the interpolation operator. @enddesc ***** misc arguments ***** @var N_dims @vdesc dimensionality of the interpolation @vtype int N_dims (must be >= 1) @endvar @var operator_handle @vdesc handle to the interpolation operator @vtype int operator_handle (must be >= 0) @endvar @var param_table_handle @vdesc handle to a key-value table giving additonal parameters for the interpolation @vtype int param_table_handle (must be >= 0) @endvar ***** arguments specifying the coordinate system ***** @var coord_origin @vdesc (pointer to) array[N_dims] of values giving the x,y,z,... coordinates which correspond to the integer input-array subscripts (0,0,0,...,0) (note there is no implication here that such a grid point need actually exist; the arrays just give the coordinates it would have if it did exist) @vtype const CCTK_REAL coord_origin[N_dims] @endvar @var coord_delta @vdesc (pointer to) array[N_dims] of values giving the coordinate spacing of the grid @vtype const CCTK_REAL coord_delta[N_dims] @endvar ***** arguments specifying the interpolation points ***** @var N_interp_points @vdesc number of interpolation points @vtype int N_interp_points (must be >= 0) @endvar @var interp_coords_type_code @vdesc one of the CCTK_VARIABLE_* codes giving the data type of the arrays pointed to by interp_coords[] @vtype int @endvar @var interp_coords @vdesc (pointer to) array[N_dims] of pointers to arrays[N_interp_points] giving x,y,z,... coordinates of interpolation points @vtype const void *const interp_coords[N_dims] @endvar ***** arguments specifying the inputs (the data to be interpolated) ***** @var N_input_arrays @vdesc number of arrays input to the interpolation @vtype int N_input_arrays (must be >= 0) @endvar @var input_array_dims @vdesc dimensions of the input arrays: unless overridden by entries in the parameter table, all input arrays are taken to have these dimensions, with [0] the most contiguous axis and [N_dims-1] the least contiguous axis, and array subscripts in the range 0 <= subscript < dims[axis] @vtype const int input_array_dims[N_dims] @endvar @var input_array_type_codes @vdesc (pointer to) array of CCTK_VARIABLE_* codes giving the data types of the input arrays @vtype const int input_array_type_codes[N_input_arrays] @endvar @var input_arrays @vdesc (pointer to) array[N_input_arrays] of pointers to input arrays @vtype const void *const input_arrays[N_input_arrays] @endvar ***** arguments specifying the outputs (the interpolation results) ***** @var N_output_arrays @vdesc number of arrays output from the interpolation @vtype int N_output_arrays (must be >= 0) @endvar @var output_array_type_codes @vdesc (pointer to) array of CCTK_VARIABLE_* codes giving the data types of the output arrays @vtype const int output_array_type_codes[N_output_arrays] @endvar @var output_arrays @vdesc (pointer to) array[N_output_arrays] of pointers to output arrays @vtype void *const output_arrays[N_output_arrays] @endvar ***** return result ***** @returntype int @returndesc 0 ==> successful, otherwise UTIL_ERROR_BAD_HANDLE ==> bad operator_handle @endreturndesc @@*/ int CCTK_InterpLocalUniform(int N_dims, int operator_handle, int param_table_handle, /***** coordinate system *****/ const CCTK_REAL coord_origin[], const CCTK_REAL coord_delta[], /***** interpolation points *****/ int N_interp_points, int interp_coords_type_code, const void *const interp_coords[], /***** input arrays *****/ int N_input_arrays, const CCTK_INT input_array_dims[], const CCTK_INT input_array_type_codes[], const void *const input_arrays[], /***** output arrays *****/ int N_output_arrays, const CCTK_INT output_array_type_codes[], void *const output_arrays[]) { const struct interp_info *p_interp_info = (struct interp_info *)Util_GetHandledData(interp_operators, operator_handle); if ( (p_interp_info == NULL) || (p_interp_info->interp_op_local_uniform == NULL) ) { CCTK_VWarn(0, __LINE__, __FILE__, "Cactus", "\n" " CCTK_InterpLocalUniform(): no interpolation operator is registered\n" " under the handle %d\n" " (did you activate LocalInterp or some other thorn\n" " providing this interpolation operator?)" , operator_handle); return UTIL_ERROR_BAD_HANDLE; } return p_interp_info->interp_op_local_uniform(N_dims, param_table_handle, /*** coordinate system ***/ coord_origin, coord_delta, /*** interpolation points ***/ N_interp_points, interp_coords_type_code, interp_coords, /***** input arrays *****/ N_input_arrays, input_array_dims, input_array_type_codes, input_arrays, /***** output arrays *****/ N_output_arrays, output_array_type_codes, output_arrays); }