From 70d53d06e35c6d7f3539f73aa80493c84b92477f Mon Sep 17 00:00:00 2001 From: jthorn Date: Mon, 3 Feb 2003 08:58:49 +0000 Subject: finish splitting up (huge) LocalInterp_InterpLocalUniform() function into smaller subfunctions, now it's "only" about 600 lines long (before it was around 1200 lines :( :( :( ) also fix a few more bugs in error reporting and add a few more comments git-svn-id: http://svn.cactuscode.org/arrangements/CactusBase/LocalInterp/trunk@132 df1f8a13-aa1d-4dd4-9681-27ded5b42416 --- .../InterpLocalUniform.c | 2073 ++++++++++++-------- .../InterpLocalUniform.h | 29 +- 2 files changed, 1221 insertions(+), 881 deletions(-) (limited to 'src') diff --git a/src/GeneralizedPolynomial-Uniform/InterpLocalUniform.c b/src/GeneralizedPolynomial-Uniform/InterpLocalUniform.c index 2c1582c..d0b4c78 100644 --- a/src/GeneralizedPolynomial-Uniform/InterpLocalUniform.c +++ b/src/GeneralizedPolynomial-Uniform/InterpLocalUniform.c @@ -1,3 +1,26 @@ +/* InterpLocalUniform.c -- driver for this interpolator */ +/* $Header: */ +/* +** *****data structures and functions local to this file ***** +** + * LocalInterp_ILU_Lagrange - external entry point for Lagrange interpolator + * LocalInterp_ILU_Hermite - external entry point for Hermite interpolator +** +** LocalInterp_InterpLocalUniform - main driver routine +** +** check_boundary_tolerances - check boundary tolerances for consistency +** get_and_decode_molecule_family - get & decode molecule_family from par table +** get_molecule_positions - get molecule_positions from parameter table +** get_Jacobian_query - get Jacobian-query info from parameter table +** set_error_info - set error information in parameter table +** set_molecule_structure - set molecule structure info in parameter table +** set_molecule_min_max_m - set molecule size info in parameter table +** +** get_and_check_INT - get and range-check CCTK_INT from parameter table +** get_INT_array - get CCTK_INT array from parameter table +** get_REAL_array - get CCTK_REAL array from parameter table + */ + /*@@ @file InterpLocalUniform.c @date 22 Oct 2001 @@ -24,7 +47,6 @@ @version $Id$ @@*/ -#include #include #include /* malloc(), free() */ #include @@ -47,6 +69,8 @@ CCTK_FILEVERSION(CactusBase_LocalInterp_src_GeneralizedPolynomialUniform_InterpL * *****data structures and functions local to this file ***** */ +/**************************************/ + /* * data structures local to this file */ @@ -59,6 +83,7 @@ enum interp_operator N_INTERP_OPERATORS /* this must be the last entry in the enum */ }; +/**************************************/ /* * prototypes for functions local to this file @@ -66,25 +91,233 @@ enum interp_operator static int LocalInterp_InterpLocalUniform(enum interp_operator interp_operator, const char* interp_operator_string, + CCTK_REAL boundary_off_cntr_tol_default, + CCTK_REAL boundary_extrap_tol_default, + /***** misc arguments *****/ int N_dims, int param_table_handle, - /***** coordinate system *****/ + /***** coordinate system *****/ const CCTK_REAL coord_origin[], const CCTK_REAL coord_delta[], - /***** interpolation points *****/ + /***** interpolation points *****/ int N_interp_points, int interp_coords_type_code, const void *const interp_coords[], - /***** input arrays *****/ + /***** 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 *****/ + /***** output arrays *****/ int N_output_arrays, const CCTK_INT output_array_type_codes[], void *const output_arrays[]); +static + void check_boundary_tolerances + (int N_boundaries, + const CCTK_REAL boundary_off_centering_tolerance[MAX_N_BOUNDARIES], + const CCTK_REAL boundary_extrapolation_tolerance[MAX_N_BOUNDARIES]); +static + int get_and_decode_molecule_family + (int param_table_handle, + int buffer_size, char molecule_family_string_buffer[], + enum molecule_family *p_molecule_family); +static + int get_molecule_positions + (int param_table_handle, + int N_dims, CCTK_INT* molecule_positions_array[MAX_N_DIMS]); +static + int get_Jacobian_info(int param_table_handle, + int N_dims, int N_output_arrays, + struct Jacobian_info* p_Jacobian_info); +static + int set_error_info(int param_table_handle, + struct error_flags* p_error_flags); +static + int set_molecule_structure + (int param_table_handle, + const struct molecule_structure_flags* p_molecule_structure_flags); +static + int set_molecule_min_max_m + (int param_table_handle, + int N_dims, + const struct molecule_min_max_m_info* p_molecule_min_max_m_info); + +static + int get_and_check_INT + (int param_table_handle, const char name[], + bool mandatory_flag, int default_value, + bool check_range_flag, int min_value, int max_value, + const char max_value_string[], + int* p_value); +static + int get_INT_array(int param_table_handle, const char name[], + bool default_flag, int default_value, + int N, CCTK_INT buffer[], + bool* p_value_not_set); +static + int get_REAL_array(int param_table_handle, const char name[], + CCTK_REAL default_value, + int N, CCTK_REAL buffer[]); + +/**************************************/ + +/* + * table of pointers to actual interpolation functions + */ + +/* + * typedef p_interp_fn_t as a function pointer pointing to an + * individual interpolator function of the sort defined by template.[ch] + */ +#undef FUNCTION_NAME +#define FUNCTION_NAME (*p_interp_fn_t) +typedef + #include "template.h" + +/* NULL (function) pointer of this type */ +#define NULL_INTERP_FN_PTR ((p_interp_fn_t) NULL) + +/* + * For each axis of this array which is marked with a "see above" comment + * in the array declaration, the array size is declared as X+1, so the + * legal subscripts are of course 0...X. But for these axes we actually + * only use 1...X. (Having the array oversize like this slightly simplifies + * the array subscripting.) + */ +static const p_interp_fn_t p_interp_fn_table[N_INTERP_OPERATORS] + [MAX_N_DIMS+1] /* see above */ + [N_MOLECULE_FAMILIES] + [MAX_ORDER+1] /* see above */ + [MAX_SMOOTHING+1] + = { + + /* in the comments on the initializers, "n.i." = "not implemented" */ + { + /* interp_operator == interp_operator_Lagrange */ + { + /* N_dims = 0 ==> unused */ + { + /* molecule_family = molecule_family_cube */ + { NULL_INTERP_FN_PTR, }, /* order=0 (unused) */ + { NULL_INTERP_FN_PTR, }, /* order=1 (unused) */ + { NULL_INTERP_FN_PTR, }, /* order=2 (unused) */ + { NULL_INTERP_FN_PTR, }, /* order=3 (unused) */ + { NULL_INTERP_FN_PTR, }, /* order=4 (unused) */ + { NULL_INTERP_FN_PTR, }, /* order=5 (unused) */ + { NULL_INTERP_FN_PTR, }, /* order=6 (unused) */ + }, + }, + + { + /* N_dims = 1 */ + { + /* molecule_family = molecule_family_cube */ + { NULL_INTERP_FN_PTR, }, /* order=0 (unused) */ + { LocalInterp_ULagrange_1d_cube10, }, /* order=1 */ + { LocalInterp_ULagrange_1d_cube20, }, /* order=2 */ + { LocalInterp_ULagrange_1d_cube30, }, /* order=3 */ + { LocalInterp_ULagrange_1d_cube40, }, /* order=4 */ + { LocalInterp_ULagrange_1d_cube50, }, /* order=5 */ + { LocalInterp_ULagrange_1d_cube60, }, /* order=6 */ + }, + }, + + { + /* N_dims = 2 */ + { + /* molecule_family = molecule_family_cube */ + { NULL_INTERP_FN_PTR, }, /* order=0 (unused) */ + { LocalInterp_ULagrange_2d_cube10, }, /* order=1 */ + { LocalInterp_ULagrange_2d_cube20, }, /* order=2 */ + { LocalInterp_ULagrange_2d_cube30, }, /* order=3 */ + { LocalInterp_ULagrange_2d_cube40, }, /* order=4 */ + { NULL_INTERP_FN_PTR, }, /* order=5 (n.i.) */ + { NULL_INTERP_FN_PTR, }, /* order=6 (n.i.) */ + }, + }, + + { + /* N_dims = 3 */ + { + /* molecule_family = molecule_family_cube */ + { NULL_INTERP_FN_PTR, }, /* order=0 (unused) */ + { LocalInterp_ULagrange_3d_cube10, }, /* order=1 */ + { LocalInterp_ULagrange_3d_cube20, }, /* order=2 */ + { LocalInterp_ULagrange_3d_cube30, }, /* order=3 */ + { LocalInterp_ULagrange_3d_cube40, }, /* order=4 */ + { NULL_INTERP_FN_PTR, }, /* order=5 (n.i.) */ + { NULL_INTERP_FN_PTR, }, /* order=6 (n.i.) */ + }, + }, + + /* end of interp_operator == interp_operator_Lagrange */ + }, + + { + /* interp_operator == interp_operator_Hermite */ + { + /* N_dims = 0 ==> unused */ + { + /* molecule_family = molecule_family_cube */ + { NULL_INTERP_FN_PTR, }, /* order=0 (unused) */ + { NULL_INTERP_FN_PTR, }, /* order=1 (unused) */ + { NULL_INTERP_FN_PTR, }, /* order=2 (unused) */ + { NULL_INTERP_FN_PTR, }, /* order=3 (unused) */ + { NULL_INTERP_FN_PTR, }, /* order=4 (unused) */ + { NULL_INTERP_FN_PTR, }, /* order=5 (unused) */ + { NULL_INTERP_FN_PTR, }, /* order=6 (unused) */ + }, + }, + + { + /* N_dims = 1 */ + { + /* molecule_family = molecule_family_cube */ + { NULL_INTERP_FN_PTR, }, /* order=0 (unused) */ + { NULL_INTERP_FN_PTR, }, /* order=1 (unused) */ + { LocalInterp_UHermite_1d_cube2, }, /* order=2 */ + { LocalInterp_UHermite_1d_cube3, }, /* order=3 */ + { LocalInterp_UHermite_1d_cube4, }, /* order=4 */ + { NULL_INTERP_FN_PTR, }, /* order=5 (n.i.) */ + { NULL_INTERP_FN_PTR, }, /* order=6 (n.i.) */ + }, + }, + + { + /* N_dims = 2 */ + { + /* molecule_family = molecule_family_cube */ + { NULL_INTERP_FN_PTR, }, /* order=0 (unused) */ + { NULL_INTERP_FN_PTR, }, /* order=1 (unused) */ + { LocalInterp_UHermite_2d_cube2, }, /* order=2 */ + { LocalInterp_UHermite_2d_cube3, }, /* order=3 */ + { NULL_INTERP_FN_PTR, }, /* order=4 (n.i.) */ + { NULL_INTERP_FN_PTR, }, /* order=5 (n.i.) */ + { NULL_INTERP_FN_PTR, }, /* order=6 (n.i.) */ + }, + }, + + { + /* N_dims = 3 */ + { + /* molecule_family = molecule_family_cube */ + { NULL_INTERP_FN_PTR, }, /* order=0 (unused) */ + { NULL_INTERP_FN_PTR, }, /* order=1 (unused) */ + { LocalInterp_UHermite_3d_cube2, }, /* order=2 */ + { LocalInterp_UHermite_3d_cube3, }, /* order=3 */ + { NULL_INTERP_FN_PTR, }, /* order=4 (n.i.) */ + { NULL_INTERP_FN_PTR, }, /* order=5 (n.i.) */ + { NULL_INTERP_FN_PTR, }, /* order=6 (n.i.) */ + }, + }, + + /* end of interp_operator == interp_operator_Hermite */ + } + + }; + /******************************************************************************/ /******************************************************************************/ /******************************************************************************/ @@ -131,6 +364,8 @@ int LocalInterp_ILU_Lagrange(int N_dims, { return LocalInterp_InterpLocalUniform(interp_operator_Lagrange, "Lagrange polynomial interpolation", + LAGRANGE_BNDRY_OFF_CNTR_TOL_DEF, + LAGRANGE_BNDRY_EXTRAP_TOL_DEF, N_dims, param_table_handle, coord_origin, coord_delta, @@ -190,6 +425,8 @@ int LocalInterp_ILU_Hermite(int N_dims, { return LocalInterp_InterpLocalUniform(interp_operator_Hermite, "Hermite polynomial interpolation", + HERMITE_BNDRY_OFF_CNTR_TOL_DEF, + HERMITE_BNDRY_EXTRAP_TOL_DEF, N_dims, param_table_handle, coord_origin, coord_delta, @@ -231,10 +468,6 @@ return LocalInterp_InterpLocalUniform(interp_operator_Hermite, and calls the appropriate subfunction to do the actual interpolation, then finally stores some results back in the parameter table. - - FIXME: - This function is huge, and should be split up into - reasonable-sized subfunctions. @enddesc @hdate 28.Jan.2003 @@ -245,7 +478,13 @@ return LocalInterp_InterpLocalUniform(interp_operator_Hermite, @var boundary_extrapolation_tolerance. @endhdesc - ***** misc arguments ***** + @hdate 1.Feb.2003 + @hauthor Jonathan Thornburg + @hdesc Change defaults to "no off-centering" for Hermite interpolator, + split off most of this (huge) function into subfunctions + @endhdesc + + ***** operator arguments ***** @var interp_operator @vdesc describes the interpolation operator @@ -258,6 +497,22 @@ return LocalInterp_InterpLocalUniform(interp_operator_Hermite, @vtype const char* interp_operator_string @endvar + @var boundary_off_cntr_tol_default + @vdesc the default value for each element of + @var boundary_off_centering_default[] + if this key isn't present in the parameter table + @vtype CCTK_REAL boundary_off_cntr_tol_default + @endvar + + @var boundary_extrap_tol_default + @vdesc the default value for each element of + @var boundary_extrapolation_default[] + if this key isn't present in the parameter table + @vtype CCTK_REAL boundary_extrap_tol_default + @endvar + + ***** misc arguments ***** + @var N_dims @vdesc dimensionality of the interpolation @vtype int N_dims (must be >= 1) @@ -358,7 +613,7 @@ return LocalInterp_InterpLocalUniform(interp_operator_Hermite, @vdesc Sets the order of the interpolating polynomial (1=linear, 2=quadratic, 3=cubic, ...). This table entry is mandatory; all others are optional. - @vtype CCTK_INT order + @vtype int order @endvar @var N_boundary_points_to_omit @@ -611,7 +866,7 @@ return LocalInterp_InterpLocalUniform(interp_operator_Hermite, UTIL_ERROR_BAD_INPUT one of the input arguments is invalid UTIL_ERROR_NO_MEMORY unable to malloc() temporary memory UTIL_ERROR_BAD_HANDLE invalid parameter table handle - CCTK_ERROR_INTERP_POINT_X_RANGE + CCTK_ERROR_INTERP_POINT_OUTSIDE interpolation point is out of range (in this case additional information is reported through the parameter table) @@ -622,21 +877,24 @@ return LocalInterp_InterpLocalUniform(interp_operator_Hermite, static int LocalInterp_InterpLocalUniform(enum interp_operator interp_operator, const char* interp_operator_string, + CCTK_REAL boundary_off_cntr_tol_default, + CCTK_REAL boundary_extrap_tol_default, + /***** misc arguments *****/ int N_dims, int param_table_handle, - /***** coordinate system *****/ + /***** coordinate system *****/ const CCTK_REAL coord_origin[], const CCTK_REAL coord_delta[], - /***** interpolation points *****/ + /***** interpolation points *****/ int N_interp_points, int interp_coords_type_code, const void *const interp_coords[], - /***** input arrays *****/ + /***** 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 *****/ + /***** output arrays *****/ int N_output_arrays, const CCTK_INT output_array_type_codes[], void *const output_arrays[]) @@ -662,7 +920,8 @@ static int N_input_arrays1 = N_input_arrays + 1; int N_output_arrays1 = N_output_arrays + 1; -int status, status1, status2, status3, status4; +int status, status1, status2; +bool value_not_set; /******************************************************************************/ @@ -707,285 +966,98 @@ if (N_dims > MAX_N_DIMS) /******************************************************************************/ /* - * interpolation order + * get the parameters from the parameter table, filling in defaults as needed */ + { CCTK_INT order; -status = Util_TableGetInt(param_table_handle, &order, "order"); -if (status == 1) - then { } /* value set from table ==> no-op here */ -else { - /* n.b. order is a mandatory parameter (no default)!*/ - CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, - __LINE__, __FILE__, CCTK_THORNSTRING, -"\n" -" CCTK_InterpLocalUniform(): bad or missing table entry for\n" -" \"order\" parameter!\n" -" (this is a mandatory parameter)\n" -" error status=%d", - status); - return status; /*** ERROR RETURN ***/ - } -if ((order < 1) || (order > MAX_ORDER)) - then { - CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, - __LINE__, __FILE__, CCTK_THORNSTRING, -"\n" -" CCTK_InterpLocalUniform(): implementation restriction: order=%d\n" -" but we only support 1 <= order <= MAX_ORDER=%d!", - order, - MAX_ORDER); - return UTIL_ERROR_BAD_INPUT; /*** ERROR RETURN ***/ - } - -/******************************************************************************/ +status = get_and_check_INT(param_table_handle, "order", + true, 0, /* this is a mandatory parameter */ + true, 1, MAX_ORDER, "MAX_ORDER", /* range check */ + &order); +if (status != 0) + then return status; /*** ERROR RETURN ***/ -/* - * out-of-range interpolation-point handling - */ { +const int N_boundaries = 2*N_dims; CCTK_INT N_boundary_points_to_omit[MAX_N_BOUNDARIES]; +status = get_INT_array(param_table_handle, "N_boundary_points_to_omit", + true, 0, /* default value */ + N_boundaries, N_boundary_points_to_omit, + NULL); +if (status != 0) + then return status; /*** ERROR RETURN ***/ + +/**************************************/ + + { CCTK_REAL boundary_off_centering_tolerance[MAX_N_BOUNDARIES]; CCTK_REAL boundary_extrapolation_tolerance[MAX_N_BOUNDARIES]; -const int N_boundaries = 2*N_dims; +status = get_REAL_array(param_table_handle, "boundary_off_centering_tolerance", + boundary_off_cntr_tol_default, + N_boundaries, boundary_off_centering_tolerance); +if (status != 0) + then return status; /*** ERROR RETURN ***/ +status = get_REAL_array(param_table_handle, "boundary_extrapolation_tolerance", + boundary_extrap_tol_default, + N_boundaries, boundary_extrapolation_tolerance); +if (status != 0) + then return status; /*** ERROR RETURN ***/ +check_boundary_tolerances(N_boundaries, + boundary_off_centering_tolerance, + boundary_extrapolation_tolerance); -status = Util_TableGetIntArray(param_table_handle, - N_boundaries, N_boundary_points_to_omit, - "N_boundary_points_to_omit"); -if (status == UTIL_ERROR_TABLE_NO_SUCH_KEY) - then { - /* set the default value */ - int ibndry; - for (ibndry = 0 ; ibndry < N_boundaries ; ++ibndry) - { - N_boundary_points_to_omit[ibndry] = 0; - } - } -else if (status == N_boundaries) - then { } /* value set from table ==> no-op here */ -else { - CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, - __LINE__, __FILE__, CCTK_THORNSTRING, -"\n" -" CCTK_InterpLocalUniform(): bad table entry for\n" -" \"N_boundary_points_to_omit\" parameter!\n" -" error status=%d", - status); - return status; /*** ERROR RETURN ***/ - } +/**************************************/ + + { +#define MOLECULE_FAMILY_BUFSIZ (MAX_MOLECULE_FAMILY_STRLEN+1) +char molecule_family_string[MOLECULE_FAMILY_BUFSIZ]; +enum molecule_family molecule_family; +status = get_and_decode_molecule_family(param_table_handle, + MOLECULE_FAMILY_BUFSIZ, molecule_family_string, + &molecule_family); +if (status != 0) + then return status; /*** ERROR RETURN ***/ + + { +CCTK_INT smoothing; +status = get_and_check_INT(param_table_handle, "smoothing", + false, 0, /* optional, default value 0 */ + true, 0, MAX_SMOOTHING, "MAX_SMOOTHING", + /* range check */ + &smoothing); +if (status != 0) + then return status; /*** ERROR RETURN ***/ /**************************************/ -status = Util_TableGetRealArray(param_table_handle, - N_boundaries, boundary_off_centering_tolerance, - "boundary_off_centering_tolerance"); -if (status == UTIL_ERROR_TABLE_NO_SUCH_KEY) + { +CCTK_INT* const input_array_offsets + = (CCTK_INT*) malloc(N_input_arrays1 * sizeof(CCTK_INT)); +if (input_array_offsets == NULL) + then return UTIL_ERROR_NO_MEMORY; /*** ERROR RETURN ***/ +status = get_INT_array(param_table_handle, "input_array_offsets", + true, 0, /* default value */ + N_input_arrays, input_array_offsets, + NULL); +if (status != 0) then { - /* set the default value */ - int ibndry; - for (ibndry = 0 ; ibndry < N_boundaries ; ++ibndry) - { - boundary_off_centering_tolerance[ibndry] - = BOUNDARY_OFF_CENTER_TOL_DEFAULT; - } - } -else if (status == N_boundaries) - then { } /* value set from table ==> no-op here */ -else { - CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, - __LINE__, __FILE__, CCTK_THORNSTRING, -"\n" -" CCTK_InterpLocalUniform(): bad table entry for\n" -" \"boundary_off_centering_tolerance\" parameter!\n" -" error status=%d", - status); + free(input_array_offsets); return status; /*** ERROR RETURN ***/ } -/**************************************/ - -status = Util_TableGetRealArray(param_table_handle, - N_boundaries, boundary_extrapolation_tolerance, - "boundary_extrapolation_tolerance"); -if (status == UTIL_ERROR_TABLE_NO_SUCH_KEY) - then { - /* set the default value */ - int ibndry; - for (ibndry = 0 ; ibndry < N_boundaries ; ++ibndry) - { - boundary_extrapolation_tolerance[ibndry] - = BOUNDARY_EXTRAP_TOL_DEFAULT; - } - } -else if (status == N_boundaries) - then { } /* value set from table ==> no-op here */ -else { - CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, - __LINE__, __FILE__, CCTK_THORNSTRING, -"\n" -" CCTK_InterpLocalUniform(): bad table entry for\n" -" \"boundary_extrapolation_tolerance\" parameter!\n" -" error status=%d", - status); - return status; /*** ERROR RETURN ***/ - } - -/* check for almost-always-a-mistake settings, warn if found */ - { -int ibndry; - for (ibndry = 0 ; ibndry < N_boundaries ; ++ibndry) - { - if ( (boundary_off_centering_tolerance[ibndry] == 0.0) - && (boundary_extrapolation_tolerance[ibndry] > 0.0) ) - then CCTK_VWarn(WARNING_MSG_SEVERITY_LEVEL, - __LINE__, __FILE__, CCTK_THORNSTRING, -"\n" -" CCTK_InterpLocalUniform():\n" -" warning: setting boundary_off_centering_tolerance[%d] == 0.0\n" -" and boundary_extrapolation_tolerance[%d] == %g > 0.0\n" -" is almost certainly a mistake\n" -" (the boundary_extrapolation_tolerance[%d] == %g\n" -" setting will have no effect because of the\n" -" boundary_off_centering_tolerance[%d] == 0.0 setting)" - , - ibndry, - ibndry, - (double) boundary_extrapolation_tolerance[ibndry], - ibndry, - (double) boundary_extrapolation_tolerance[ibndry], - ibndry); - } - } - -/******************************************************************************/ - -/* - * molecule family - */ - { -#define MOLECULE_FAMILY_BUFSIZ (MAX_MOLECULE_FAMILY_STRLEN+1) -char molecule_family_strbuf[MOLECULE_FAMILY_BUFSIZ]; -char* molecule_family_str = molecule_family_strbuf; -enum molecule_family molecule_family; -status = Util_TableGetString(param_table_handle, - MOLECULE_FAMILY_BUFSIZ, molecule_family_strbuf, - "molecule_family"); -if (status == UTIL_ERROR_TABLE_NO_SUCH_KEY) - then { - /* set the default value for code here*/ - /* (need enum for main code, string for error messages) */ - molecule_family = molecule_family_cube; - molecule_family_str = "cube"; - - /* set this key in the parameter table */ - /* to give the (default) molecule family we're going to use */ - status = Util_TableSetString(param_table_handle, - molecule_family_str, - "molecule_family"); - if (status < 0) - then { - CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, - __LINE__, __FILE__, CCTK_THORNSTRING, -"\n" -" CCTK_InterpLocalUniform(): error setting default molecule family\n" -" in parameter table!\n" -" error status=%d", - status); - return status; /*** ERROR RETURN ***/ - } - } -else if (status > 0) - then { - /* decode molecule family string */ - if (strcmp(molecule_family_strbuf, "cube") == 0) - then molecule_family = molecule_family_cube; - else { - CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, - __LINE__, __FILE__, CCTK_THORNSTRING, -"CCTK_InterpLocalUniform(): unknown molecule_family=\"%s\"!", - molecule_family_strbuf); - return UTIL_ERROR_BAD_INPUT; /*** ERROR RETURN ***/ - } - } -else { - CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, - __LINE__, __FILE__, CCTK_THORNSTRING, -"\n" -" CCTK_InterpLocalUniform(): bad table entry for\n" -" \"molecule_family\" parameter!\n" -" error status=%d", - status); - return status; /*** ERROR RETURN ***/ - } - -/**************************************/ - -/* - * smoothing - */ - { -CCTK_INT smoothing; -status = Util_TableGetInt(param_table_handle, &smoothing, "smoothing"); -if (status == UTIL_ERROR_TABLE_NO_SUCH_KEY) - then smoothing = 0; /* set default value */ -else if (status == 1) - then { } /* value set from table ==> no-op here */ -else { - CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, - __LINE__, __FILE__, CCTK_THORNSTRING, -"\n" -" CCTK_InterpLocalUniform(): bad table entry for\n" -" \"smoothing\" parameter!\n" -" error status=%d", - status); - return status; /*** ERROR RETURN ***/ - } -if ((smoothing < 0) || (smoothing > MAX_SMOOTHING)) - then return UTIL_ERROR_BAD_INPUT; /*** ERROR RETURN ***/ - -/******************************************************************************/ - -/* - * input array offsets - */ { -CCTK_INT* const input_array_offsets - = (CCTK_INT*) malloc(N_input_arrays1 * sizeof(CCTK_INT)); -if (input_array_offsets == NULL) - then return UTIL_ERROR_NO_MEMORY; /*** ERROR RETURN ***/ -status = Util_TableGetIntArray(param_table_handle, - N_input_arrays, input_array_offsets, - "input_array_offsets"); -if (status == UTIL_ERROR_TABLE_NO_SUCH_KEY) +CCTK_INT input_array_strides[MAX_N_DIMS]; +status = get_INT_array(param_table_handle, "input_array_strides", + false, 0, /* don't set default value */ + N_dims, input_array_strides, + &value_not_set); +if (status != 0) then { - /* default offset = 0 along each axis */ - LocalInterp_zero_int_array(N_input_arrays, input_array_offsets); - } -else if (status == N_input_arrays) - then { } /* ok ==> no-op here */ -else { free(input_array_offsets); - CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, - __LINE__, __FILE__, CCTK_THORNSTRING, -"\n" -" CCTK_InterpLocalUniform(): bad table entry for\n" -" \"input_array_offsets\" parameter!\n" -" error status=%d", - status); return status; /*** ERROR RETURN ***/ } - -/**************************************/ - -/* - * input array strides - */ - { -CCTK_INT input_array_strides[MAX_N_DIMS]; -status = Util_TableGetIntArray(param_table_handle, - N_dims, input_array_strides, - "input_array_strides"); -if (status == UTIL_ERROR_TABLE_NO_SUCH_KEY) +if (value_not_set) then { /* default stride = Fortran storage order, */ /* determined from input_array_dims[] */ @@ -1009,54 +1081,30 @@ if (status == UTIL_ERROR_TABLE_NO_SUCH_KEY) } } } -else if (status == N_dims) - then { } /* ok ==> no-op here */ -else { - free(input_array_offsets); - CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, - __LINE__, __FILE__, CCTK_THORNSTRING, -"\n" -" CCTK_InterpLocalUniform(): bad table entry for\n" -" \"input_array_offsets\" parameter!\n" -" error status=%d", - status); - return status; /*** ERROR RETURN ***/ - } -/**************************************/ - -/* - * input array min/max subscripts - */ { CCTK_INT input_array_min_subscripts[MAX_N_DIMS]; -status = Util_TableGetIntArray(param_table_handle, - N_dims, input_array_min_subscripts, - "input_array_min_subscripts"); -if (status == UTIL_ERROR_TABLE_NO_SUCH_KEY) +CCTK_INT input_array_max_subscripts[MAX_N_DIMS]; +status = get_INT_array(param_table_handle, "input_array_min_subscripts", + true, 0, /* default value */ + N_dims, input_array_min_subscripts, + NULL); +if (status != 0) then { - /* default min subscript = 0 along each axis */ - LocalInterp_zero_int_array(N_dims, input_array_min_subscripts); + free(input_array_offsets); + return status; /*** ERROR RETURN ***/ } -else if (status == N_dims) - then { } /* ok ==> no-op here */ -else { + +status = get_INT_array(param_table_handle, "input_array_max_subscripts", + false, 0, /* don't set default value */ + N_dims, input_array_max_subscripts, + &value_not_set); +if (status != 0) + then { free(input_array_offsets); - CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, - __LINE__, __FILE__, CCTK_THORNSTRING, -"\n" -" CCTK_InterpLocalUniform(): bad table entry for\n" -" \"input_array_min_subscripts\" parameter!\n" -" error status=%d", - status); return status; /*** ERROR RETURN ***/ } - { -CCTK_INT input_array_max_subscripts[MAX_N_DIMS]; -status = Util_TableGetIntArray(param_table_handle, - N_dims, input_array_max_subscripts, - "input_array_max_subscripts"); -if (status == UTIL_ERROR_TABLE_NO_SUCH_KEY) +if (value_not_set) then { /* default max subscript = input_array_dims[]-1 along each axis */ if (input_array_dims == NULL) @@ -1077,25 +1125,9 @@ if (status == UTIL_ERROR_TABLE_NO_SUCH_KEY) } } } -else if (status == N_dims) - then { } /* ok ==> no-op here */ -else { - free(input_array_offsets); - CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, - __LINE__, __FILE__, CCTK_THORNSTRING, -"\n" -" CCTK_InterpLocalUniform(): bad table entry for\n" -" \"input_array_max_subscripts\" parameter!\n" -" error status=%d", - status); - return status; /*** ERROR RETURN ***/ - } -/******************************************************************************/ +/**************************************/ -/* - * operand indices - */ { CCTK_INT* const operand_indices = (CCTK_INT*) malloc(N_output_arrays1 * sizeof(CCTK_INT)); @@ -1104,23 +1136,30 @@ if (operand_indices == NULL) free(input_array_offsets); return UTIL_ERROR_NO_MEMORY; /*** ERROR RETURN ***/ } -status = Util_TableGetIntArray(param_table_handle, - N_output_arrays, operand_indices, - "operand_indices"); -if (status == UTIL_ERROR_TABLE_NO_SUCH_KEY) +status = get_INT_array(param_table_handle, "operand_indices", + false, 0, /* don't set default value */ + N_output_arrays, operand_indices, + &value_not_set); +if (status != 0) + then { + free(input_array_offsets); + free(operand_indices); + return UTIL_ERROR_NO_MEMORY; /*** ERROR RETURN ***/ + } +if (value_not_set) then { /* default operand will use each input exactly once, */ /* but this only makes since if N_input_arrays == N_output_arrays */ if (N_input_arrays != N_input_arrays) then { - free(operand_indices); free(input_array_offsets); + free(operand_indices); CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, __LINE__, __FILE__, CCTK_THORNSTRING, "\n" " CCTK_InterpLocalUniform(): can't default \"operand_indices\"\n" " with N_input_arrays=%d != N_output_arrays=%d!" -, + , N_input_arrays, N_output_arrays); return UTIL_ERROR_BAD_INPUT; /*** ERROR RETURN ***/ } @@ -1132,102 +1171,40 @@ if (status == UTIL_ERROR_TABLE_NO_SUCH_KEY) } } } -else if (status == N_output_arrays) - then { - /* check that all operands are within range */ - int out; - for (out = 0 ; out < N_output_arrays ; ++out) - { - if ( (operand_indices[out] < 0) - || (operand_indices[out] >= N_input_arrays) ) - then { - free(operand_indices); - free(input_array_offsets); - CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, - __LINE__, __FILE__, CCTK_THORNSTRING, -"\n" -" CCTK_InterpLocalUniform(): operand_indices[out=%d]=%d specifies\n" -" nonexistent input array!\n" -" (valid range is 0 <= value < N_input_arrays=%d)" -, - out, operand_indices[out], - N_input_arrays); - return UTIL_ERROR_BAD_INPUT; /*** ERROR RETURN ***/ - } - } - } -else { - free(operand_indices); - free(input_array_offsets); - CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, - __LINE__, __FILE__, CCTK_THORNSTRING, -"\n" -" CCTK_InterpLocalUniform(): bad table entry for\n" -" \"operand_indices\" parameter!\n" -" error status=%d", - status); - return status; /*** ERROR RETURN ***/ - } /**************************************/ -/* - * operation_codes - */ { CCTK_INT* const operation_codes = (CCTK_INT*) malloc(N_output_arrays1 * sizeof(CCTK_INT)); if (operation_codes == NULL) then { - free(operand_indices); free(input_array_offsets); + free(operand_indices); return UTIL_ERROR_NO_MEMORY; /*** ERROR RETURN ***/ } -status = Util_TableGetIntArray(param_table_handle, - N_output_arrays, operation_codes, - "operation_codes"); -if (status == UTIL_ERROR_TABLE_NO_SUCH_KEY) +status = get_INT_array(param_table_handle, "operation_codes", + true, 0, /* default value */ + N_output_arrays, operation_codes, + NULL); +if (status != 0) then { - /* default operation_codes = all 0 */ - /* but this only makes since if N_input_arrays == N_output_arrays */ - if (N_input_arrays != N_input_arrays) - then { - free(operation_codes); - free(operand_indices); - free(input_array_offsets); - CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, - __LINE__, __FILE__, CCTK_THORNSTRING, -"\n" -" CCTK_InterpLocalUniform(): can't default \"operation_codes\"\n" -" with N_input_arrays=%d != N_output_arrays=%d!" -, - N_input_arrays, N_output_arrays); - return UTIL_ERROR_BAD_INPUT; /*** ERROR RETURN ***/ - } - LocalInterp_zero_int_array(N_output_arrays, operation_codes); - } -else if (status == N_output_arrays) - then { } /* ok ==> no-op here */ -else { - free(operation_codes); - free(operand_indices); free(input_array_offsets); - CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, - __LINE__, __FILE__, CCTK_THORNSTRING, -"\n" -" CCTK_InterpLocalUniform(): bad table entry for\n" -" \"operation_codes\" parameter!\n" -" error status=%d", - status); - return status; /*** ERROR RETURN ***/ + free(operand_indices); + free(operation_codes); + return UTIL_ERROR_NO_MEMORY; /*** ERROR RETURN ***/ } /******************************************************************************/ /* - * set up for any molecule min/max m and position queries + * set up for any molecule min/max m and/or position queries */ + { +struct molecule_min_max_m_info* p_molecule_min_max_m_info = NULL; +struct molecule_min_max_m_info molecule_min_max_m_info; + /* molecule min/max m */ status1 = Util_TableQueryValueInfo(param_table_handle, NULL, NULL, @@ -1237,60 +1214,54 @@ status2 = Util_TableQueryValueInfo(param_table_handle, "molecule_max_m"); if ((status1 < 0) || (status2 < 0)) then { - free(operation_codes); - free(operand_indices); free(input_array_offsets); + free(operand_indices); + free(operation_codes); CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, __LINE__, __FILE__, CCTK_THORNSTRING, "\n" -" CCTK_InterpLocalUniform(): bad \"interp_molecule_{min,max}_m\"\n" -" table entry/entries!\n" +" CCTK_InterpLocalUniform(): bad \"molecule_{min,max}_m\"\n" +" table entry/entries in query!\n" " error status1=%d status2=%d" -, + , status1, status2); return (status1 < 0) ? status1 : status2; /*** ERROR RETURN ***/ } - { -struct molecule_min_max_m_info molecule_min_max_m_info; -struct molecule_min_max_m_info *p_molecule_min_max_m_info - = (status1 && status2) ? &molecule_min_max_m_info : NULL; +if (status1 && status2) + then p_molecule_min_max_m_info = &molecule_min_max_m_info; -/* molecule positions */ { -CCTK_INT* const* p_molecule_positions = NULL; -CCTK_POINTER molecule_positions_temp [MAX_N_DIMS]; -CCTK_INT* molecule_positions_array[MAX_N_DIMS]; -status = Util_TableGetPointerArray(param_table_handle, - N_dims, molecule_positions_temp, - "molecule_positions"); -if (status == UTIL_ERROR_TABLE_NO_SUCH_KEY) - then { - /* no querying of interpolator molecule positions */ - } -else if (status == N_dims) +CCTK_INT** p_molecule_positions = NULL; +CCTK_INT* molecule_positions_array[MAX_N_DIMS]; + +/* are we doing a molecule-positions query? */ +status = Util_TableQueryValueInfo(param_table_handle, + NULL, NULL, + "molecule_positions"); +if (status < 0) then { - /* type-cast CCTK_POINTER pointers to CCTK_INT* pointers */ - /* (which point to where the query results will be stored) */ - int axis; - for (axis = 0 ; axis < N_dims ; ++axis) - { - molecule_positions_array[axis] - = (CCTK_INT *) molecule_positions_temp[axis]; - } - p_molecule_positions = molecule_positions_array; - } -else { - free(operation_codes); - free(operand_indices); free(input_array_offsets); + free(operand_indices); + free(operation_codes); CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, __LINE__, __FILE__, CCTK_THORNSTRING, "\n" -" CCTK_InterpLocalUniform(): bad \"molecule_positions\" table entry!\n" -" error status=%d", +" CCTK_InterpLocalUniform(): bad \"molecule_positions\"\n" +" table entry in query!\n" +" error status=%d" + , status); return status; /*** ERROR RETURN ***/ } +if (status) + then { + /* yes, we're doing a molecule-positions query */ + status = get_molecule_positions(param_table_handle, + N_dims, molecule_positions_array); + if (status != 0) + then return status; /*** ERROR RETURN ***/ + p_molecule_positions = molecule_positions_array; + } /******************************************************************************/ @@ -1299,470 +1270,611 @@ else { */ { -struct Jacobian_info *p_Jacobian_info = NULL; -struct Jacobian_info Jacobian_info; +struct Jacobian_info* p_Jacobian_info = NULL; +struct Jacobian_info Jacobian_info; Jacobian_info.Jacobian_pointer = NULL; Jacobian_info.Jacobian_offset = NULL; -/* first find out if we are doing Jacobian queries at all */ +/* are we doing a Jacobian query? */ status = Util_TableQueryValueInfo(param_table_handle, NULL, NULL, "Jacobian_pointer"); if (status < 0) then { - free(operation_codes); - free(operand_indices); free(input_array_offsets); + free(operand_indices); + free(operation_codes); CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, __LINE__, __FILE__, CCTK_THORNSTRING, "\n" " CCTK_InterpLocalUniform(): bad \"Jacobian_pointer\" table entry!\n" " (preliminary check)\n" " error status=%d" -, + , status); return status; /*** ERROR RETURN ***/ } if (status) then { - /* yes, we're doing Jacobian queries */ - Jacobian_info.Jacobian_pointer - = (CCTK_REAL**) malloc(N_output_arrays1 * sizeof(CCTK_REAL *)); - if (Jacobian_info.Jacobian_pointer == NULL) + /* yes, we're doing a Jacobian query */ + status = get_Jacobian_info(param_table_handle, + N_dims, N_output_arrays, + &Jacobian_info); + if (status != 0) then { - free(operation_codes); - free(operand_indices); free(input_array_offsets); - return UTIL_ERROR_NO_MEMORY; /*** ERROR RETURN ***/ - } - { - CCTK_POINTER* Jacobian_pointer_temp - = (CCTK_POINTER*) malloc(N_output_arrays1 * sizeof(CCTK_POINTER)); - if (Jacobian_pointer_temp == NULL) - then { - free(operation_codes); free(operand_indices); - free(input_array_offsets); - free(Jacobian_info.Jacobian_pointer); - return UTIL_ERROR_NO_MEMORY; /*** ERROR RETURN ***/ - } - status = Util_TableGetPointerArray(param_table_handle, - N_output_arrays, - Jacobian_pointer_temp, - "Jacobian_pointer"); - if (status == N_output_arrays) - then { - /* type-cast CCTK_POINTER pointers to CCTK_REAL* pointers */ - /* (which point to where the query results will be stored) */ - int out; - for (out = 0 ; out < N_output_arrays ; ++out) - { - Jacobian_info.Jacobian_pointer[out] - = (CCTK_REAL *) Jacobian_pointer_temp[out]; - } - free(Jacobian_pointer_temp); - } - else { free(operation_codes); - free(operand_indices); - free(input_array_offsets); - free(Jacobian_info.Jacobian_pointer); - free(Jacobian_pointer_temp); - CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, - __LINE__, __FILE__, CCTK_THORNSTRING, -"\n" -" CCTK_InterpLocalUniform(): bad \"Jacobian_pointer\" table entry!\n" -" error status=%d" -, - status); return status; /*** ERROR RETURN ***/ } - } + p_Jacobian_info = & Jacobian_info; + } + +/******************************************************************************/ + +/* + * decode (interp_operator, N_dims, molecule_family, order, smoothing) + * and call the appropriate subfunction to do the actual interpolation + */ + + +/* firewall: make sure all the lookup indices are in range */ +/* ... commented-out comparisons are always true since enums are unsigned */ +if ( /*(interp_operator >= 0) &&*/ (interp_operator < N_INTERP_OPERATORS) + && (N_dims >= 0) && (N_dims <= MAX_N_DIMS) + && /*(molecule_family >= 0) &&*/ (molecule_family < N_MOLECULE_FAMILIES) + && (order >= 0) && (order <= MAX_ORDER) + && (smoothing >= 0) && (smoothing <= MAX_SMOOTHING) ) + then { + /* ok ==> no-op here */ + } + else { + CCTK_VWarn(BUG_MSG_SEVERITY_LEVEL, + __LINE__, __FILE__, CCTK_THORNSTRING, +"\n" +" CCTK_InterpLocalUniform():\n" +" internal error (interpolator bug):\n" +" one or more function-pointer-table indices is out of range!\n" +" interp_operator=(int)%d, should be in [0,%d)\n" +" N_dims = %d, should be in [0,%d]\n" +" molecule_family=(int)%d, should be in [0,%d)\n" +" order = %d, should be in [0,%d]\n" +" smoothing = %d, should be in [0,%d]\n" + , + (int) interp_operator, (int) N_INTERP_OPERATORS, + (int) N_dims , (int) MAX_N_DIMS, + (int) molecule_family, (int) N_MOLECULE_FAMILIES, + (int) order , (int) MAX_ORDER, + (int) smoothing , (int) MAX_SMOOTHING); + /*NOTREACHED*/ + CCTK_Abort(NULL, BUG_ABORT_CODE); /*NOTREACHED*/ + } - Jacobian_info.Jacobian_offset - = (CCTK_INT*) malloc(N_output_arrays1 * sizeof(CCTK_INT)); - if (Jacobian_info.Jacobian_offset == NULL) +/* look up the subfunction to do the interpolation */ + { +const p_interp_fn_t p_interp_fn = p_interp_fn_table[interp_operator] + [N_dims] + [molecule_family] + [order] + [smoothing]; +if (p_interp_fn == NULL_INTERP_FN_PTR) + then { + free(input_array_offsets); + free(operand_indices); + free(operation_codes); + if (p_Jacobian_info != NULL) then { - free(operation_codes); - free(operand_indices); - free(input_array_offsets); + free(Jacobian_info.Jacobian_offset); free(Jacobian_info.Jacobian_pointer); - return UTIL_ERROR_NO_MEMORY; /*** ERROR RETURN ***/ - } - status = Util_TableGetIntArray(param_table_handle, - N_output_arrays, - Jacobian_info.Jacobian_offset, - "Jacobian_offset"); - if (status == N_output_arrays) - then { - /* ok ==> no-op here */ } - else if (status == UTIL_ERROR_TABLE_NO_SUCH_KEY) + CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, + __LINE__, __FILE__, CCTK_THORNSTRING, +"\n" +" CCTK_InterpLocalUniform():\n" +" interpolation not implemented for the combination\n" +" interp_operator=\"%s\", N_dims=%d\n" +" molecule_family=\"%s\", order=%d, smoothing=%d", + interp_operator_string, N_dims, + molecule_family_string, order, smoothing); + return UTIL_ERROR_BAD_INPUT; /*** ERROR RETURN ***/ + } + +/* call the subfunction to actually do the interpolation */ + { +struct error_flags error_flags; +struct molecule_structure_flags molecule_structure_flags; +const int return_code + = (*p_interp_fn)(coord_origin, coord_delta, + N_interp_points, + interp_coords_type_code, interp_coords, + N_boundary_points_to_omit, + boundary_off_centering_tolerance, + boundary_extrapolation_tolerance, + N_input_arrays, + input_array_offsets, input_array_strides, + input_array_min_subscripts, + input_array_max_subscripts, + input_array_type_codes, input_arrays, + N_output_arrays, + output_array_type_codes, output_arrays, + operand_indices, operation_codes, + &error_flags, + &molecule_structure_flags, + p_molecule_min_max_m_info, + p_molecule_positions, + p_Jacobian_info); + +/******************************************************************************/ + +/* + * is an interpolation point out of range? + */ +if (return_code == CCTK_ERROR_INTERP_POINT_OUTSIDE) + then { + status = set_error_info(param_table_handle, + &error_flags); + if (status != 0) then { - /* default offset = all 0 */ - LocalInterp_zero_int_array(N_output_arrays, - Jacobian_info.Jacobian_offset); - } - else { - free(operation_codes); - free(operand_indices); free(input_array_offsets); - free(Jacobian_info.Jacobian_pointer); - free(Jacobian_info.Jacobian_offset); - CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, - __LINE__, __FILE__, CCTK_THORNSTRING, -"\n" -" CCTK_InterpLocalUniform(): bad \"Jacobian_offset\" table entry!\n" -" error status=%d" -, - status); + free(operand_indices); + free(operation_codes); + if (p_Jacobian_info != NULL) + then { + free(Jacobian_info.Jacobian_offset); + free(Jacobian_info.Jacobian_pointer); + } return status; /*** ERROR RETURN ***/ } + } + +/******************************************************************************/ + +/* + * store query results + * + * ... only molecule structure and min/max m have to be stored explicitly; + * the other queries are stored directly into their final destinatios + * by the interpolation subfunction + */ - status1 = Util_TableGetInt(param_table_handle, - & Jacobian_info.Jacobian_interp_point_stride, - "Jacobian_interp_point_stride"); - status2 = Util_TableGetIntArray(param_table_handle, - N_dims, - Jacobian_info.Jacobian_m_strides, - "Jacobian_m_strides"); - if ((status1 < 0) || (status2 != N_dims)) +status = set_molecule_structure(param_table_handle, + &molecule_structure_flags); +if (status != 0) + then { + free(input_array_offsets); + free(operand_indices); + free(operation_codes); + if (p_Jacobian_info != NULL) then { - free(operation_codes); - free(operand_indices); - free(input_array_offsets); - free(Jacobian_info.Jacobian_pointer); free(Jacobian_info.Jacobian_offset); - CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, - __LINE__, __FILE__, CCTK_THORNSTRING, -"\n" -" CCTK_InterpLocalUniform(): bad \"Jacobian_interp_point_stride\" and/or\n" -" \"Jacobian_m_strides\" table entry/entries!\n" -" error status1=%d status2=%d" -, - status1, status2); - return (status1 < 0) ? status1 : status2; /*** ERROR RETURN ***/ + free(Jacobian_info.Jacobian_pointer); } + return status; /*** ERROR RETURN ***/ + } - status = Util_TableGetInt(param_table_handle, - & Jacobian_info.Jacobian_part_stride, - "Jacobian_part_stride"); - if (status == 1) +if (p_molecule_min_max_m_info != NULL) + then { + status = set_molecule_min_max_m(param_table_handle, + N_dims, p_molecule_min_max_m_info); + if (status != 0) then { - /* ok ==> no-op here */ - } - else if (status == UTIL_ERROR_TABLE_NO_SUCH_KEY) - then Jacobian_info.Jacobian_part_stride = 0; /* default */ - else { - free(operation_codes); - free(operand_indices); free(input_array_offsets); - free(Jacobian_info.Jacobian_pointer); - free(Jacobian_info.Jacobian_offset); - CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, - __LINE__, __FILE__, CCTK_THORNSTRING, -"\n" -" CCTK_InterpLocalUniform(): bad \"Jacobian_part_stride\" table entry!\n" -" error status=%d" -, - status); + free(operand_indices); + free(operation_codes); + if (p_Jacobian_info != NULL) + then { + free(Jacobian_info.Jacobian_offset); + free(Jacobian_info.Jacobian_pointer); + } return status; /*** ERROR RETURN ***/ } - - p_Jacobian_info = & Jacobian_info; } /******************************************************************************/ -/* - * decode (interp_operator, N_dims, molecule_family, order, smoothing) - * and call the appropriate subfunction to do the actual interpolation - */ - { -/* - * typedef p_interp_fn_t as a function pointer pointing to an - * individual interpolator function of the sort defined by template.[ch] - */ -#undef FUNCTION_NAME -#define FUNCTION_NAME (*p_interp_fn_t) -typedef - #include "template.h" +free(input_array_offsets); +free(operand_indices); +free(operation_codes); +if (p_Jacobian_info != NULL) + then { + free(Jacobian_info.Jacobian_offset); + free(Jacobian_info.Jacobian_pointer); + } -/* NULL (function) pointer of this type */ -#define NULL_INTERP_FN_PTR ((p_interp_fn_t) NULL) +return return_code; /*** NORMAL RETURN ***/ + } + } + } + } + } + } + } + } + } + } + } + } + } + } + } +} + +/******************************************************************************/ +/******************************************************************************/ +/******************************************************************************/ /* - * For each axis of this array which is marked with a "see above" comment - * in the array declaration, the array size is declared as X+1, so the - * legal subscripts are of course 0...X. But for these axes we actually - * only use 1...X. (Having the array oversize like this slightly simplifies - * the array subscripting.) + * This function checks + * CCTK_REAL boundary_off_centering_tolerance[N_boundaries] + * CCTK_REAL boundary_extrapolation_tolerance[N_boundaries] + * for validity, warning about any dubious settings. */ -static const p_interp_fn_t p_interp_fn_table[N_INTERP_OPERATORS] - [MAX_N_DIMS+1] /* see above */ - [N_MOLECULE_FAMILIES] - [MAX_ORDER+1] /* see above */ - [MAX_SMOOTHING+1] - = { - - /* in the comments on the initializers, "n.i." = "not implemented" */ +static + void check_boundary_tolerances + (int N_boundaries, + const CCTK_REAL boundary_off_centering_tolerance[MAX_N_BOUNDARIES], + const CCTK_REAL boundary_extrapolation_tolerance[MAX_N_BOUNDARIES]) +{ +int ibndry; + for (ibndry = 0 ; ibndry < N_boundaries ; ++ibndry) { - /* interp_operator == interp_operator_Lagrange */ - { - /* N_dims = 0 ==> unused */ - { - /* molecule_family = molecule_family_cube */ - { NULL_INTERP_FN_PTR, }, /* order=0 (unused) */ - { NULL_INTERP_FN_PTR, }, /* order=1 (unused) */ - { NULL_INTERP_FN_PTR, }, /* order=2 (unused) */ - { NULL_INTERP_FN_PTR, }, /* order=3 (unused) */ - { NULL_INTERP_FN_PTR, }, /* order=4 (unused) */ - { NULL_INTERP_FN_PTR, }, /* order=5 (unused) */ - { NULL_INTERP_FN_PTR, }, /* order=6 (unused) */ - }, - }, - - { - /* N_dims = 1 */ - { - /* molecule_family = molecule_family_cube */ - { NULL_INTERP_FN_PTR, }, /* order=0 (unused) */ - { LocalInterp_ULagrange_1d_cube10, }, /* order=1 */ - { LocalInterp_ULagrange_1d_cube20, }, /* order=2 */ - { LocalInterp_ULagrange_1d_cube30, }, /* order=3 */ - { LocalInterp_ULagrange_1d_cube40, }, /* order=4 */ - { LocalInterp_ULagrange_1d_cube50, }, /* order=5 */ - { LocalInterp_ULagrange_1d_cube60, }, /* order=6 */ - }, - }, + if ( (boundary_off_centering_tolerance[ibndry] == 0.0) + && (boundary_extrapolation_tolerance[ibndry] > 0.0) ) + then CCTK_VWarn(WARNING_MSG_SEVERITY_LEVEL, + __LINE__, __FILE__, CCTK_THORNSTRING, +"\n" +" CCTK_InterpLocalUniform():\n" +" warning: setting boundary_off_centering_tolerance[%d] == 0.0\n" +" and boundary_extrapolation_tolerance[%d] == %g > 0.0\n" +" is almost certainly a mistake\n" +" (the boundary_extrapolation_tolerance[%d] == %g\n" +" setting will have no effect because of the\n" +" boundary_off_centering_tolerance[%d] == 0.0 setting)" + , + ibndry, + ibndry, + (double) boundary_extrapolation_tolerance[ibndry], + ibndry, + (double) boundary_extrapolation_tolerance[ibndry], + ibndry); + } +} - { - /* N_dims = 2 */ - { - /* molecule_family = molecule_family_cube */ - { NULL_INTERP_FN_PTR, }, /* order=0 (unused) */ - { LocalInterp_ULagrange_2d_cube10, }, /* order=1 */ - { LocalInterp_ULagrange_2d_cube20, }, /* order=2 */ - { LocalInterp_ULagrange_2d_cube30, }, /* order=3 */ - { LocalInterp_ULagrange_2d_cube40, }, /* order=4 */ - { NULL_INTERP_FN_PTR, }, /* order=5 (n.i.) */ - { NULL_INTERP_FN_PTR, }, /* order=6 (n.i.) */ - }, - }, +/******************************************************************************/ - { - /* N_dims = 3 */ - { - /* molecule_family = molecule_family_cube */ - { NULL_INTERP_FN_PTR, }, /* order=0 (unused) */ - { LocalInterp_ULagrange_3d_cube10, }, /* order=1 */ - { LocalInterp_ULagrange_3d_cube20, }, /* order=2 */ - { LocalInterp_ULagrange_3d_cube30, }, /* order=3 */ - { LocalInterp_ULagrange_3d_cube40, }, /* order=4 */ - { NULL_INTERP_FN_PTR, }, /* order=5 (n.i.) */ - { NULL_INTERP_FN_PTR, }, /* order=6 (n.i.) */ - }, - }, +/* + * This function gets + * const char molecule_family[] + * from the parameter table, or sets the default value "cube" if this + * key isn't present. This function also decodes molecule_family into + * an enum molecule_family . + * + * Results: + * This function returns 0 for ok, or the (nonzero) return code for error. + */ +static + int get_and_decode_molecule_family + (int param_table_handle, + int buffer_size, char molecule_family_string_buffer[], + enum molecule_family *p_molecule_family) +{ +enum molecule_family molecule_family; +int status; - /* end of interp_operator == interp_operator_Lagrange */ - }, +status = Util_TableGetString(param_table_handle, + buffer_size, molecule_family_string_buffer, + "molecule_family"); - { - /* interp_operator == interp_operator_Hermite */ - { - /* N_dims = 0 ==> unused */ - { - /* molecule_family = molecule_family_cube */ - { NULL_INTERP_FN_PTR, }, /* order=0 (unused) */ - { NULL_INTERP_FN_PTR, }, /* order=1 (unused) */ - { NULL_INTERP_FN_PTR, }, /* order=2 (unused) */ - { NULL_INTERP_FN_PTR, }, /* order=3 (unused) */ - { NULL_INTERP_FN_PTR, }, /* order=4 (unused) */ - { NULL_INTERP_FN_PTR, }, /* order=5 (unused) */ - { NULL_INTERP_FN_PTR, }, /* order=6 (unused) */ - }, - }, +if (status == UTIL_ERROR_TABLE_NO_SUCH_KEY) + then { + /* set the default value */ + LocalInterp_Strlcpy(molecule_family_string_buffer, "cube", buffer_size); + molecule_family = molecule_family_cube; - { - /* N_dims = 1 */ - { - /* molecule_family = molecule_family_cube */ - { NULL_INTERP_FN_PTR, }, /* order=0 (unused) */ - { NULL_INTERP_FN_PTR, }, /* order=1 (unused) */ - { LocalInterp_UHermite_1d_cube2, }, /* order=2 */ - { LocalInterp_UHermite_1d_cube3, }, /* order=3 */ - { LocalInterp_UHermite_1d_cube4, }, /* order=4 */ - { NULL_INTERP_FN_PTR, }, /* order=5 (n.i.) */ - { NULL_INTERP_FN_PTR, }, /* order=6 (n.i.) */ - }, - }, + /* set this key in the parameter table */ + /* to give the (default) molecule family we're going to use */ + status = Util_TableSetString(param_table_handle, + molecule_family_string_buffer, + "molecule_family"); + if (status < 0) + then { + CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, + __LINE__, __FILE__, CCTK_THORNSTRING, +"\n" +" CCTK_InterpLocalUniform(): error setting default molecule family\n" +" in parameter table!\n" +" error status=%d", + status); + return status; /*** ERROR RETURN ***/ + } + } +else if (status > 0) + then { + /* decode molecule family string */ + if (strcmp(molecule_family_string_buffer, "cube") == 0) + then molecule_family = molecule_family_cube; + else { + CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, + __LINE__, __FILE__, CCTK_THORNSTRING, +"CCTK_InterpLocalUniform(): unknown molecule_family string \"%s\"!", + molecule_family_string_buffer); + return UTIL_ERROR_BAD_INPUT; /*** ERROR RETURN ***/ + } + } +else { + CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, + __LINE__, __FILE__, CCTK_THORNSTRING, +"\n" +" CCTK_InterpLocalUniform(): bad table entry for\n" +" \"molecule_family\" parameter!\n" +" error status=%d", + status); + return status; /*** ERROR RETURN ***/ + } - { - /* N_dims = 2 */ - { - /* molecule_family = molecule_family_cube */ - { NULL_INTERP_FN_PTR, }, /* order=0 (unused) */ - { NULL_INTERP_FN_PTR, }, /* order=1 (unused) */ - { LocalInterp_UHermite_2d_cube2, }, /* order=2 */ - { LocalInterp_UHermite_2d_cube3, }, /* order=3 */ - { NULL_INTERP_FN_PTR, }, /* order=4 (n.i.) */ - { NULL_INTERP_FN_PTR, }, /* order=5 (n.i.) */ - { NULL_INTERP_FN_PTR, }, /* order=6 (n.i.) */ - }, - }, +*p_molecule_family = molecule_family; +return 0; /*** NORMAL RETURN ***/ +} - { - /* N_dims = 3 */ - { - /* molecule_family = molecule_family_cube */ - { NULL_INTERP_FN_PTR, }, /* order=0 (unused) */ - { NULL_INTERP_FN_PTR, }, /* order=1 (unused) */ - { LocalInterp_UHermite_3d_cube2, }, /* order=2 */ - { LocalInterp_UHermite_3d_cube3, }, /* order=3 */ - { NULL_INTERP_FN_PTR, }, /* order=4 (n.i.) */ - { NULL_INTERP_FN_PTR, }, /* order=5 (n.i.) */ - { NULL_INTERP_FN_PTR, }, /* order=6 (n.i.) */ - }, - }, +/******************************************************************************/ - /* end of interp_operator == interp_operator_Hermite */ - } +/* + * This function gets + * CCTK_INT* molecule_positions[N_dims] + * from the parameter table. + * + * Results: + * This function returns 0 for ok, or the (nonzero) return code for error. + */ +static + int get_molecule_positions + (int param_table_handle, + int N_dims, CCTK_INT* molecule_positions_array[MAX_N_DIMS]) +{ +CCTK_POINTER molecule_positions_temp[MAX_N_DIMS]; +int status; - }; +status = Util_TableGetPointerArray(param_table_handle, + N_dims, molecule_positions_temp, + "molecule_positions"); -/* look up the subfunction to do the interpolation */ -p_interp_fn_t p_interp_fn = p_interp_fn_table[interp_operator] - [N_dims] - [molecule_family] - [order] - [smoothing]; -if (p_interp_fn == NULL_INTERP_FN_PTR) +if (status == N_dims) then { - free(operation_codes); - free(operand_indices); - free(input_array_offsets); - if (p_Jacobian_info != NULL) - then { - free(Jacobian_info.Jacobian_pointer); - free(Jacobian_info.Jacobian_offset); + /* type-cast CCTK_POINTER pointers to CCTK_INT* pointers */ + /* (which point to where the query results will be stored) */ + int axis; + for (axis = 0 ; axis < N_dims ; ++axis) + { + molecule_positions_array[axis] + = (CCTK_INT *) molecule_positions_temp[axis]; } + } +else { CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, __LINE__, __FILE__, CCTK_THORNSTRING, "\n" -" CCTK_InterpLocalUniform():\n" -" interpolation not implemented for the combination\n" -" interp_operator=\"%s\", N_dims=%d\n" -" molecule_family=\"%s\", order=%d, smoothing=%d", - interp_operator_string, N_dims, - molecule_family_str, order, smoothing); - return UTIL_ERROR_BAD_INPUT; /*** ERROR RETURN ***/ +" CCTK_InterpLocalUniform(): bad \"molecule_positions\" table entry!\n" +" error status=%d", + status); + return status; /*** ERROR RETURN ***/ } +return 0; /*** NORMAL RETURN ***/ +} + +/******************************************************************************/ + /* - * call the subfunction to actually do the interpolation + * This function gets the Jacobian-query parameters + * CCTK_REAL* Jacobian_pointer[N_output_arrays] + * CCTK_INT Jacobian_offset [N_output_arrays] # optional, default=all 0 + * CCTK_INT Jacobian_interp_point_stride + * CCTK_INT Jacobian_m_point_strides[N_dims] + * CCTK_INT Jacobian_part_stride # optional, default=1 + * and stores them in a struct Jacobian_info . + * + * Results: + * This function returns 0 for ok, or the (nonzero) return code for error. */ +static + int get_Jacobian_info(int param_table_handle, + int N_dims, int N_output_arrays, + struct Jacobian_info* p_Jacobian_info) +{ +/* padded array size, cf. LocalInterp_InterpLocalUniform() header comments */ +const int N_output_arrays1 = N_output_arrays + 1; + +int status; + +p_Jacobian_info->Jacobian_pointer + = (CCTK_REAL**) malloc(N_output_arrays1 * sizeof(CCTK_REAL *)); +if (p_Jacobian_info->Jacobian_pointer == NULL) + then return UTIL_ERROR_NO_MEMORY; /*** ERROR RETURN ***/ + { -struct error_flags error_flags; -struct molecule_structure_flags molecule_structure_flags; -const int return_code - = (*p_interp_fn)(coord_origin, coord_delta, - N_interp_points, - interp_coords_type_code, interp_coords, - N_boundary_points_to_omit, - boundary_off_centering_tolerance, - boundary_extrapolation_tolerance, - N_input_arrays, - input_array_offsets, input_array_strides, - input_array_min_subscripts, - input_array_max_subscripts, - input_array_type_codes, input_arrays, - N_output_arrays, - output_array_type_codes, output_arrays, - operand_indices, operation_codes, - &error_flags, - &molecule_structure_flags, - p_molecule_min_max_m_info, - p_molecule_positions, - p_Jacobian_info); +CCTK_POINTER* Jacobian_pointer_temp + = (CCTK_POINTER*) malloc(N_output_arrays1 * sizeof(CCTK_POINTER)); +if (Jacobian_pointer_temp == NULL) + then { + free(p_Jacobian_info->Jacobian_pointer); + return UTIL_ERROR_NO_MEMORY; /*** ERROR RETURN ***/ + } +status = Util_TableGetPointerArray(param_table_handle, + N_output_arrays, Jacobian_pointer_temp, + "Jacobian_pointer"); +if (status == N_output_arrays) + then { + /* type-cast CCTK_POINTER pointers to CCTK_REAL* pointers */ + /* (which point to where the query results will be stored) */ + int out; + for (out = 0 ; out < N_output_arrays ; ++out) + { + p_Jacobian_info->Jacobian_pointer[out] + = (CCTK_REAL *) Jacobian_pointer_temp[out]; + } + free(Jacobian_pointer_temp); + } +else { + free(p_Jacobian_info->Jacobian_pointer); + free(Jacobian_pointer_temp); + CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, + __LINE__, __FILE__, CCTK_THORNSTRING, +"\n" +" CCTK_InterpLocalUniform(): bad \"Jacobian_pointer\" table entry!\n" +" error status=%d" + , + status); + return status; /*** ERROR RETURN ***/ + } + } + +p_Jacobian_info->Jacobian_offset + = (CCTK_INT*) malloc(N_output_arrays1 * sizeof(CCTK_INT)); +if (p_Jacobian_info->Jacobian_offset == NULL) + then { + free(p_Jacobian_info->Jacobian_pointer); + return UTIL_ERROR_NO_MEMORY; /*** ERROR RETURN ***/ + } +status = get_INT_array(param_table_handle, "Jacobian_offset", + true, 0, /* default value */ + N_output_arrays, p_Jacobian_info->Jacobian_offset, + NULL); +if (status != 0) + then { + free(p_Jacobian_info->Jacobian_pointer); + free(p_Jacobian_info->Jacobian_offset); + return status; /*** ERROR RETURN ***/ + } + +status = get_and_check_INT(param_table_handle, "Jacobian_interp_point_stride", + true, 0, /* this is a mandatory parameter */ + /* if we're executing this function */ + false, 0, 0, NULL, /* no range check */ + & p_Jacobian_info->Jacobian_interp_point_stride); +if (status != 0) + then { + free(p_Jacobian_info->Jacobian_offset); + free(p_Jacobian_info->Jacobian_pointer); + return status; /*** ERROR RETURN ***/ + } + +status = get_INT_array(param_table_handle, "Jacobian_m_strides", + false, 0, /* don't set default value */ + N_dims, p_Jacobian_info->Jacobian_m_strides, + NULL); /* this is a mandatory parameter */ + /* if we're executing this function */ +if (status != 0) + then { + free(p_Jacobian_info->Jacobian_pointer); + free(p_Jacobian_info->Jacobian_offset); + return status; /*** ERROR RETURN ***/ + } + +status = get_and_check_INT(param_table_handle, "Jacobian_part_stride", + false, 1, /* default value */ + false, 0, 0, NULL, /* no range check */ + & p_Jacobian_info->Jacobian_part_stride); +if (status != 0) + then { + free(p_Jacobian_info->Jacobian_pointer); + free(p_Jacobian_info->Jacobian_offset); + return status; /*** ERROR RETURN ***/ + } + +return 0; /*** NORMAL RETURN ***/ +} /******************************************************************************/ /* - * is an interpolation point out of range? + * This function sets the error-info entries + * CCTK_INT error_pt + * CCTK_INT error_ibndry + * CCTK_INT error_axis + * CCTK_INT error_direction + * in the parameter table. + * + * Results: + * This function returns 0 for ok, or the (nonzero) return code for error. */ -if (return_code == CCTK_ERROR_INTERP_POINT_X_RANGE) +static + int set_error_info(int param_table_handle, + struct error_flags* p_error_flags) +{ +const int status1 = Util_TableSetInt(param_table_handle, + p_error_flags->error_pt, + "error_pt"); +const int status2 = Util_TableSetInt(param_table_handle, + p_error_flags->error_ibndry, + "error_ibndry"); +const int status3 = Util_TableSetInt(param_table_handle, + p_error_flags->error_axis, + "error_axis"); +const int status4 = Util_TableSetInt(param_table_handle, + p_error_flags->error_direction, + "error_direction"); +if ((status1 < 0) || (status2 < 0) || (status3 < 0) || (status4 < 0)) then { - /* store details about the error in the parameter table */ - - status1 = Util_TableSetInt(param_table_handle, - error_flags.error_pt, - "out_of_range_pt"); - status2 = Util_TableSetInt(param_table_handle, - error_flags.error_axis, - "out_of_range_axis"); - status3 = Util_TableSetInt(param_table_handle, - error_flags.error_end, - "out_of_range_end"); - if ((status1 < 0) || (status2 < 0) || (status3 < 0)) - then { - free(operation_codes); - free(operand_indices); - free(input_array_offsets); - if (p_Jacobian_info != NULL) - then { - free(Jacobian_info.Jacobian_pointer); - free(Jacobian_info.Jacobian_offset); - } - CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, - __LINE__, __FILE__, CCTK_THORNSTRING, + CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, + __LINE__, __FILE__, CCTK_THORNSTRING, "\n" -" CCTK_InterpLocalUniform(): error setting out-of-range table entry/entries!" - " error status1=%d status2=%d status3=%d", - status1, status2, status3); - return (status1 < 0) ? status1 /*** ERROR RETURN ***/ - : (status2 < 0) ? status2 /*** ERROR RETURN ***/ - : status3; /*** ERROR RETURN ***/ - } +" CCTK_InterpLocalUniform():\n" +" error setting error-info table entry/entries!\n" +" error status1=%d status2=%d status3=%d status4=%d" + , + status1, status2, status3, status4); + return (status1 < 0) ? status1 /*** ERROR RETURN ***/ + : (status2 < 0) ? status2 /*** ERROR RETURN ***/ + : (status3 < 0) ? status3 /*** ERROR RETURN ***/ + : status4; /*** ERROR RETURN ***/ } +return 0; /*** NORMAL RETURN ***/ +} + /******************************************************************************/ /* - * store molecule structure info + * This function sets the molecule structure flags + * CCTK_INT MSS_is_fn_of_interp_coords + * CCTK_INT MSS_is_fn_of_which_operation + * CCTK_INT MSS_is_fn_of_input_array_values + * CCTK_INT Jacobian_is_fn_of_input_array_values + * in the parameter table. + * + * Results: + * This function returns 0 for ok, or the (nonzero) return code for error. */ -status1 = Util_TableSetInt(param_table_handle, - molecule_structure_flags.MSS_is_fn_of_interp_coords, - "MSS_is_fn_of_interp_coords"); -status2 = Util_TableSetInt(param_table_handle, - molecule_structure_flags - .MSS_is_fn_of_which_operation, - "MSS_is_fn_of_which_operation"); -status3 = Util_TableSetInt(param_table_handle, - molecule_structure_flags - .MSS_is_fn_of_input_array_values, - "MSS_is_fn_of_input_array_values"); -status4 = Util_TableSetInt(param_table_handle, - molecule_structure_flags - .Jacobian_is_fn_of_input_array_values, - "Jacobian_is_fn_of_input_array_values"); +static + int set_molecule_structure + (int param_table_handle, + const struct molecule_structure_flags* p_molecule_structure_flags) +{ +const int status1 + = Util_TableSetInt(param_table_handle, + p_molecule_structure_flags->MSS_is_fn_of_interp_coords, + "MSS_is_fn_of_interp_coords"); +const int status2 + = Util_TableSetInt(param_table_handle, + p_molecule_structure_flags->MSS_is_fn_of_which_operation, + "MSS_is_fn_of_which_operation"); +const int status3 + = Util_TableSetInt(param_table_handle, + p_molecule_structure_flags->MSS_is_fn_of_input_array_values, + "MSS_is_fn_of_input_array_values"); +const int status4 + = Util_TableSetInt(param_table_handle, + p_molecule_structure_flags->Jacobian_is_fn_of_input_array_values, + "Jacobian_is_fn_of_input_array_values"); if ((status1 < 0) || (status2 < 0) || (status3 < 0) || (status4 < 0)) then { - free(operation_codes); - free(operand_indices); - free(input_array_offsets); - if (p_Jacobian_info != NULL) - then { - free(Jacobian_info.Jacobian_pointer); - free(Jacobian_info.Jacobian_offset); - } CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, __LINE__, __FILE__, CCTK_THORNSTRING, "\n" -" CCTK_InterpLocalUniform(): error setting Jacobian structure\n" -" table entry/entries!" -" error status1=%d status2=%d\n" -" status3=%d status4=%d" +" CCTK_InterpLocalUniform():\n" +" error setting molecule structure flags table entry/entries!" +" error status1=%d status2=%d status3=%d status4=%d" , status1, status2, status3, status4); return (status1 < 0) ? status1 /*** ERROR RETURN ***/ @@ -1771,72 +1883,289 @@ if ((status1 < 0) || (status2 < 0) || (status3 < 0) || (status4 < 0)) : status4; /*** ERROR RETURN ***/ } +return 0; /*** NORMAL RETURN ***/ +} + /******************************************************************************/ /* - * store results from molecule min/max m queries + * This function sets the molecule size parameters + * CCTK_INT molecule_min_m[N_dims] + * CCTK_INT molecule_max_m[N_dims] + * in the parameter table. * - * ... molecule position queries are stored directly into - * their final destination by subfunction - * ==> no further store or cleanup needed here + * Results: + * This function returns 0 for ok, or the (nonzero) return code for error. */ +static + int set_molecule_min_max_m + (int param_table_handle, + int N_dims, + const struct molecule_min_max_m_info* p_molecule_min_max_m_info) +{ +const int status1 + = Util_TableSetIntArray(param_table_handle, + N_dims, p_molecule_min_max_m_info->molecule_min_m, + "molecule_min_m"); +const int status2 + = Util_TableSetIntArray(param_table_handle, + N_dims, p_molecule_min_max_m_info->molecule_max_m, + "molecule_max_m"); +if ((status1 < 0) || (status2 < 0)) + then { + CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, + __LINE__, __FILE__, CCTK_THORNSTRING, +"\n" +" CCTK_InterpLocalUniform(): error setting\n" +" \"molecule_{min,max}_m\" table entry/entries!" +" error status1=%d status2=%d", + status1, status2); + return (status1 < 0) ? status1 : status2; /*** ERROR RETURN ***/ + } -if (p_molecule_min_max_m_info != NULL) +return 0; /*** NORMAL RETURN ***/ +} + +/******************************************************************************/ +/******************************************************************************/ +/******************************************************************************/ + +/* + * This function gets and range-checks a CCTK_INT from the parameter table. + * + * Arguments: + * param_table_handle = handle to the Cactus key-value table + * name = the character-string key in the table + * mandatory_flag = true ==> the value must be present in the parameter + * table (default_value is ignored) + * false ==> the value is optional + * default_value = the default value (for each array element) if the + * key isn't in the parameter table and mandatory_flag == false + * check_range_flag = true ==> check that the value satisfies + * min_value <= value <= max_value + * false ==> don't do a range check, and ignore + * min_value, max_value, and max_value_string + * {min,max}_value = the inclusive range of valid values + * max_value_string = the character-string name of max_value (this is used + * only in formatting the out-of-range error message) + * p_value --> where we should store the value + * + * Results: + * This function returns 0 for ok, or the (nonzero) return code for error. + */ +static + int get_and_check_INT + (int param_table_handle, const char name[], + bool mandatory_flag, int default_value, + bool check_range_flag, int min_value, int max_value, + const char max_value_string[], + CCTK_INT* p_value) +{ +CCTK_INT value; + +const int status = Util_TableGetInt(param_table_handle, &value, name); + +if (status == 1) + then { } /* value set from table ==> no-op here */ +else if (status == UTIL_ERROR_TABLE_NO_SUCH_KEY) then { - status1 = Util_TableSetIntArray(param_table_handle, - N_dims, p_molecule_min_max_m_info - ->molecule_min_m, - "molecule_min_m"); - status2 = Util_TableSetIntArray(param_table_handle, - N_dims, p_molecule_min_max_m_info - ->molecule_max_m, - "molecule_max_m"); - if ((status1 < 0) || (status2 < 0)) + if (mandatory_flag) then { - free(operation_codes); - free(operand_indices); - free(input_array_offsets); - if (p_Jacobian_info != NULL) - then { - free(Jacobian_info.Jacobian_pointer); - free(Jacobian_info.Jacobian_offset); - } CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, __LINE__, __FILE__, CCTK_THORNSTRING, "\n" -" CCTK_InterpLocalUniform(): error setting\n" -" \"molecule_{min,max}_m\" table entry/entries!" -" error status1=%d status2=%d", - status1, status2); - return (status1 < 0) ? status1 : status2; /*** ERROR RETURN ***/ +" CCTK_InterpLocalUniform(): missing table entry for\n" +" \"%s\" parameter!\n" +" (this is a mandatory parameter)!\n" + , + name); + return status; /*** ERROR RETURN ***/ } + else value = default_value; + } +else { + CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, + __LINE__, __FILE__, CCTK_THORNSTRING, +"\n" +" CCTK_InterpLocalUniform(): bad table entry for\n" +" \"%s\" parameter!\n" +" error status=%d" + , + name, + status); + return status; /*** ERROR RETURN ***/ } +if (check_range_flag) + then { + if ((value < min_value) || (value > max_value)) + then { + CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, + __LINE__, __FILE__, CCTK_THORNSTRING, +"\n" +" CCTK_InterpLocalUniform(): %s=%d is invalid!\n" +" valid range is %d <= %s <= %s=%d", + name, value, + min_value, name, max_value_string, max_value); + return UTIL_ERROR_BAD_INPUT; /*** ERROR RETURN ***/ + } + } + +*p_value = value; +return 0; /*** NORMAL RETURN ***/ +} + /******************************************************************************/ -free(operation_codes); -free(operand_indices); -free(input_array_offsets); -if (p_Jacobian_info != NULL) +/* + * This function gets an array of CCTK_INT values from the parameter table. + * + * Arguments: + * param_table_handle = handle to the Cactus key-value table + * name = the character-string key in the table + * default_flag = false ==> ignore default_value and set + * *p_value_not_set to true if the + * key is not present in the parameter table, + * false if it is present + * true ==> if the key is not present in the parameter table, + * set each array element to default_value + * default_value = the default value (for each array element) if the + * key isn't in the parameter table + * N = the size of the array + * buffer --> a buffer in which to store the array + * p_value_not_set = NULL or a pointer to a flag to be set if and only if + * this function did *not* set values in the buffer. + * More precisely, if default_flag is set, this argument + * is ignored. Otherwise, the logic is this: + * if (this function set values in the buffer) + * then { + * if (p_value_not_set is non-NULL) + * then *p_value_not_set = false + * } + * else { + * if (p_value_not_set is non-NULL) + * then *p_value_not_set = true + * else give an error message that + * a mandatory value is missing + * from the parameter table, then + * return UTIL_ERROR_TABLE_NO_SUCH_KEY + * } + * + * Results: + * This function returns 0 for ok, or the desired (nonzero) error return + * code if something is wrong. Note that error return codes may be either + * positive or negative! + */ +static + int get_INT_array(int param_table_handle, const char name[], + bool default_flag, int default_value, + int N, CCTK_INT buffer[], + bool* p_value_not_set) +{ +const int status + = Util_TableGetIntArray(param_table_handle, + N, buffer, + name); + +if (status == UTIL_ERROR_TABLE_NO_SUCH_KEY) then { - free(Jacobian_info.Jacobian_pointer); - free(Jacobian_info.Jacobian_offset); + /* we didn't set values in the buffer */ + /* ==> either set the default value ourself, */ + /* or print an error message, */ + /* or notify our caller that we didn't do so */ + if (default_flag) + then { + int i; + for (i = 0 ; i < N ; ++i) + { + buffer[i] = default_value; + } + } + else { + if (p_value_not_set == NULL) + then { + CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, + __LINE__, __FILE__, CCTK_THORNSTRING, +"\n" +" CCTK_InterpLocalUniform(): missing table entry for\n" +" \"%s\" parameter!\n" +" (this is a mandatory parameter)!\n" + , + name); + return UTIL_ERROR_TABLE_NO_SUCH_KEY; /*** ERROR RETURN ***/ + } + else *p_value_not_set = true; + } + } +else if (status == N) + then { + /* we did set values in the buffer */ + if (! default_flag) + then if (p_value_not_set != NULL) + then *p_value_not_set = false; + } +else { + CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, + __LINE__, __FILE__, CCTK_THORNSTRING, +"\n" +" CCTK_InterpLocalUniform(): bad table entry for\n" +" \"%s\" parameter!\n" +" error status=%d", + name, status); + return status; /*** ERROR RETURN ***/ } -return return_code; - } - } - } - } - } - } - } - } - } - } - } - } - } - } - } +return 0; /*** NORMAL RETURN ***/ +} + +/******************************************************************************/ + +/* + * This function gets an array of CCTK_REAL values from the parameter table. + * + * Arguments: + * param_table_handle = handle to the Cactus key-value table + * name = the character-string key in the table + * default_value = the default value (for each array element) if the + * key isn't in the parameter table + * N = the size of the array + * buffer --> a buffer in which to store the array + * + * Results: + * This function returns 0 for ok, or the desired (nonzero) error return + * code if something is wrong. Note that error return codes may be either + * positive or negative! + */ +static + int get_REAL_array(int param_table_handle, const char name[], + CCTK_REAL default_value, + int N, CCTK_REAL buffer[]) +{ +const int status + = Util_TableGetRealArray(param_table_handle, + N, buffer, + name); +if (status == UTIL_ERROR_TABLE_NO_SUCH_KEY) + then { + /* set default value */ + int i; + for (i = 0 ; i < N ; ++i) + { + buffer[i] = default_value; + } + } +else if (status == N) + then { } /* ok ==> no-op here */ +else { + CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, + __LINE__, __FILE__, CCTK_THORNSTRING, +"\n" +" CCTK_InterpLocalUniform(): bad table entry for\n" +" \"%s\" parameter!\n" +" error status=%d", + name, status); + return status; /*** ERROR RETURN ***/ + } + +return 0; /*** NORMAL RETURN ***/ } diff --git a/src/GeneralizedPolynomial-Uniform/InterpLocalUniform.h b/src/GeneralizedPolynomial-Uniform/InterpLocalUniform.h index 851378a..8d9c551 100644 --- a/src/GeneralizedPolynomial-Uniform/InterpLocalUniform.h +++ b/src/GeneralizedPolynomial-Uniform/InterpLocalUniform.h @@ -109,15 +109,26 @@ enum molecule_family * other compile-time settings */ -/* default for each element of boundary_off_centering_tolerance[] */ -#define BOUNDARY_OFF_CENTER_TOL_DEFAULT 999.0 - -/* default for each element of boundary_extrapolation_tolerance[] */ -#define BOUNDARY_EXTRAP_TOL_DEFAULT 1.0e-10 +/* defaults for boundary_{off_centering,extrapolation}_tolerance[] */ +#ifdef NOT_YET + #define LAGRANGE_BNDRY_OFF_CNTR_TOL_DEF 999.0 + #define LAGRANGE_BNDRY_EXTRAP_TOL_DEF 1.0e-10 + #define HERMITE_BNDRY_OFF_CNTR_TOL_DEF 1.0e-10 + #define HERMITE_BNDRY_EXTRAP_TOL_DEF 0.0 +#else + #define LAGRANGE_BNDRY_OFF_CNTR_TOL_DEF 999.0 + #define LAGRANGE_BNDRY_EXTRAP_TOL_DEF 1.0e-10 + #define HERMITE_BNDRY_OFF_CNTR_TOL_DEF 999.0 + #define HERMITE_BNDRY_EXTRAP_TOL_DEF 1.0e-10 +#endif /* CCTK_VWarn() severity level for error/warning messages */ -#define ERROR_MSG_SEVERITY_LEVEL 0 -#define WARNING_MSG_SEVERITY_LEVEL 1 +#define BUG_MSG_SEVERITY_LEVEL 0 +#define ERROR_MSG_SEVERITY_LEVEL 0 +#define WARNING_MSG_SEVERITY_LEVEL 1 + +/* CCTK_Abort() exit code for internal error (interpolator bug) aborts */ +#define BUG_ABORT_CODE 42 /******************************************************************************/ @@ -133,7 +144,7 @@ struct error_flags int error_pt; int error_ibndry; int error_axis; - int error_end; + int error_direction; }; struct molecule_structure_flags @@ -230,5 +241,5 @@ int LocalInterp_molecule_posn(fp grid_origin, fp grid_delta, int* i_center, fp* x_rel); /* functions in util.c */ -void LocalInterp_zero_int_array(int N, CCTK_INT array[]); int LocalInterp_decode_N_parts(int type_code); +size_t LocalInterp_Strlcpy(char* dst, const char* src, size_t dst_size); -- cgit v1.2.3