/*@@ @file InterpLocalUniform.c @date 22 Oct 2001 @author Jonathan Thornburg @desc Generalized Interpolation of processor-local arrays. This code interpolates a set of functions defined on a uniform N-dimensional grid, to a set of interpolation points. It can also do differentiation and/or smoothing in the process of the interpolation. Conceptually, the generalized interpolation is done by least-squares fitting a polynomial of the specified order to the data values, applying the differentiation (if any) to it, then evaluating the result at the interpolation points. Since this ultimately yields _some_ linear combination of the data values, we precompute the coefficients for efficiency. This is done with separate Maple-generated formulas for each combination of number of interpolation operator, number of dimensions, molecule family, order, amount of Savitzky-Golay smoothing to be done, and differentiation operator. @enddesc @version $Id$ @@*/ #include #include #include /* malloc(), free() */ #include #include "cctk.h" #include "util_ErrorCodes.h" #include "util_Table.h" #include "InterpLocalUniform.h" #include "Lagrange/all_prototypes.h" #include "Hermite/all_prototypes.h" /* the rcs ID and its dummy function to use it */ static const char *rcsid = "$Header$"; CCTK_FILEVERSION(CactusBase_LocalInterp_src_GeneralizedPolynomialUniform_InterpLocalUniform_c) /******************************************************************************/ /* * *****data structures and functions local to this file ***** */ /* * data structures local to this file */ /* this enum describes which interpolation operator we're using */ enum interp_operator { interp_operator_Lagrange, /* "Lagrange polynomial interpolation" */ interp_operator_Hermite, /* "Hermite polynomial interpolation" */ N_INTERP_OPERATORS /* this must be the last entry in the enum */ }; /* * prototypes for functions local to this file */ static int LocalInterp_InterpLocalUniform(enum interp_operator interp_operator, const char* interp_operator_string, int N_dims, 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[]); /******************************************************************************/ /******************************************************************************/ /******************************************************************************/ /*@@ @routine LocalInterp_ILU_Lagrange @date 22 Oct 2001 @author Jonathan Thornburg @desc This function is a top-level driver for LocalInterp::CCTK_InterpLocalUniform(). It does Lagrange interpolation of a set of input arrays defined on a uniform N-dimensional tensor-product grid, to a set of interpolation points. It can also do differentiation and/or smoothing in the process of the interpolation. It can also return information about the Jacobian of the interpolated values with respect to the input arrays. This function does nothing except pass all its arguments down to LocalInterp_InterpLocalUniform() with an extra 2 arguments added to indicate the operator handle. See that function for details on this function's arguments/results. @enddesc @@*/ int LocalInterp_ILU_Lagrange(int N_dims, 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[]) { return LocalInterp_InterpLocalUniform(interp_operator_Lagrange, "Lagrange polynomial interpolation", N_dims, param_table_handle, coord_origin, coord_delta, N_interp_points, interp_coords_type_code, interp_coords, N_input_arrays, input_array_dims, input_array_type_codes, input_arrays, N_output_arrays, output_array_type_codes, output_arrays); } /******************************************************************************/ /*@@ @routine LocalInterp_ILU_Hermite @date 22 Oct 2001 @author Jonathan Thornburg @desc This function is a top-level driver for LocalInterp::CCTK_InterpLocalUniform(). It does Hermite interpolation of a set of input arrays defined on a uniform N-dimensional tensor-product grid, to a set of interpolation points. It can also do differentiation and/or smoothing in the process of the interpolation. It can also return information about the Jacobian of the interpolated values with respect to the input arrays. This function does nothing except pass all its arguments down to LocalInterp_InterpLocalUniform() with an extra 2 arguments added to indicate the operator handle. See that function for details on this function's arguments/results. @enddesc @@*/ int LocalInterp_ILU_Hermite(int N_dims, 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[]) { return LocalInterp_InterpLocalUniform(interp_operator_Hermite, "Hermite polynomial interpolation", N_dims, param_table_handle, coord_origin, coord_delta, N_interp_points, interp_coords_type_code, interp_coords, N_input_arrays, input_array_dims, input_array_type_codes, input_arrays, N_output_arrays, output_array_type_codes, output_arrays); } /******************************************************************************/ /******************************************************************************/ /******************************************************************************/ /*@@ @routine LocalInterp_InterpLocalUniform @date 22 Oct 2001 @author Jonathan Thornburg @desc This function is the main driver for LocalInterp::CCTK_InterpLocalUniform(). It interpolates a set of input arrays defined on a uniform N-dimensional tensor-product grid, to a set of interpolation points. It can also do differentiation and/or smoothing in the process of the interpolation. It can also return information about the Jacobian of the interpolated values with respect to the input arrays. This function is just a driver: it validates arguments, extracts optional arguments from the parameter table, then decodes (interp_operator, N_dims, molecule_family, order, smoothing) and calls the appropriate subfunction to do the actual interpolation, then finally stores some results back in the parameter table. @enddesc ***** misc arguments ***** @var interp_operator @vdesc describes the interpolation operator @vtype enum interp_operator interp_operator @endvar @var interp_operator_string @vdesc gives the character-string name of the interpolation operator (this is used only for printing error messages) @vtype const char* interp_operator_string @endvar @var N_dims @vdesc dimensionality of the interpolation @vtype int N_dims (must be >= 1) @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 ***** input arguments passed in the parameter table ***** @var order @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 @endvar @var out_of_range_tolerance @vdesc Specifies how out-of-range interpolation points should be handled. The array elements are matched up with the axes and minimum/maximum ends of the grid in the order [x_min, x_max, y_min, y_max, z_min, z_max, ...]. An array value TOL is interpreted as follows: If TOL >= 0.0, then an interpolation point is considered to be "out of range" if and only if the interpolation point is > TOL * coord_delta[axis] outside the grid in this coordinate direction. If TOL == -1.0, then an interpolation point is considered to be "out of range" if and only if a centered molecule (or more generally, a molecule whose centering is chosen pretending that the grid is of infinite extent), would require data from outside the grid in this direction. Other values of TOL are illegal. @vtype const CCTK_REAL out_of_range_tolerance[2*N_dims] @endvar @var input_array_offsets @vdesc If this key is present, the value should be a pointer to an array giving an "offset" for each input array, use in the subscripting computation: for input array number in, this computation (given for 3D; the generalization to other numbers of dimensions is obvious) is input_arrays[in][offset + i*istride + j*jstride + k*kstride] where offset = input_array_offsets[in] or is 0 if input_array_offsets is absent (istride,jstride,kstride,...) = input_array_stride[] and where (i,j,k,...) run from input_array_min_subscripts[] to input_array_max_subscripts[] inclusive. @vtype const CCTK_INT input_array_offsets[N_input_arrays] @endvar @var input_array_strides @vdesc (pointer to) array giving strides of input arrays (this is shared by all input arrays) @vtype const CCTK_INT input_array_strides[N_dims] @endvar @var input_array_min_subscripts @vdesc (pointer to) array giving minimum subscripts of input arrays (this is shared by all input arrays) @vtype const CCTK_INT input_array_min_subscripts[N_dims] @endvar @var input_array_max_subscripts @vdesc (pointer to) array giving maximum subscripts of input arrays (this is shared by all input arrays) @vtype const CCTK_INT input_array_max_subscripts[N_dims] @endvar @var operand_indices @vdesc (pointer to) array of integer operand indices specifying which input array (0-origin indexing into input_arrays[]) is to be generalized-interpolated to produce this output @vtype const CCTK_INT operand_indices[N_output_arrays] @endvar @var operation_codes @vdesc (pointer to) array of integer operation codes specifying what (if any) derivatives are to be takin in the interpolation process: 0=no derivative, 1=d/dx, 2=d/dy, 3=d/dz (for coords (x,y,z)), 11=d^2/dx^2, 12=21=d^2/dxdy, etc etc. @vtype const CCTK_INT operation_codes[N_output_arrays] @endvar ***** output arguments passed in the parameter table ***** @var out_of_range_pt @vdesc If an out-of-range point is detected, this table entry is set to its pt value. @vtype CCTK_INT out_of_range_pt @vio out @endvar @var out_of_range_axis @vdesc If an out-of-range point is detected, this table entry is set to the axis value which is out of range. @vtype CCTK_INT out_of_range_axis @vio out @endvar @var out_of_range_end @vdesc If an out-of-range point is detected, this table entry is set to -1/+1 if the point is out of range on the min/max end of the axis. @vtype CCTK_INT out_of_range_end @vio out @endvar @var MSS_is_fn_of_interp_coords @vdesc If the interpolation molecule size and/or shape vary with the interpolation coordinates, this table entry is set to 1. Otherwise (i.e. if the interpolation molecule size and shape are independent of the interpolation coordinates), it's set to 0. @vtype CCTK_INT MSS_is_fn_of_interp_coords @vio out @endvar @var MSS_is_fn_of_which_operation @vdesc If the interpolator supports computing derivatives, and the interpolation molecule size and/or shape vary from one operation_code[] value to another, this table entry is set to 1. Otherwise (i.e. if the interpolator doesn't support computing derivatives, or if the interpolator does support computing derivatives but the interpolation molecule size and shape are independent of the operation_code[] values), it's set to 0. Note that this query tests whether the molecule size and/or shape depend on operation_codes[] in general, independent of whether there are in fact any distinct values (or even any values at all) passed in operation_codes[] in this particular interpolator call. In other words, this is a query about the basic design of this interpolator, not about this particular call. @vtype CCTK_INT MSS_is_fn_of_which_operation @vio out @endvar @var MSS_is_fn_of_input_array_values @vdesc If the interpolation molecule size and/or shape vary with the actual floating-point values of the input arrays, this table entry is set to 1. Otherwise (i.e. if the interpolation molecule size and shape are independent of the input array values; this is a necessary, but not sufficient, condition for the interpolation to be linear), it's set to 0. @vtype CCTK_INT MSS_is_fn_of_input_array_values @vio out @endvar @var Jacobian_is_fn_of_input_array_values @vdesc If the actual floating-point values of the Jacobian (*) (see the discussion of Jacobian_pointer below) depends on the actual floating-point values of the input arrays (i.e. if the interpolation is nonlinear), this table entry is set to 1. Otherwise (i.e. if the interpolation is linear), it's set to 0. @vtype CCTK_INT Jacobian_is_fn_of_input_array_values @vio out @endvar @var molecule_min_m @vdesc If this key is present (the value doesn't matter), then the value is (re)set to an array of N_dims integers giving the minimum molecule m coordinate in each axis @vtype CCTK_INT molecule_min_m[N_dims] @vio out @endvar @var molecule_max_m @vdesc If this key is present (the value doesn't matter), then the value is (re)set to an array of N_dims integers giving the maximum molecule m coordinate in each axis @vtype CCTK_INT molecule_max_m[N_dims] @vio out @endvar @var molecule_positions @vdesc If this key is present, then the value should be an array of N_dims pointers to (caller-supplied) arrays of N_interp_points CCTK_INTs each; the interpolator will store the molecule positions in the pointed-to arrays. @vtype CCTK_INT *const molecule_positions[N_dims] @vio out @endvar @var Jacobian_pointer @vdesc If this key is present, then the value should be an array of N_output_arrays pointers to (caller-supplied) arrays of CCTK_INTs; the interpolator will store the Jacobian elements in the pointed-to arrays. For output array number out, the subscripting computation (given for 3D; the generalization to other numbers of dimensions is obvious) is Jacobian_pointer[out][offset + pt*Jacobian_interp_point_stride + mi*stride_i + mj*stride_j + mk*stride_k + part*Jacobian_part_stride] where offset = Jacobian_offset[out] or is 0 if Jacobian_offset[] is absent pt = the point number (stride_i,stride_j,stride_k) = Jacobian_m_strides[] part = 0 for real values, 0 for the real parts of complex values, 1 for the imaginary parts of complex values, or 0 if Jacobian_part_stride is absent @vtype CCTK_REAL *const Jacobian_pointer[N_output_arrays] @endvar @var Jacobian_offset @vdesc If this key is present, then the value should be a pointer to an array of N_output_arrays CCTK_INTs giving offsets to use in the Jacobian subscripting computation. If this key is absent, the effect is as if a default array of all 0s had been supplied. @vtype const CCTK_INT Jacobian_offset[N_output_arrays] @endvar @var Jacobian_interp_point_stride @vdesc The pt-stride for the Jacobian subscripting computation. @vtype const CCTK_INT Jacobian_interp_point_stride @endvar @var Jacobian_m_strides @vdesc An array of N_dims CCTK_INTs giving the m strides along each axis for the Jacobian subscripting computation @vtype const CCTK_INT Jacobian_m_strides[N_dims] @endvar @var Jacobian_part_stride @vdesc If this key is present, it gives the part stride for the Jacobian subscripting computation. If this key is absent, the default value 0 is used (n.b. this is suitable only for real data arrays) @vtype const CCTK_INT Jacobian_m_strides[N_dims] @endvar ***** return result ***** @returntype int @returndesc 0 successful interpolation 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 interpolation point is out of range (in this case additional information is reported through the parameter table) or any error code returned by one of the Util_TableGet*() functions called by this function @endreturndesc @@*/ static int LocalInterp_InterpLocalUniform(enum interp_operator interp_operator, const char* interp_operator_string, int N_dims, 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[]) { /* * Implementation Note: * * We malloc() several scratch arrays, some with sizes determined by * N_{input,output}_arrays. Thus if N_{input,output}_arrays == 0, with * the obvious code we would malloc(0). Alas, the C standard permits * malloc(0) to return a NULL pointer, which the usual malloc() idiom * CCTK_INT *const p = malloc(N * sizeof(CCTK_INT)); * if (p == NULL) * then return UTIL_ERROR_NO_MEMORY * would falsely detect as an out-of-memory condition. * * As a workaround, we pad all our malloc request sizes, i.e. * CCTK_INT *const p = malloc((N+1) * sizeof(CCTK_INT)); * if (p == NULL) * then return UTIL_ERROR_NO_MEMORY * This is a kludge, but so are the other plausible solutions. :( :( */ int N_input_arrays1 = N_input_arrays + 1; int N_output_arrays1 = N_output_arrays + 1; int status, status1, status2, status3, status4; /******************************************************************************/ /* * basic sanity checks on input parameters */ if ( (N_dims <= 0) || (param_table_handle < 0) || ((N_interp_points > 0) && (coord_origin == NULL)) || ((N_interp_points > 0) && (coord_delta == NULL)) || (N_interp_points < 0) || ((N_interp_points > 0) && (interp_coords == NULL)) || (N_input_arrays < 0) /* no check here on input_array_dims, */ /* since it may be NULL if overridden by parameter-table stuff */ || ((N_input_arrays > 0) && (input_array_type_codes == NULL)) || ((N_input_arrays > 0) && (input_arrays == NULL)) || (N_output_arrays < 0) || ((N_output_arrays > 0) && (output_array_type_codes == NULL)) || ((N_output_arrays > 0) && (output_arrays == NULL)) ) then { CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, __LINE__, __FILE__, CCTK_THORNSTRING, "InterpLocalUniform(): invalid arguments\n" " (N_dims <= 0, param_table_handle < 0, N_input_arrays < 0,\n" " N_output_arrays < 0, and/or NULL pointers-that-shouldn't-be-NULL)!"); return UTIL_ERROR_BAD_INPUT; /*** ERROR RETURN ***/ } if (N_dims > MAX_N_DIMS) then { CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, __LINE__, __FILE__, CCTK_THORNSTRING, "\n" " CCTK_InterpLocalUniform(): implementation restriction: N_dims=%d\n" " but we only support N_dims <= MAX_N_DIMS=%d!", N_dims, MAX_N_DIMS); return UTIL_ERROR_BAD_INPUT; /*** ERROR RETURN ***/ } /******************************************************************************/ /* * interpolation order */ { CCTK_INT order; status = Util_TableGetInt(param_table_handle, &order, "order"); if (status == 1) then { } /* ok ==> 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\n" " for \"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 ***/ } /* * out-of-range interpolation-point handling */ { CCTK_REAL out_of_range_tolerance[2*MAX_N_DIMS]; const int N_tolerances = 2*N_dims; status = Util_TableGetRealArray(param_table_handle, N_tolerances, out_of_range_tolerance, "out_of_range_tolerance"); if (status == UTIL_ERROR_TABLE_NO_SUCH_KEY) then { /* default */ int tol_index; for (tol_index = 0 ; tol_index < N_tolerances ; ++tol_index) { out_of_range_tolerance[tol_index] = OUT_OF_RANGE_TOLERANCE_DEFAULT; } } else if (status == N_tolerances) then { /* check that all values are valid */ int tol_index; for (tol_index = 0 ; tol_index < N_tolerances ; ++tol_index) { if (! ( (out_of_range_tolerance[tol_index] >= 0.0) || (out_of_range_tolerance[tol_index] == -1.0) ) ) then { CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, __LINE__, __FILE__, CCTK_THORNSTRING, "\n" " CCTK_InterpLocalUniform():\n" " invalid out_of_range_tolerance[tol_index=%d] = %g!\n" " (valid values are -1.0 or >= 0.0)", tol_index, out_of_range_tolerance[tol_index]); return UTIL_ERROR_BAD_INPUT; /*** ERROR RETURN ***/ } } } else { CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, __LINE__, __FILE__, CCTK_THORNSTRING, "\n" " CCTK_InterpLocalUniform(): bad table entry\n" " for \"out_of_range_tolerance\" parameter!\n" " error status=%d", status); return status; /*** ERROR RETURN ***/ } /******************************************************************************/ /* * molecule family */ { #define MOLECULE_FAMILY_BUFSIZ (MAX_MOLECULE_FAMILY_STRLEN+1) char molecule_family_string[MOLECULE_FAMILY_BUFSIZ]; enum molecule_family molecule_family; status = Util_TableGetString(param_table_handle, MOLECULE_FAMILY_BUFSIZ, molecule_family_string, "molecule_family"); if (status == UTIL_ERROR_TABLE_NO_SUCH_KEY) then { /* default */ /* (need enum for main code, string for error messages) */ strcpy(molecule_family_string, "cube"); molecule_family = molecule_family_cube; } else if (status == 0) then { /* this is a query of our default molecule family */ status = Util_TableSetString(param_table_handle, "cube", "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, "cube") == 0) then molecule_family = molecule_family_cube; else { CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, __LINE__, __FILE__, CCTK_THORNSTRING, "InterpLocalUniform(): unknown molecule_family=\"%s\"!", molecule_family_string); return UTIL_ERROR_BAD_INPUT; /*** ERROR RETURN ***/ } } else { CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, __LINE__, __FILE__, CCTK_THORNSTRING, "\n" " CCTK_InterpLocalUniform(): bad table entry\n" " for \"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; /* default */ 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 \"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 = 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) 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\n" " for \"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) then { /* default stride = Fortran storage order, */ /* determined from input_array_dims[] */ if (input_array_dims == NULL) then { free(input_array_offsets); CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, __LINE__, __FILE__, CCTK_THORNSTRING, "\n" " CCTK_InterpLocalUniform(): can't default \"input_array_strides\"\n" " with input_array_dims == NULL!"); return UTIL_ERROR_BAD_INPUT; /*** ERROR RETURN ***/ } { int stride = 1; int axis; for (axis = 0 ; axis < N_dims ; ++axis) { input_array_strides[axis] = stride; stride *= input_array_dims[axis]; } } } 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\n" " for \"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) then { /* default min subscript = 0 along each axis */ LocalInterp_zero_int_array(N_dims, input_array_min_subscripts); } 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\n" " for \"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) then { /* default max subscript = input_array_dims[]-1 along each axis */ if (input_array_dims == NULL) then { free(input_array_offsets); CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, __LINE__, __FILE__, CCTK_THORNSTRING, "\n" " CCTK_InterpLocalUniform(): can't default \"input_array_max_subscripts\"\n" " with input_array_dims == NULL!"); return UTIL_ERROR_BAD_INPUT; /*** ERROR RETURN ***/ } { int axis; for (axis = 0 ; axis < N_dims ; ++axis) { input_array_max_subscripts[axis] = input_array_dims[axis] - 1; } } } 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\n" " for \"input_array_max_subscripts\" parameter!\n" " error status=%d", status); return status; /*** ERROR RETURN ***/ } /******************************************************************************/ /* * operand indices */ { CCTK_INT *const operand_indices = malloc(N_output_arrays1 * sizeof(CCTK_INT)); if (operand_indices == NULL) then { 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) 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); 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 ***/ } { int out; for (out = 0 ; out < N_output_arrays ; ++out) { operand_indices[out] = out; } } } 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\n" " for \"operand_indices\" parameter!\n" " error status=%d", status); return status; /*** ERROR RETURN ***/ } /* * operation_codes */ { CCTK_INT *const operation_codes = malloc(N_output_arrays1 * sizeof(CCTK_INT)); if (operation_codes == NULL) then { free(operand_indices); free(input_array_offsets); 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) 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\n" " for \"operation_codes\" parameter!\n" " error status=%d", status); return status; /*** ERROR RETURN ***/ } /******************************************************************************/ /* * set up for any molecule min/max m and position queries */ /* molecule min/max m */ status1 = Util_TableQueryValueInfo(param_table_handle, NULL, NULL, "molecule_min_m"); status2 = Util_TableQueryValueInfo(param_table_handle, NULL, NULL, "molecule_max_m"); if ((status1 < 0) || (status2 < 0)) then { free(operation_codes); free(operand_indices); free(input_array_offsets); CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, __LINE__, __FILE__, CCTK_THORNSTRING, "\n" " CCTK_InterpLocalUniform(): bad \"interp_molecule_{min,max}_m\"\n" " table entry/entries!\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; /* 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) 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); CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, __LINE__, __FILE__, CCTK_THORNSTRING, "\n" " CCTK_InterpLocalUniform(): bad \"molecule_positions\" table entry!\n" " error status=%d", status); return status; /*** ERROR RETURN ***/ } /******************************************************************************/ /* * set up for any Jacobian queries */ { 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 */ status = Util_TableQueryValueInfo(param_table_handle, NULL, NULL, "Jacobian_pointer"); if (status < 0) then { free(operation_codes); free(operand_indices); free(input_array_offsets); 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) 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 ***/ } } Jacobian_info.Jacobian_offset = (CCTK_INT *) malloc(N_output_arrays1 * sizeof(CCTK_INT)); if (Jacobian_info.Jacobian_offset == 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_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) 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); return status; /*** ERROR RETURN ***/ } 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)) 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 ***/ } status = Util_TableGetInt(param_table_handle, & Jacobian_info.Jacobian_part_stride, "Jacobian_part_stride"); if (status == 1) 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); 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" /* 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 */ } }; /* 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) 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():\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, out_of_range_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_X_RANGE) 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, "\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 ***/ } } /******************************************************************************/ /* * store molecule structure info */ 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"); 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" , status1, status2, status3, status4); return (status1 < 0) ? status1 /*** ERROR RETURN ***/ : (status2 < 0) ? status2 /*** ERROR RETURN ***/ : (status3 < 0) ? status3 /*** ERROR RETURN ***/ : status4; /*** ERROR RETURN ***/ } /******************************************************************************/ /* * store results from molecule min/max m queries * * ... molecule position queries are stored directly into * their final destination by subfunction * ==> no further store or cleanup needed here */ if (p_molecule_min_max_m_info != NULL) 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)) 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 ***/ } } /******************************************************************************/ 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); } return return_code; } } } } } } } } } } } } } } } }