diff options
Diffstat (limited to 'src/GeneralizedPolynomial-Uniform/InterpLocalUniform.c')
| -rw-r--r-- | src/GeneralizedPolynomial-Uniform/InterpLocalUniform.c | 293 |
1 files changed, 198 insertions, 95 deletions
diff --git a/src/GeneralizedPolynomial-Uniform/InterpLocalUniform.c b/src/GeneralizedPolynomial-Uniform/InterpLocalUniform.c index 56f386a..1f63a03 100644 --- a/src/GeneralizedPolynomial-Uniform/InterpLocalUniform.c +++ b/src/GeneralizedPolynomial-Uniform/InterpLocalUniform.c @@ -38,7 +38,7 @@ #include "Hermite/all_prototypes.h" /* the rcs ID and its dummy function to use it */ -static const char *rcsid = "$Header$"; +static const char* rcsid = "$Header$"; CCTK_FILEVERSION(CactusBase_LocalInterp_src_GeneralizedPolynomialUniform_InterpLocalUniform_c) /******************************************************************************/ @@ -231,8 +231,20 @@ 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 + @hauthor Jonathan Thornburg <jthorn@aei.mpg.de> + @hdesc Support parameter-table entries + @var N_boundary_points_to_omit, + @var boundary_off_centering_tolerance, and + @var boundary_extrapolation_tolerance. + @endhdesc + ***** misc arguments ***** @var interp_operator @@ -349,26 +361,43 @@ return LocalInterp_InterpLocalUniform(interp_operator_Hermite, @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 + @var N_boundary_points_to_omit + @vdesc If this key is present, then this array specifies how many + input grid points to omit (not use for the interpolation) + on each grid boundary. 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] + If this key isn't present, it defaults to all zeros, i.e. + to use all the input grid points in the interpolation. + @vtype const CCTK_INT N_boundary_points_to_omit[2*N_dims] + @endvar + + @var boundary_off_centering_tolerance + @vdesc If this key is present, then the interpolator allows + interpolation points to be up to (<=) this many grid + spacings outside the default-centering region on each + grid boundary, off-centering the interpolation molecules + as necessary. + 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, ...]. + If this key isn't present, it defaults to all infinities, + i.e. to place no restriction on the interpolation point. + @vtype const CCTK_REAL boundary_off_centering_tolerance[2*N_dims] + @endvar + + @var boundary_extrapolation_tolerance + @vdesc If this key is present, then the interpolator allows + interpolation points to be up to (<=) this many grid + spacings outside the valid-data region on each grid + boundary, doing extrapolation instead of interpolation + as necessary. + 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, ...]. + If this key isn't present, it defaults to all 1.0e-10, + i.e. to allow up to 1e-10 grid spacings of extrapolation. + @vtype const CCTK_REAL boundary_extrapolation_tolerance[2*N_dims] @endvar @var input_array_offsets @@ -625,7 +654,7 @@ static * 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)); + * CCTK_INT* const p = (CCTK_INT*) 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. :( :( @@ -657,7 +686,7 @@ if ( (N_dims <= 0) then { CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, __LINE__, __FILE__, CCTK_THORNSTRING, -"InterpLocalUniform(): invalid arguments\n" +"CCTK_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 ***/ @@ -684,14 +713,14 @@ if (N_dims > MAX_N_DIMS) CCTK_INT order; status = Util_TableGetInt(param_table_handle, &order, "order"); if (status == 1) - then { } /* ok ==> no-op here */ + 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\n" -" for \"order\" parameter!\n" +" CCTK_InterpLocalUniform(): bad or missing table entry for\n" +" \"order\" parameter!\n" " (this is a mandatory parameter)\n" " error status=%d", status); @@ -708,59 +737,122 @@ if ((order < 1) || (order > MAX_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"); +CCTK_INT N_boundary_points_to_omit[MAX_N_BOUNDARIES]; +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 = 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 { - /* default */ - int tol_index; - for (tol_index = 0 ; tol_index < N_tolerances ; ++tol_index) + /* set the default value */ + int ibndry; + for (ibndry = 0 ; ibndry < N_boundaries ; ++ibndry) { - out_of_range_tolerance[tol_index] - = OUT_OF_RANGE_TOLERANCE_DEFAULT; + N_boundary_points_to_omit[ibndry] = 0; } } -else if (status == N_tolerances) +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 ***/ + } + +/**************************************/ + +status = Util_TableGetRealArray(param_table_handle, + N_boundaries, boundary_off_centering_tolerance, + "boundary_off_centering_tolerance"); +if (status == UTIL_ERROR_TABLE_NO_SUCH_KEY) then { - /* check that all values are valid */ - int tol_index; - for (tol_index = 0 ; tol_index < N_tolerances ; ++tol_index) + /* set the default value */ + int ibndry; + for (ibndry = 0 ; ibndry < N_boundaries ; ++ibndry) { - 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, + 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():\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 ***/ - } +" CCTK_InterpLocalUniform(): bad table entry for\n" +" \"boundary_off_centering_tolerance\" parameter!\n" +" error status=%d", + status); + 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\n" -" for \"out_of_range_tolerance\" parameter!\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(): warning: setting\n" +" boundary_off_centering_tolerance[%d] == 0.0\n" +" and\n" +" boundary_extrapolation_tolerance[%d] > 0.0\n" +" is almost certainly a mistake\n" +" (the boundary_extrapolation_tolerance[%d]\n" +" setting will have no effect because of the\n" +" boundary_off_centering_tolerance[%d] setting)" + , + ibndry, ibndry, ibndry, ibndry); + } + } + /******************************************************************************/ /* @@ -768,23 +860,23 @@ else { */ { #define MOLECULE_FAMILY_BUFSIZ (MAX_MOLECULE_FAMILY_STRLEN+1) -char molecule_family_string[MOLECULE_FAMILY_BUFSIZ]; +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_string, + MOLECULE_FAMILY_BUFSIZ, molecule_family_strbuf, "molecule_family"); if (status == UTIL_ERROR_TABLE_NO_SUCH_KEY) then { - /* default */ + /* set the default value for code here*/ /* (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 */ + 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, - "cube", + molecule_family_str, "molecule_family"); if (status < 0) then { @@ -801,13 +893,13 @@ else if (status == 0) else if (status > 0) then { /* decode molecule family string */ - if (strcmp(molecule_family_string, "cube") == 0) + if (strcmp(molecule_family_strbuf, "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); +"CCTK_InterpLocalUniform(): unknown molecule_family=\"%s\"!", + molecule_family_strbuf); return UTIL_ERROR_BAD_INPUT; /*** ERROR RETURN ***/ } } @@ -815,13 +907,15 @@ else { CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, __LINE__, __FILE__, CCTK_THORNSTRING, "\n" -" CCTK_InterpLocalUniform(): bad table entry\n" -" for \"molecule_family\" parameter!\n" +" CCTK_InterpLocalUniform(): bad table entry for\n" +" \"molecule_family\" parameter!\n" " error status=%d", status); return status; /*** ERROR RETURN ***/ } +/**************************************/ + /* * smoothing */ @@ -829,14 +923,15 @@ else { CCTK_INT smoothing; status = Util_TableGetInt(param_table_handle, &smoothing, "smoothing"); if (status == UTIL_ERROR_TABLE_NO_SUCH_KEY) - then smoothing = 0; /* default */ + 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 \"smoothing\" parameter!\n" +" CCTK_InterpLocalUniform(): bad table entry for\n" +" \"smoothing\" parameter!\n" " error status=%d", status); return status; /*** ERROR RETURN ***/ @@ -850,8 +945,8 @@ if ((smoothing < 0) || (smoothing > MAX_SMOOTHING)) * input array offsets */ { -CCTK_INT *const input_array_offsets - = malloc(N_input_arrays1 * sizeof(CCTK_INT)); +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, @@ -869,13 +964,15 @@ else { CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, __LINE__, __FILE__, CCTK_THORNSTRING, "\n" -" CCTK_InterpLocalUniform(): bad table entry\n" -" for \"input_array_offsets\" parameter!\n" +" CCTK_InterpLocalUniform(): bad table entry for\n" +" \"input_array_offsets\" parameter!\n" " error status=%d", status); return status; /*** ERROR RETURN ***/ } +/**************************************/ + /* * input array strides */ @@ -915,13 +1012,15 @@ else { CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, __LINE__, __FILE__, CCTK_THORNSTRING, "\n" -" CCTK_InterpLocalUniform(): bad table entry\n" -" for \"input_array_offsets\" parameter!\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 */ @@ -942,8 +1041,8 @@ else { CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, __LINE__, __FILE__, CCTK_THORNSTRING, "\n" -" CCTK_InterpLocalUniform(): bad table entry\n" -" for \"input_array_min_subscripts\" parameter!\n" +" CCTK_InterpLocalUniform(): bad table entry for\n" +" \"input_array_min_subscripts\" parameter!\n" " error status=%d", status); return status; /*** ERROR RETURN ***/ @@ -981,8 +1080,8 @@ else { CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, __LINE__, __FILE__, CCTK_THORNSTRING, "\n" -" CCTK_InterpLocalUniform(): bad table entry\n" -" for \"input_array_max_subscripts\" parameter!\n" +" CCTK_InterpLocalUniform(): bad table entry for\n" +" \"input_array_max_subscripts\" parameter!\n" " error status=%d", status); return status; /*** ERROR RETURN ***/ @@ -994,8 +1093,8 @@ else { * operand indices */ { -CCTK_INT *const operand_indices - = malloc(N_output_arrays1 * sizeof(CCTK_INT)); +CCTK_INT* const operand_indices + = (CCTK_INT*) malloc(N_output_arrays1 * sizeof(CCTK_INT)); if (operand_indices == NULL) then { free(input_array_offsets); @@ -1059,19 +1158,21 @@ else { CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, __LINE__, __FILE__, CCTK_THORNSTRING, "\n" -" CCTK_InterpLocalUniform(): bad table entry\n" -" for \"operand_indices\" parameter!\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 - = malloc(N_output_arrays1 * sizeof(CCTK_INT)); +CCTK_INT* const operation_codes + = (CCTK_INT*) malloc(N_output_arrays1 * sizeof(CCTK_INT)); if (operation_codes == NULL) then { free(operand_indices); @@ -1110,8 +1211,8 @@ else { CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, __LINE__, __FILE__, CCTK_THORNSTRING, "\n" -" CCTK_InterpLocalUniform(): bad table entry\n" -" for \"operation_codes\" parameter!\n" +" CCTK_InterpLocalUniform(): bad table entry for\n" +" \"operation_codes\" parameter!\n" " error status=%d", status); return status; /*** ERROR RETURN ***/ @@ -1223,7 +1324,7 @@ if (status) then { /* yes, we're doing Jacobian queries */ Jacobian_info.Jacobian_pointer - = (CCTK_REAL **) malloc(N_output_arrays1 * sizeof(CCTK_REAL *)); + = (CCTK_REAL**) malloc(N_output_arrays1 * sizeof(CCTK_REAL *)); if (Jacobian_info.Jacobian_pointer == NULL) then { free(operation_codes); @@ -1232,8 +1333,8 @@ if (status) return UTIL_ERROR_NO_MEMORY; /*** ERROR RETURN ***/ } { - CCTK_POINTER *Jacobian_pointer_temp - = (CCTK_POINTER *) malloc(N_output_arrays1 * sizeof(CCTK_POINTER)); + CCTK_POINTER* Jacobian_pointer_temp + = (CCTK_POINTER*) malloc(N_output_arrays1 * sizeof(CCTK_POINTER)); if (Jacobian_pointer_temp == NULL) then { free(operation_codes); @@ -1276,7 +1377,7 @@ if (status) } Jacobian_info.Jacobian_offset - = (CCTK_INT *) malloc(N_output_arrays1 * sizeof(CCTK_INT)); + = (CCTK_INT*) malloc(N_output_arrays1 * sizeof(CCTK_INT)); if (Jacobian_info.Jacobian_offset == NULL) then { free(operation_codes); @@ -1550,7 +1651,7 @@ if (p_interp_fn == NULL_INTERP_FN_PTR) " interp_operator=\"%s\", N_dims=%d\n" " molecule_family=\"%s\", order=%d, smoothing=%d", interp_operator_string, N_dims, - molecule_family_string, order, smoothing); + molecule_family_str, order, smoothing); return UTIL_ERROR_BAD_INPUT; /*** ERROR RETURN ***/ } @@ -1564,7 +1665,9 @@ 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_boundary_points_to_omit, + boundary_off_centering_tolerance, + boundary_extrapolation_tolerance, N_input_arrays, input_array_offsets, input_array_strides, input_array_min_subscripts, |