From 1e1d195e8148ec50b9ff287139322ad135d84d0b Mon Sep 17 00:00:00 2001 From: jthorn Date: Wed, 29 Jan 2003 13:34:53 +0000 Subject: Code changes * support the same off-centering and boundary options as the new global interpolator (which needs the support here to work): this means changing the former parameter-table entry off_centering_tolerance to the new entries boundary_off_centering_tolerance boundary_extrapolation_tolerance * implement the N_boundary_points_to_omit parameter-table entry, (previously this was documented but not implemented) * completely rewrite the LocalInterp_molecule_posn() function and its test driver to support the above changes * change error return code from old CCTK_ERROR_INTERP_POINT_X_RANGE to new CCTK_ERROR_INTERP_POINT_OUTSIDE to avoid confusion with the old name suggesting that this was specific to the X coordinate * fix documentation of interpolation operator names for older (uniform Cartesian) CCTK_InterpLocal() interpolator -- this closes Erik Schnetter's bug report CactusBase/1367 Documentation changes * document the above code changes * much expanded description of molecule centering and boundary handling, including 2 new diagrams (done in latex picture mode :) * update README files to list all code files * move caching parameters (not implemented, and probably won't be implemented any time soon) out of documentation.tex into a new file future-ideas.tex * general cleanup of documentation.tex (not quite done, but 80% or so...) to more clearly describe what this thorn actually implements, rather than a generic spec of which this thorn only implements some pieces git-svn-id: http://svn.cactuscode.org/arrangements/CactusBase/LocalInterp/trunk@125 df1f8a13-aa1d-4dd4-9681-27ded5b42416 --- .../InterpLocalUniform.c | 293 +++++++++----- .../InterpLocalUniform.h | 39 +- src/GeneralizedPolynomial-Uniform/README | 29 +- src/GeneralizedPolynomial-Uniform/molecule_posn.c | 177 ++++---- src/GeneralizedPolynomial-Uniform/template.c | 200 ++++----- src/GeneralizedPolynomial-Uniform/template.h | 8 +- .../test_molecule_posn.c | 447 +++++++++++---------- 7 files changed, 654 insertions(+), 539 deletions(-) (limited to 'src') 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 + @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, diff --git a/src/GeneralizedPolynomial-Uniform/InterpLocalUniform.h b/src/GeneralizedPolynomial-Uniform/InterpLocalUniform.h index d2990df..851378a 100644 --- a/src/GeneralizedPolynomial-Uniform/InterpLocalUniform.h +++ b/src/GeneralizedPolynomial-Uniform/InterpLocalUniform.h @@ -70,11 +70,17 @@ typedef int bool; */ #define MAX_N_DIMS 3 +/* a "boundary" is the combination of a dimension and a min/max "side" */ +#define MAX_N_BOUNDARIES (2*MAX_N_DIMS) + /* * if molecule_family_string is a C string specifying a molecule family * (i.e. value associated with the "molecule_family" key in the parameter * table), we must have * strlen(molecule_family_string) <= MAX_MOLECULE_FAMILY_STRLEN + * n.b. exceeding this won't cause a buffer overflow, it will "just" + * cause the string to be truncated (and probably not recognized + * by the interpolator) */ #define MAX_MOLECULE_FAMILY_STRLEN 20 @@ -103,12 +109,15 @@ enum molecule_family * other compile-time settings */ -/* default for each element of out_of_range_tolerance[] */ -/* FIXME: this is based on C "double", which may not match CCTK_REAL */ -#define OUT_OF_RANGE_TOLERANCE_DEFAULT (10000.0*DBL_EPSILON) +/* 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 -/* CCTK_VWarn() severity level for error messages */ -#define ERROR_MSG_SEVERITY_LEVEL 1 +/* CCTK_VWarn() severity level for error/warning messages */ +#define ERROR_MSG_SEVERITY_LEVEL 0 +#define WARNING_MSG_SEVERITY_LEVEL 1 /******************************************************************************/ @@ -122,6 +131,7 @@ enum molecule_family struct error_flags { int error_pt; + int error_ibndry; int error_axis; int error_end; }; @@ -153,20 +163,16 @@ struct Jacobian_info /* * error codes for LocalInterp_molecule_posn() - * ... these are defined in terms of INT_MIN from */ /* x < minimum allowable x in grid */ -#define MOLECULE_POSN_ERROR_X_LT_MIN (INT_MIN + 0) +#define MOLECULE_POSN_ERROR_X_LT_MIN (-1) /* x > maximum allowable x in grid */ -#define MOLECULE_POSN_ERROR_X_GT_MAX (INT_MIN + 1) +#define MOLECULE_POSN_ERROR_X_GT_MAX (-2) /* grid is smaller than molecule */ -#define MOLECULE_POSN_ERROR_GRID_TINY (INT_MIN + 2) - -/* is a given integer an error code? */ -#define IS_MOLECULE_POSN_ERROR_CODE(x) (x <= MOLECULE_POSN_ERROR_GRID_TINY) +#define MOLECULE_POSN_ERROR_GRID_TINY (-3) /******************************************************************************/ @@ -216,11 +222,12 @@ int LocalInterp_ILU_Hermite(int N_dims, int LocalInterp_molecule_posn(fp grid_origin, fp grid_delta, int grid_i_min, int grid_i_max, int molecule_size, - fp out_of_range_tolerance_min, - fp out_of_range_tolerance_max, + fp boundary_off_centering_tolerance_min, + fp boundary_off_centering_tolerance_max, + fp boundary_extrapolation_tolerance_min, + fp boundary_extrapolation_tolerance_max, fp x, - fp *x_rel, - int *molecule_m_min, int *molecule_m_max); + int* i_center, fp* x_rel); /* functions in util.c */ void LocalInterp_zero_int_array(int N, CCTK_INT array[]); diff --git a/src/GeneralizedPolynomial-Uniform/README b/src/GeneralizedPolynomial-Uniform/README index 0963854..b022489 100644 --- a/src/GeneralizedPolynomial-Uniform/README +++ b/src/GeneralizedPolynomial-Uniform/README @@ -7,7 +7,8 @@ under the name The source code files are as follows: -* startup.c registers the interpolation operator +* startup.c registers the interpolation operators +* InterpLocalUniform.h is an overall header file for the whole interpolator * InterpLocalUniform.c is the top-level driver: it gets the various options from the parameter table, then decodes (N_dims, molecule_family, order, smoothing) @@ -16,6 +17,7 @@ The source code files are as follows: * [123]d.cube.order*.smooth*.c define the individual interpolation subfunctions. Each of them just #defines a whole bunch of macros, then #includes template.c (which has the actual code). +* template.h defines a prototype for the function defined by template.c * template.c is the actual interpolation code. It is written in terms of a large number of macros, which should be #defined before #including template.c. There's a long block comment @@ -26,6 +28,10 @@ The source code files are as follows: Maple-generated code fragments in the [123]d.coeffs/directories; template.c uses various macros to tell it which fragments to #include. +* molecule_posn.c contains the LocalInterp_molecule_posn() function + to compute where in the grid each (an) interpolation molecule should + be centered for each (a given) interpolation point. +* util.c contains some low-level utility routines * [123]d.maple are the top-level Maple code files; they call various functions in interpolate.maple and util.maple to do the actual work. @@ -33,7 +39,19 @@ The source code files are as follows: computing an interpolant and manipulating/printing it in various ways * util.maple contains low-level utility routines +* makefile is a makefile for... +* test_molecule_posn.c is a standalone test driver for the code + in molecule_posn.c . By default it runs a large set of (around 100) + test cases stored in a table in the code. +The subdirectories are as follows: +* [123]d.coeffs/ are empty directory trees left-over from previous + versions of this interpolator (CVS doesn't have any easy way to + remove directories) +* Lagrange/ contains the code for the Lagrange polynomial interpolator. +* Hermite/ contains the code for the Hermite polynomial interpolator. +* common/ contains low-level code common to both the Lagrange and + the Hermite interpolators. To add a new combination of (N_dims, molecule_family, order, smoothing), @@ -60,17 +78,10 @@ to this interpolator, you need to -Other makefile targets: -test_molecule_posn - This makes a standalone test driver for the - LocalInterp_molecule_posn() - function defined in molecule_posn.c - - Copyright ========= -This interpolator is copyright (C) 2002-2003 +This interpolator is copyright (C) 2001-2003 by Jonathan Thornburg . This interpolator is free software; you can redistribute it and/or modify diff --git a/src/GeneralizedPolynomial-Uniform/molecule_posn.c b/src/GeneralizedPolynomial-Uniform/molecule_posn.c index e007396..6d42ef3 100644 --- a/src/GeneralizedPolynomial-Uniform/molecule_posn.c +++ b/src/GeneralizedPolynomial-Uniform/molecule_posn.c @@ -8,7 +8,6 @@ @@*/ #include -#include #include #include @@ -73,6 +72,18 @@ static const char *rcsid = "$Header$"; center, i.e. the above diagram has columns labelled with [i+m]. @enddesc + @hdate 27.Jan.2003 + @hauthor Jonathan Thornburg + @hdesc Complete rewrite: now supports + @var boundary_off_centering_tolerance_min, + @var boundary_off_centering_tolerance_max, + @var boundary_extrapolation_tolerance_min, + @var boundary_extrapolation_tolerance_max, + also change to returning status code, + also drop returning @var min_m and @var max_m + since they were never used. + @endhdesc + @var grid_origin @vdesc The floating-point coordinate x of the grid point i=0. @vtype fp grid_origin @@ -98,50 +109,53 @@ static const char *rcsid = "$Header$"; @vtype int molecule_size @endvar - @var out_of_range_tolerance_min, out_of_range_tolerance_max - @vdesc Specifies how out-of-range interpolation points should - be handled for the {minimum,maximum} ends of the grid - respectively. - If out_of_range_tolerance >= 0.0, - then an interpolation point is considered to be - "out of range" if and only if the interpolation - point is > out_of_range_tolerance * grid_delta - outside the grid in any coordinate. - If out_of_range_tolerance == -1.0, - then an interpolation point is considered to be - "out of range" if and only if a centered molecule - would require data from outside the grid. - Other values of out_of_range_tolerance are illegal. - @vtype fp out_of_range_tolerance_min, out_of_range_tolerance_max + @var boundary_off_centering_tolerance_min, + boundary_off_centering_tolerance_max + @vdesc Specifies how many grid spacings the interpolation point + may be beyond the default-centering region before being + declared "out of range" on the {minimum,maximum} boundaries + of the grid respectively. + @vtype fp boundary_off_centering_tolerance_min, + boundary_off_centering_tolerance_max + @endvar + + @var boundary_extrapolation_tolerance_min, + boundary_extrapolation_tolerance_max + @vdesc Specifies how many grid spacings the interpolation point + may be beyond the grid boundary before being declared + "out of range" on the {minimum,maximum} boundaries + of the grid respectively. + @vtype fp boundary_extrapolation_tolerance_min, + boundary_extrapolation_tolerance_max @endvar @var x - @vdesc The floating-point coordinate x at which we would like to center - the molecule. + @vdesc The floating-point coordinate of the interpolation point. @vtype fp x @endvar - @var x_rel - @vdesc A pointer to a floating-point value where this function should - store the x coordinate relative to the molecule center, measured - in units of the grid spacing; or NULL to skip storing this. - @vtype fp *x_rel + @var i_center + @vdesc A pointer to an value where this function should + store the integer coordinate of the molecule center, + or NULL to skip storing this + @vtype int *i_center @vio pointer to out @endvar - @var molecule_m_min, molecule_m_max - @vdesc A pointer to an int where this function should store the - {minimum,maximum} molecule coordinate m of the molecule; + @var x_rel + @vdesc A pointer to where this function should store the + interpolation point's x coordinate relative to the + molecule center, measured in units of the grid spacing; or NULL to skip storing this. - @vtype int *molecule_m_min, *molecule_m_max + @vtype fp *x_rel @vio pointer to out @endvar @returntype int @returndesc - This function returns the integer coordinate i_center at which - the molecule should be centered, or one of the error codes defined - in "InterpLocalUniform.h": + This function returns 0 if the interpolation point is "in range", + or one of the (negative) error codes defined in "InterpLocalUniform.h" + if the interpolation point is "out of range": MOLECULE_POSN_ERROR_X_LT_MIN if x is out-of-range on the min end of the grid (i.e. x < the minimum allowable x) @@ -155,73 +169,70 @@ static const char *rcsid = "$Header$"; int LocalInterp_molecule_posn(fp grid_origin, fp grid_delta, int grid_i_min, int grid_i_max, int molecule_size, - fp out_of_range_tolerance_min, - fp out_of_range_tolerance_max, + fp boundary_off_centering_tolerance_min, + fp boundary_off_centering_tolerance_max, + fp boundary_extrapolation_tolerance_min, + fp boundary_extrapolation_tolerance_max, fp x, - fp *x_rel, - int *molecule_m_min, int *molecule_m_max) + int* i_center, fp* x_rel) { -/* min/max molecule coordinates */ -const int m_max = (molecule_size >> 1); /* a positive number */ -const int m_min = m_max - molecule_size + 1; /* a negative number */ +/* molecule radia (inherently positive) in +/- directions */ +const int mr_plus = (molecule_size >> 1); +const int mr_minus = molecule_size - mr_plus - 1; -/* integer coordinate i of point x, as a floating-point number */ -const fp fp_i = (x - grid_origin) / grid_delta; +/* range of x_rel for which this molecule is centered, cf. diagram in header comment */ +const fp centered_min_x_rel = IS_EVEN(molecule_size) ? 0.0 : -0.5; +const fp centered_max_x_rel = IS_EVEN(molecule_size) ? 1.0 : +0.5; -/* is point x out-of-range? */ -if (out_of_range_tolerance_min >= 0.0) - then { - const fp fp_effective_grid_i_min - = ((fp) grid_i_min) - out_of_range_tolerance_min; - if (fp_i < fp_effective_grid_i_min) - then return INT_MIN; /*** ERROR RETURN ***/ - } -if (out_of_range_tolerance_max >= 0.0) - then { - const fp fp_effective_grid_i_max - = ((fp) grid_i_max) + out_of_range_tolerance_max; - if (fp_i > fp_effective_grid_i_max) - then return INT_MAX; /*** ERROR RETURN ***/ - } +/* range of i where we *could* center the molecule, as floating-point numbers */ +const fp fp_centered_min_possible_i = grid_i_min + mr_minus + centered_min_x_rel; +const fp fp_centered_max_possible_i = grid_i_max - mr_plus + centered_max_x_rel; -/* integer coordinate i where we will center the molecule */ -/* (see diagram in header comment for explanation) */ -/* ... as a floating-point number */ - { -const fp fp_i_center = IS_EVEN(molecule_size) ? floor(fp_i) : JT_ROUND(fp_i); -/* ... as an integer */ -int i_center = (int) fp_i_center; +/* integer coordinate i of interpolation point, as a floating-point number */ +const fp fp_i = (x - grid_origin) / grid_delta; -/* min/max integer coordinates in molecule assuming it's fully centered */ -const int centered_i_min = i_center + m_min; -const int centered_i_max = i_center + m_max; +/* is the molecule larger than the grid? */ +if (molecule_size > HOW_MANY_IN_RANGE(grid_i_min,grid_i_max)) + then return MOLECULE_POSN_ERROR_GRID_TINY; /*** ERROR RETURN ***/ -/* check if off-centered molecules are forbidden? */ -if ((out_of_range_tolerance_min == -1.0) && (centered_i_min < grid_i_min)) +/* is interpolation point x beyond the extrapolation tolerance? */ +if (fp_i < (fp)grid_i_min - boundary_extrapolation_tolerance_min) then return MOLECULE_POSN_ERROR_X_LT_MIN; /*** ERROR RETURN ***/ -if ((out_of_range_tolerance_max == -1.0) && (centered_i_max > grid_i_max)) +if (fp_i > (fp)grid_i_max + boundary_extrapolation_tolerance_max) then return MOLECULE_POSN_ERROR_X_GT_MAX; /*** ERROR RETURN ***/ -/* off-center as needed if we're close to the edge of the grid */ +/* is interpolation point x beyond the off-centering tolerance? */ +if (fp_i < fp_centered_min_possible_i - boundary_off_centering_tolerance_min) + then return MOLECULE_POSN_ERROR_X_LT_MIN; /*** ERROR RETURN ***/ +if (fp_i > fp_centered_max_possible_i + boundary_off_centering_tolerance_max) + then return MOLECULE_POSN_ERROR_X_GT_MAX; /*** ERROR RETURN ***/ + +/* now choose the actual molecule position/centering: */ +/* first set up a centered molecule */ { -int shift = 0; -if (molecule_size > HOW_MANY_IN_RANGE(grid_i_min,grid_i_max)) - then return MOLECULE_POSN_ERROR_GRID_TINY; /*** ERROR RETURN ***/ -if (centered_i_min < grid_i_min) - then shift = grid_i_min - centered_i_min; /* a positive number */ -if (centered_i_max > grid_i_max) - then shift = grid_i_max - centered_i_max; /* a negative number */ -i_center += shift; +fp fp_i_center = IS_EVEN(molecule_size) /* ... as a floating-point number */ + ? floor(fp_i) + : JT_ROUND(fp_i); +int int_i_center = (int) fp_i_center; /* ... as an integer */ + +/* then clamp the molecule at the grid boundaries */ +if (int_i_center - mr_minus < grid_i_min) + then { + int_i_center = grid_i_min + mr_minus; + fp_i_center = (fp) int_i_center; + } +if (int_i_center + mr_plus > grid_i_max) + then { + int_i_center = grid_i_max - mr_plus; + fp_i_center = (fp) int_i_center; + } /* store the results */ +if (i_center != NULL) + then *i_center = int_i_center; if (x_rel != NULL) - then *x_rel = fp_i - ((fp) i_center); -if (molecule_m_min != NULL) - then *molecule_m_min = m_min; -if (molecule_m_max != NULL) - then *molecule_m_max = m_max; + then *x_rel = fp_i - fp_i_center; -return i_center; - } +return 0; /*** NORMAL RETURN ***/ } } diff --git a/src/GeneralizedPolynomial-Uniform/template.c b/src/GeneralizedPolynomial-Uniform/template.c index 0bb409b..a567119 100644 --- a/src/GeneralizedPolynomial-Uniform/template.c +++ b/src/GeneralizedPolynomial-Uniform/template.c @@ -20,6 +20,13 @@ 3-D molecules. @endhdesc + @hdate 28.Jan.2003 + @hauthor Jonathan Thornburg + @hdesc Support parameter-table entries + @var boundary_off_centering_tolerance and + @var boundary_extrapolation_tolerance. + @endhdesc + @desc The following header files must be #included before #including this file: @@ -301,6 +308,9 @@ only briefly; see the InterpLocalUniform() documentation and/or the thorn guide for this thorn for further details. + If you change the arguments for this function, note that you + must also change the prototype in "template.h". + @var error_flags @vdesc If we encounter an out-of-range point and this pointer is non-NULL, we store a description of the out-of-range @@ -346,13 +356,15 @@ @@*/ int FUNCTION_NAME(/***** coordinate system *****/ - const CCTK_REAL coord_system_origin[], - const CCTK_REAL grid_spacing[], + 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[], - const CCTK_REAL out_of_range_tolerance[], + const CCTK_INT N_boundary_points_to_omit[], + const CCTK_REAL boundary_off_centering_tolerance[], + const CCTK_REAL boundary_extrapolation_tolerance[], /***** input arrays *****/ int N_input_arrays, const CCTK_INT input_array_offsets[], @@ -691,8 +703,8 @@ int out; */ { #if N_DIMS >= 1 - const fp origin_x = coord_system_origin[X_AXIS]; - const fp delta_x = grid_spacing[X_AXIS]; + const fp origin_x = coord_origin[X_AXIS]; + const fp delta_x = coord_delta[X_AXIS]; #if defined(HAVE_OP_DX) \ || defined(HAVE_OP_DXY) || defined(HAVE_OP_DXZ) \ || defined(HAVE_OP_DXX) @@ -700,8 +712,8 @@ int out; #endif #endif #if N_DIMS >= 2 - const fp origin_y = coord_system_origin[Y_AXIS]; - const fp delta_y = grid_spacing[Y_AXIS]; + const fp origin_y = coord_origin[Y_AXIS]; + const fp delta_y = coord_delta[Y_AXIS]; #if defined(HAVE_OP_DY) \ || defined(HAVE_OP_DXY) || defined(HAVE_OP_DYZ) \ || defined(HAVE_OP_DYY) @@ -709,8 +721,8 @@ int out; #endif #endif #if N_DIMS >= 3 - const fp origin_z = coord_system_origin[Z_AXIS]; - const fp delta_z = grid_spacing[Z_AXIS]; + const fp origin_z = coord_origin[Z_AXIS]; + const fp delta_z = coord_delta[Z_AXIS]; #if defined(HAVE_OP_DZ) \ || defined(HAVE_OP_DXZ) || defined(HAVE_OP_DYZ) \ || defined(HAVE_OP_DZZ) @@ -764,6 +776,7 @@ int out; /* * interpolate at each point */ +int return_code = 0; int pt; for (pt = 0 ; pt < N_interp_points ; ++pt) { @@ -890,124 +903,84 @@ int pt; * * n.b. we need the final answers in variables with the magic * names (x,y,z) (machine-generated code uses these names), - * but we use temp variables as intermediates for (likely) - * better performance: the temp variables have their addresses - * taken and so may not be register-allocated, whereas the - * final (x,y,z) are "clean" and thus more likely to be - * register-allocated + * but we use temp variables as intermediates for these and + * for center_[ijk] for (likely) better performance: + * the temp variables have their addresses taken and so + * may not be register-allocated, whereas the final variables + * are "clean" and thus more likely to be register-allocated */ { - #if (N_DIMS >= 1) - fp x_temp; - const int center_i - = LocalInterp_molecule_posn(origin_x, delta_x, - input_array_min_subscripts[X_AXIS], - input_array_max_subscripts[X_AXIS], - MOLECULE_SIZE, - out_of_range_tolerance[X_AXIS_MIN], - out_of_range_tolerance[X_AXIS_MAX], - interp_coords_fp[X_AXIS], - &x_temp, - (int *) NULL, (int *) NULL); - const fp x = x_temp; - #endif - #if (N_DIMS >= 2) - fp y_temp; - const int center_j - = LocalInterp_molecule_posn(origin_y, delta_y, - input_array_min_subscripts[Y_AXIS], - input_array_max_subscripts[Y_AXIS], - MOLECULE_SIZE, - out_of_range_tolerance[Y_AXIS_MIN], - out_of_range_tolerance[Y_AXIS_MAX], - interp_coords_fp[Y_AXIS], - &y_temp, - (int *) NULL, (int *) NULL); - const fp y = y_temp; - #endif - #if (N_DIMS >= 3) - fp z_temp; - const int center_k - = LocalInterp_molecule_posn(origin_z, delta_z, - input_array_min_subscripts[Z_AXIS], - input_array_max_subscripts[Z_AXIS], - MOLECULE_SIZE, - out_of_range_tolerance[Z_AXIS_MIN], - out_of_range_tolerance[Z_AXIS_MAX], - interp_coords_fp[Z_AXIS], - &z_temp, - (int *) NULL, (int *) NULL); - const fp z = z_temp; - #endif - #if (N_DIMS > 3) - #error "N_DIMS may not be > 3!" - #endif - - /* is the interpolation point out-of-range? */ - #if (N_DIMS >= 1) - if (IS_MOLECULE_POSN_ERROR_CODE(center_i)) - then { - if (center_i == MOLECULE_POSN_ERROR_GRID_TINY) - then return CCTK_ERROR_INTERP_GRID_TOO_TINY; - if (error_flags != NULL) - then { - error_flags->error_pt = pt; - error_flags->error_axis = X_AXIS; - error_flags->error_end - = (center_i == MOLECULE_POSN_ERROR_X_GT_MAX) ? +1 : -1; - } - return CCTK_ERROR_INTERP_POINT_X_RANGE; /*** ERROR RETURN ***/ - } - #endif - #if (N_DIMS >= 2) - if (IS_MOLECULE_POSN_ERROR_CODE(center_j)) - then { - if (center_j == MOLECULE_POSN_ERROR_GRID_TINY) - then return CCTK_ERROR_INTERP_GRID_TOO_TINY; - if (error_flags != NULL) - then { - error_flags->error_pt = pt; - error_flags->error_axis = Y_AXIS; - error_flags->error_end - = (center_j == MOLECULE_POSN_ERROR_X_GT_MAX) ? +1 : -1; - } - return CCTK_ERROR_INTERP_POINT_X_RANGE; /*** ERROR RETURN ***/ - } - #endif - #if (N_DIMS >= 3) - if (IS_MOLECULE_POSN_ERROR_CODE(center_k)) - then { - if (center_k == MOLECULE_POSN_ERROR_GRID_TINY) - then return CCTK_ERROR_INTERP_GRID_TOO_TINY; - if (error_flags != NULL) - then { - error_flags->error_pt = pt; - error_flags->error_axis = Z_AXIS; - error_flags->error_end - = (center_k == MOLECULE_POSN_ERROR_X_GT_MAX) ? +1 : -1; - } - return CCTK_ERROR_INTERP_POINT_X_RANGE; /*** ERROR RETURN ***/ - } - #endif + int center_ijk_temp[MAX_N_DIMS]; + fp xyz_temp[MAX_N_DIMS]; + int axis; + for (axis = 0 ; axis < N_DIMS ; ++axis) + { + const int ibndry_min = 2*axis; + const int ibndry_max = 2*axis + 1; + const int status = LocalInterp_molecule_posn + (coord_origin[axis], coord_delta[axis], + input_array_min_subscripts[axis] + + N_boundary_points_to_omit[ibndry_min], + input_array_max_subscripts[axis] + - N_boundary_points_to_omit[ibndry_max], + MOLECULE_SIZE, + boundary_off_centering_tolerance[ibndry_min], + boundary_off_centering_tolerance[ibndry_max], + boundary_extrapolation_tolerance[ibndry_min], + boundary_extrapolation_tolerance[ibndry_max], + interp_coords_fp[axis], + ¢er_ijk_temp[axis], &xyz_temp[axis]); + if (status < 0) + then { + if (status == MOLECULE_POSN_ERROR_GRID_TINY) + then return CCTK_ERROR_INTERP_GRID_TOO_TINY; + /*** ERROR RETURN ***/ + if (error_flags != NULL) + then { + error_flags->error_pt = pt; + error_flags->error_axis = axis; + if (status == MOLECULE_POSN_ERROR_X_LT_MIN) + then { + error_flags->error_ibndry = ibndry_min; + error_flags->error_end = -1; + } + else { + error_flags->error_ibndry = ibndry_max; + error_flags->error_end = +1; + } + } + return CCTK_ERROR_INTERP_POINT_OUTSIDE; /*** ERROR RETURN ***/ + } + } + { + #if (N_DIMS >= 1) + const int center_i = center_ijk_temp[X_AXIS]; + const fp x = xyz_temp [X_AXIS]; + #endif + #if (N_DIMS >= 2) + const int center_j = center_ijk_temp[Y_AXIS]; + const fp y = xyz_temp [Y_AXIS]; + #endif + #if (N_DIMS >= 3) + const int center_k = center_ijk_temp[Z_AXIS]; + const fp z = xyz_temp [Z_AXIS]; + #endif /* * compute 1-d position of molecule center in input arrays */ { - #if (N_DIMS == 1) + #if (N_DIMS == 1) const int molecule_center_posn = stride_i*center_i; - #endif - #if (N_DIMS == 2) + #elif (N_DIMS == 2) const int molecule_center_posn = stride_i*center_i + stride_j*center_j; - #endif - #if (N_DIMS == 3) + #elif (N_DIMS == 3) const int molecule_center_posn = stride_i*center_i + stride_j*center_j + stride_k*center_k; - #endif - #if (N_DIMS > 3) + #else #error "N_DIMS may not be > 3!" #endif @@ -1566,6 +1539,7 @@ CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, } } } + } /* end of for (pt = ...) loop */ } diff --git a/src/GeneralizedPolynomial-Uniform/template.h b/src/GeneralizedPolynomial-Uniform/template.h index 3bc2a13..b5eb5b4 100644 --- a/src/GeneralizedPolynomial-Uniform/template.h +++ b/src/GeneralizedPolynomial-Uniform/template.h @@ -8,13 +8,15 @@ @version $Header$ @@*/ int FUNCTION_NAME(/***** coordinate system *****/ - const CCTK_REAL coord_system_origin[], - const CCTK_REAL grid_spacing[], + 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[], - const CCTK_REAL out_of_range_tolerance[], + const CCTK_INT N_boundary_points_to_omit[], + const CCTK_REAL boundary_off_centering_tolerance[], + const CCTK_REAL boundary_extrapolation_tolerance[], /***** input arrays *****/ int N_input_arrays, const CCTK_INT input_array_offsets[], diff --git a/src/GeneralizedPolynomial-Uniform/test_molecule_posn.c b/src/GeneralizedPolynomial-Uniform/test_molecule_posn.c index bd3a1f9..4f6a5f5 100644 --- a/src/GeneralizedPolynomial-Uniform/test_molecule_posn.c +++ b/src/GeneralizedPolynomial-Uniform/test_molecule_posn.c @@ -7,11 +7,27 @@ * Usage: * test_molecule_posn # run a preset set of tests * test_molecule_posn molecule_size \ - * out_of_range_tolerance_min \ - * out_of_range_tolerance_max \ + * boundary_off_centering_tolerance_min \ + * boundary_off_centering_tolerance_max \ + * boundary_extrapolation_tolerance_min \ + * boundary_extrapolation_tolerance_max \ + * boundary_off_centering_tolerance_max \ * x # do a single test as specified */ +const char* help_msg = +"usage:\n" +"# run a preset series of tests:\n" +" test_molecule_posn\n" +"# run a single test as specified:\n" +" test_molecule_posn molecule_size \\\n" +" boundary_off_centering_tolerance_min \\\n" +" boundary_off_centering_tolerance_max \\\n" +" boundary_extrapolation_tolerance_min \\\n" +" boundary_extrapolation_tolerance_max \\\n" +" x\n" +; + #include #include #include @@ -39,8 +55,10 @@ const int grid_i_max = 105; /* x_max = 13.6 */ */ static void run_interactive_test(int molecule_size, - fp out_of_range_tolerance_min, - fp out_of_range_tolerance_max, + fp boundary_off_centering_tolerance_min, + fp boundary_off_centering_tolerance_max, + fp boundary_extrapolation_tolerance_min, + fp boundary_extrapolation_tolerance_max, fp x); static int run_batch_tests(void); @@ -56,152 +74,135 @@ struct args_results { /* args */ int molecule_size; - fp out_of_range_tolerance_min; - fp out_of_range_tolerance_max; + fp boundary_off_centering_tolerance_min; + fp boundary_off_centering_tolerance_max; + fp boundary_extrapolation_tolerance_min; + fp boundary_extrapolation_tolerance_max; fp x; /* results */ + int status; int i_center; fp x_rel; - int m_min, m_max; }; /* test data */ +#define OK 0 +#define X_LT_MIN MOLECULE_POSN_ERROR_X_LT_MIN +#define X_GT_MAX MOLECULE_POSN_ERROR_X_GT_MAX static const struct args_results test_data[] = { - /* molecule size 2 */ - { 2, 1.0, 1.0e-12, 7.19, INT_MIN, -1.1, 0,+1 }, - { 2, 1.0, 1.0e-12, 7.21, 42, -0.9, 0,+1 }, - { 2, 1.0, 1.0e-12, 7.24, 42, -0.6, 0,+1 }, - { 2, 1.0, 1.0e-12, 7.26, 42, -0.4, 0,+1 }, - { 2, 1.0, 1.0e-12, 7.29, 42, -0.1, 0,+1 }, - { 2, -1.0, 1.0e-12, 7.29, INT_MIN, -0.1, 0,+1 }, - { 2, 1.0e-12, 1.0e-12, 7.3, 42, 0.0, 0,+1 },/* grid_x_min */ - { 2, 1.0e-12, 1.0e-12, 7.31, 42, +0.1, 0,+1 }, - { 2, -1.0 , 1.0e-12, 7.31, 42, +0.1, 0,+1 }, - { 2, 1.0e-12, 1.0e-12, 7.35, 42, +0.5, 0,+1 }, - { 2, 1.0e-12, 1.0e-12, 7.39, 42, +0.9, 0,+1 }, - { 2, 1.0e-12, 1.0e-12, 7.41, 43, 0.1, 0,+1 }, - { 2, 1.0e-12, 1.0e-12, 9.81, 67, +0.1, 0,+1 }, - { 2, 1.0e-12, 1.0e-12, 9.85, 67, +0.5, 0,+1 }, - { 2, 1.0e-12, 1.0e-12, 9.89, 67, +0.9, 0,+1 }, - { 2, 1.0e-12, 1.0e-12, 13.45, 103, +0.5, 0,+1 }, - { 2, 1.0e-12, 1.0e-12, 13.51, 104, +0.1, 0,+1 }, - { 2, 1.0e-12, 1.0e-12, 13.55, 104, +0.5, 0,+1 }, - { 2, 1.0e-12, 1.0e-12, 13.59, 104, +0.9, 0,+1 }, - { 2, 1.0e-12, -1.0, 13.59, 104, +0.9, 0,+1 }, - { 2, 1.0e-12, 1.0e-12, 13.6, 104, +1.0, 0,+1 },/* grid_x_max */ - { 2, 1.0e-12, 1.0, 13.61, 104, +1.1, 0,+1 }, - { 2, 1.0e-12, -1.0, 13.61, INT_MAX, +1.1, 0,+1 }, - { 2, 1.0e-12, 1.0, 13.65, 104, +1.5, 0,+1 }, - { 2, 1.0e-12, 1.0, 13.69, 104, +1.9, 0,+1 }, - { 2, 1.0e-12, 1.0, 13.71, INT_MAX, +2.1, 0,+1 }, - - /* molecule size 3 */ - { 3, 1.0, 1.0e-12, 7.19, INT_MIN, -2.1, -1,+1 }, - { 3, 1.0, 1.0e-12, 7.21, 43, -1.9, -1,+1 }, - { 3, 1.0, 1.0e-12, 7.24, 43, -1.6, -1,+1 }, - { 3, 1.0, 1.0e-12, 7.26, 43, -1.4, -1,+1 }, - { 3, 1.0, 1.0e-12, 7.29, 43, -1.1, -1,+1 }, - { 3, -1.0, 1.0e-12, 7.29, INT_MIN, -1.1, -1,+1 }, - { 3, 1.0e-12, 1.0e-12, 7.3, 43, -1.0, -1,+1 },/* grid_x_min */ - { 3, 1.0e-12, 1.0e-12, 7.31, 43, -0.9, -1,+1 }, - { 3, 1.0e-12, 1.0e-12, 7.34, 43, -0.6, -1,+1 }, - { 3, -1.0, 1.0e-12, 7.34, INT_MIN, -0.6, -1,+1 }, - { 3, 1.0e-12, 1.0e-12, 7.36, 43, -0.4, -1,+1 }, - { 3, -1.0, 1.0e-12, 7.36, 43, -0.4, -1,+1 }, - { 3, 1.0e-12, 1.0e-12, 7.39, 43, -0.1, -1,+1 }, - { 3, 1.0e-12, 1.0e-12, 7.4, 43, 0.0, -1,+1 }, - { 3, 1.0e-12, 1.0e-12, 7.44, 43, +0.4, -1,+1 }, - { 3, 1.0e-12, 1.0e-12, 7.46, 44, -0.4, -1,+1 }, - { 3, 1.0e-12, 1.0e-12, 9.8, 67, 0.0, -1,+1 }, - { 3, 1.0e-12, 1.0e-12, 9.81, 67, +0.1, -1,+1 }, - { 3, 1.0e-12, 1.0e-12, 9.84, 67, +0.4, -1,+1 }, - { 3, 1.0e-12, 1.0e-12, 9.86, 68, -0.4, -1,+1 }, - { 3, 1.0e-12, 1.0e-12, 9.89, 68, -0.1, -1,+1 }, - { 3, 1.0e-12, 1.0e-12, 9.9, 68, 0.0, -1,+1 }, - { 3, 1.0e-12, 1.0e-12, 13.44, 103, +0.4, -1,+1 }, - { 3, 1.0e-12, 1.0e-12, 13.46, 104, -0.4, -1,+1 }, - { 3, 1.0e-12, 1.0e-12, 13.5, 104, 0.0, -1,+1 }, - { 3, 1.0e-12, 1.0e-12, 13.51, 104, +0.1, -1,+1 }, - { 3, 1.0e-12, 1.0e-12, 13.54, 104, +0.4, -1,+1 }, - { 3, 1.0e-12, -1.0, 13.54, 104, +0.4, -1,+1 }, - { 3, 1.0e-12, 1.0e-12, 13.56, 104, +0.6, -1,+1 }, - { 3, 1.0e-12, -1.0, 13.56, INT_MAX, +0.6, -1,+1 }, - { 3, 1.0e-12, 1.0e-12, 13.59, 104, +0.9, -1,+1 }, - { 3, 1.0e-12, 1.0e-12, 13.6, 104, +1.0, -1,+1 },/* grid_x_max */ - { 3, 1.0e-12, 1.0, 13.61, 104, +1.1, -1,+1 }, - { 3, 1.0e-12, 1.0, 13.65, 104, +1.5, -1,+1 }, - { 3, 1.0e-12, 1.0, 13.69, 104, +1.9, -1,+1 }, - { 3, 1.0e-12, 1.0, 13.71, INT_MAX, +2.1, -1,+1 }, - - /* molecule size 4 */ - { 4, 0.2, 1.0e-12, 7.27, INT_MIN, -1.3, -1,+2 }, - { 4, 0.2, 1.0e-12, 7.29, 43, -1.1, -1,+2 }, - { 4, 1.0e-12, 1.0e-12, 7.3, 43, -1.0, -1,+2 },/* grid_x_min */ - { 4, 1.0e-12, 1.0e-12, 7.33, 43, -0.7, -1,+2 }, - { 4, 1.0e-12, 1.0e-12, 7.39, 43, -0.1, -1,+2 }, - { 4, -1.0, 1.0e-12, 7.39, INT_MIN, -0.1, -1,+2 }, - { 4, 1.0e-12, 1.0e-12, 7.4, 43, 0.0, -1,+2 }, - { 4, -1.0, 1.0e-12, 7.41, 43, +0.1, -1,+2 }, - { 4, 1.0e-12, 1.0e-12, 7.42, 43, +0.2, -1,+2 }, - { 4, 1.0e-12, 1.0e-12, 7.48, 43, +0.8, -1,+2 }, - { 4, 1.0e-12, 1.0e-12, 7.51, 44, +0.1, -1,+2 }, - { 4, 1.0e-12, 1.0e-12, 9.81, 67, +0.1, -1,+2 }, - { 4, 1.0e-12, 1.0e-12, 9.85, 67, +0.5, -1,+2 }, - { 4, 1.0e-12, 1.0e-12, 9.89, 67, +0.9, -1,+2 }, - { 4, 1.0e-12, 1.0e-12, 13.39, 102, +0.9, -1,+2 }, - { 4, 1.0e-12, 1.0e-12, 13.41, 103, +0.1, -1,+2 }, - { 4, 1.0e-12, 1.0e-12, 13.48, 103, +0.8, -1,+2 }, - { 4, 1.0e-12, -1.0, 13.48, 103, +0.8, -1,+2 }, - { 4, 1.0e-12, 1.0e-12, 13.5, 103, +1.0, -1,+2 }, - { 4, 1.0e-12, 1.0e-12, 13.51, 103, +1.1, -1,+2 }, - { 4, 1.0e-12, -1.0, 13.51, INT_MAX, +1.1, -1,+2 }, - { 4, 1.0e-12, 1.0e-12, 13.55, 103, +1.5, -1,+2 }, - { 4, 1.0e-12, 1.0e-12, 13.59, 103, +1.9, -1,+2 }, - { 4, 1.0e-12, 1.0e-12, 13.6, 103, +2.0, -1,+2 },/* grid_x_max */ - { 4, 1.0e-12, 2.0, 13.79, 103, +3.9, -1,+2 }, - { 4, 1.0e-12, 2.0, 13.81, INT_MAX, +4.1, -1,+2 }, - - /* molecule size 5 */ - { 5, 3.0, 1.0e-12, 6.99, INT_MIN, -5.1, -2,+2 }, - { 5, 3.0, 1.0e-12, 7.01, 44, -4.9, -2,+2 }, - { 5, 1.0e-12, 1.0e-12, 7.3, 44, -2.0, -2,+2 },/* grid_x_min */ - { 5, 1.0e-12, 1.0e-12, 7.4, 44, -1.0, -2,+2 }, - { 5, 1.0e-12, 1.0e-12, 7.44, 44, -0.6, -2,+2 }, - { 5, -1.0, 1.0e-12, 7.44, INT_MIN, -0.6, -2,+2 }, - { 5, 1.0e-12, 1.0e-12, 7.46, 44, -0.4, -2,+2 }, - { 5, -1.0, 1.0e-12, 7.46, 44, -0.4, -2,+2 }, - { 5, 1.0e-12, 1.0e-12, 7.49, 44, -0.1, -2,+2 }, - { 5, 1.0e-12, 1.0e-12, 7.5, 44, 0.0, -2,+2 }, - { 5, 1.0e-12, 1.0e-12, 7.54, 44, +0.4, -2,+2 }, - { 5, 1.0e-12, 1.0e-12, 7.56, 45, -0.4, -2,+2 }, - { 5, 1.0e-12, 1.0e-12, 7.6, 45, 0.0, -2,+2 }, - { 5, 1.0e-12, 1.0e-12, 7.64, 45, +0.4, -2,+2 }, - { 5, 1.0e-12, 1.0e-12, 9.8, 67, 0.0, -2,+2 }, - { 5, 1.0e-12, 1.0e-12, 9.81, 67, +0.1, -2,+2 }, - { 5, 1.0e-12, 1.0e-12, 9.84, 67, +0.4, -2,+2 }, - { 5, 1.0e-12, 1.0e-12, 9.86, 68, -0.4, -2,+2 }, - { 5, 1.0e-12, 1.0e-12, 9.89, 68, -0.1, -2,+2 }, - { 5, 1.0e-12, 1.0e-12, 9.9, 68, 0.0, -2,+2 }, - { 5, 1.0e-12, 1.0e-12, 13.34, 102, +0.4, -2,+2 }, - { 5, 1.0e-12, 1.0e-12, 13.36, 103, -0.4, -2,+2 }, - { 5, 1.0e-12, 1.0e-12, 13.39, 103, -0.1, -2,+2 }, - { 5, 1.0e-12, 1.0e-12, 13.41, 103, +0.1, -2,+2 }, - { 5, 1.0e-12, 1.0e-12, 13.44, 103, +0.4, -2,+2 }, - { 5, 1.0e-12, -1.0, 13.44, 103, +0.4, -2,+2 }, - { 5, 1.0e-12, 1.0e-12, 13.46, 103, +0.6, -2,+2 }, - { 5, 1.0e-12, -1.0, 13.46, INT_MAX, +0.6, -2,+2 }, - { 5, 1.0e-12, 1.0e-12, 13.48, 103, +0.8, -2,+2 }, - { 5, 1.0e-12, 1.0e-12, 13.5, 103, +1.0, -2,+2 }, - { 5, 1.0e-12, 1.0e-12, 13.51, 103, +1.1, -2,+2 }, - { 5, 1.0e-12, 1.0e-12, 13.54, 103, +1.4, -2,+2 }, - { 5, 1.0e-12, -1.0, 13.54, INT_MAX, +1.4, -2,+2 }, - { 5, 1.0e-12, 1.0e-12, 13.56, 103, +1.6, -2,+2 }, - { 5, 1.0e-12, -1.0, 13.56, INT_MAX, +1.6, -2,+2 }, - { 5, 1.0e-12, 1.0e-12, 13.59, 103, +1.9, -2,+2 }, - { 5, 1.0e-12, 1.0e-12, 13.6, 103, +2.0, -2,+2 },/* grid_x_max */ - { 5, 1.0e-12, 1.5, 13.74, 103, +3.4, -2,+2 }, - { 5, 1.0e-12, 1.5, 13.76, INT_MAX, +3.6, -2,+2 }, + /*** molecule size 2 ***/ + /* status */ + /* off_centering extrapolation i_center */ + /* min max min max x x_rel */ + { 2, 999.0, 999.0, 1.0, 999.0, 7.19, X_LT_MIN, 0, 0.0 }, + { 2, 999.0, 999.0, 1.0, 999.0, 7.21, OK, 42, -0.9 }, + { 2, 999.0, 999.0, 1.0, 999.0, 7.24, OK, 42, -0.6 }, + { 2, 999.0, 999.0, 1.0, 999.0, 7.26, OK, 42, -0.4 }, + { 2, 0.0, 0.0, 1.0, 999.0, 7.26, X_LT_MIN, 0, 0.0 }, + { 2, 1.0e-6, 999.0, 0.2, 999.0, 7.29, X_LT_MIN, 0, 0.0 }, + { 2, 0.2, 999.0, 1.0e-6, 999.0, 7.29, X_LT_MIN, 0, 0.0 }, + { 2, 0.2, 999.0, 0.2, 999.0, 7.29, OK, 42, -0.1 }, + { 2, 0.0, 0.0, 1.0e-6, 999.0, 7.31, OK, 42, +0.1 }, + { 2, 0.0, 0.0, 1.0e-6, 999.0, 7.34, OK, 42, +0.4 }, + { 2, 0.0, 0.0, 1.0e-6, 999.0, 7.36, OK, 42, +0.6 }, + { 2, 0.0, 0.0, 1.0e-6, 999.0, 7.39, OK, 42, +0.9 }, + { 2, 0.0, 0.0, 1.0e-6, 999.0, 7.41, OK, 43, +0.1 }, + { 2, 0.0, 0.0, 1.0e-6, 1.0e-6, 9.81, OK, 67, +0.1 }, + { 2, 0.0, 0.0, 1.0e-6, 1.0e-6, 9.84, OK, 67, +0.4 }, + { 2, 0.0, 0.0, 1.0e-6, 1.0e-6, 9.85, OK, 67, +0.5 }, + { 2, 0.0, 0.0, 1.0e-6, 1.0e-6, 9.86, OK, 67, +0.6 }, + { 2, 0.0, 0.0, 1.0e-6, 1.0e-6, 9.89, OK, 67, +0.9 }, + { 2, 0.0, 0.0, 1.0e-6, 1.0e-6, 13.45, OK, 103, +0.5 }, + { 2, 0.0, 0.0, 1.0e-6, 1.0e-6, 13.51, OK, 104, +0.1 }, + { 2, 0.0, 0.0, 1.0e-6, 1.0e-6, 13.59, OK, 104, +0.9 }, + { 2, 0.0, 0.2, 1.0e-6, 1.0e-6, 13.61, X_GT_MAX, 0, 0.0 }, + { 2, 0.0, 1.0e-6, 1.0e-6, 1.0e-6, 13.61, X_GT_MAX, 0, 0.0 }, + { 2, 0.0, 0.2, 0.0, 0.2, 13.61, OK, 104, +1.1 }, + + /*** molecule size 3 ***/ + /* status */ + /* off_centering extrapolation i_center */ + /* min max min max x x_rel */ + { 3, 0.5, 999.0, 999.0, 999.0, 7.29, X_LT_MIN, 0, 0.0 }, + { 3, 999.0, 999.0, 0.0, 999.0, 7.29, X_LT_MIN, 0, 0.0 }, + { 3, 0.7, 999.0, 0.2, 999.0, 7.29, OK, 43, -1.1 }, + { 3, 0.3, 999.0, 0.0, 999.0, 7.31, X_LT_MIN, 0, 0.0 }, + { 3, 0.5, 999.0, 0.0, 999.0, 7.31, OK, 43, -0.9 }, + { 3, 0.0, 0.0, 0.0, 0.0, 7.36, OK, 43, -0.4 }, + { 3, 0.0, 0.0, 0.0, 0.0, 7.44, OK, 43, +0.4 }, + { 3, 0.0, 0.0, 0.0, 0.0, 7.46, OK, 44, -0.4 }, + { 3, 0.0, 0.0, 0.0, 0.0, 9.81, OK, 67, +0.1 }, + { 3, 0.0, 0.0, 0.0, 0.0, 9.84, OK, 67, +0.4 }, + { 3, 0.0, 0.0, 0.0, 0.0, 9.86, OK, 68, -0.4 }, + { 3, 0.0, 0.0, 0.0, 0.0, 9.89, OK, 68, -0.1 }, + { 3, 0.0, 0.0, 0.0, 0.0, 13.44, OK, 103, +0.4 }, + { 3, 0.0, 0.0, 0.0, 0.0, 13.46, OK, 104, -0.4 }, + { 3, 0.0, 0.0, 0.0, 0.0, 13.54, OK, 104, +0.4 }, + { 3, 0.0, 0.2, 0.0, 0.0, 13.56, OK, 104, +0.6 }, + { 3, 999.0, 0.0, 999.0, 999.0, 13.56, X_GT_MAX, 0, 0.0 }, + { 3, 0.0, 0.5, 0.0, 0.0, 13.59, OK, 104, +0.9 }, + { 3, 999.0, 0.3, 999.0, 999.0, 13.59, X_GT_MAX, 0, 0.0 }, + { 3, 999.0, 0.5, 999.0, 999.0, 13.61, X_GT_MAX, 0, 0.0 }, + { 3, 999.0, 999.0, 999.0, 0.0, 13.61, X_GT_MAX, 0, 0.0 }, + { 3, 999.0, 0.7, 999.0, 0.2, 13.61, OK, 104, +1.1 }, + { 3, 999.0, 999.0, 999.0, 0.2, 13.63, X_GT_MAX, 0, 0.0 }, + { 3, 999.0, 0.7, 999.0, 999.0, 13.63, X_GT_MAX, 0, 0.0 }, + { 3, 999.0, 0.9, 999.0, 0.4, 13.63, OK, 104, +1.3 }, + + /*** molecule size 4 ***/ + /* status */ + /* off_centering extrapolation i_center */ + /* min max min max x x_rel */ + { 4, 1.0, 999.0, 999.0, 999.0, 7.29, X_LT_MIN, 0, 0.0 }, + { 4, 999.0, 999.0, 0.0, 999.0, 7.29, X_LT_MIN, 0, 0.0 }, + { 4, 1.2, 999.0, 0.2, 999.0, 7.29, OK, 43, -1.1 }, + { 4, 0.5, 999.0, 0.0, 999.0, 7.31, X_LT_MIN, 0, 0.0 }, + { 4, 1.0, 999.0, 0.0, 999.0, 7.31, OK, 43, -0.9 }, + { 4, 0.0, 999.0, 0.0, 999.0, 7.39, X_LT_MIN, 0, 0.0 }, + { 4, 0.2, 999.0, 0.0, 999.0, 7.39, OK, 43, -0.1 }, + { 4, 0.0, 0.0, 0.0, 999.0, 7.41, OK, 43, +0.1 }, + { 4, 0.0, 0.0, 0.0, 999.0, 7.49, OK, 43, +0.9 }, + { 4, 0.0, 0.0, 0.0, 999.0, 7.51, OK, 44, +0.1 }, + { 4, 0.0, 0.0, 0.0, 0.0, 9.81, OK, 67, +0.1 }, + { 4, 0.0, 0.0, 0.0, 0.0, 9.84, OK, 67, +0.4 }, + { 4, 0.0, 0.0, 0.0, 0.0, 9.85, OK, 67, +0.5 }, + { 4, 0.0, 0.0, 0.0, 0.0, 9.86, OK, 67, +0.6 }, + { 4, 0.0, 0.0, 0.0, 0.0, 9.89, OK, 67, +0.9 }, + { 4, 0.0, 0.0, 0.0, 0.0, 13.44, OK, 103, +0.4 }, + { 4, 0.0, 0.0, 0.0, 0.0, 13.49, OK, 103, +0.9 }, + { 4, 0.0, 0.2, 0.0, 0.0, 13.51, OK, 103, +1.1 }, + { 4, 999.0, 0.0, 999.0, 999.0, 13.51, X_GT_MAX, 0, 0.0 }, + { 4, 999.0, 0.8, 999.0, 0.0, 13.59, X_GT_MAX, 0, 0.0 }, + { 4, 0.0, 1.0, 0.0, 0.0, 13.59, OK, 103, +1.9 }, + { 4, 999.0, 999.0, 999.0, 0.0, 13.61, X_GT_MAX, 0, 0.0 }, + { 4, 0.0, 1.0, 999.0, 999.0, 13.61, X_GT_MAX, 0, 0.0 }, + { 4, 0.0, 1.2, 0.0, 0.2, 13.61, OK, 103, +2.1 }, + + /*** molecule size 5 ***/ + /* status */ + /* off_centering extrapolation i_center */ + /* min max min max x x_rel */ + { 5, 1.5, 999.0, 999.0, 999.0, 7.29, X_LT_MIN, 0, 0.0 }, + { 5, 999.0, 999.0, 0.0, 999.0, 7.29, X_LT_MIN, 0, 0.0 }, + { 5, 1.7, 999.0, 0.2, 999.0, 7.29, OK, 44, -2.1 }, + { 5, 0.9, 999.0, 0.0, 999.0, 7.35, X_LT_MIN, 0, 0.0 }, + { 5, 1.1, 999.0, 0.0, 999.0, 7.35, OK, 44, -1.5 }, + { 5, 0.0, 999.0, 999.0, 999.0, 7.44, X_LT_MIN, 0, 0.0 }, + { 5, 0.2, 999.0, 0.0, 999.0, 7.44, OK, 44, -0.6 }, + { 5, 0.0, 0.0, 0.0, 0.0, 7.46, OK, 44, -0.4 }, + { 5, 0.0, 0.0, 0.0, 0.0, 9.81, OK, 67, +0.1 }, + { 5, 0.0, 0.0, 0.0, 0.0, 9.84, OK, 67, +0.4 }, + { 5, 0.0, 0.0, 0.0, 0.0, 9.86, OK, 68, -0.4 }, + { 5, 0.0, 0.0, 0.0, 0.0, 9.89, OK, 68, -0.1 }, + { 5, 0.0, 0.0, 0.0, 0.0, 13.44, OK, 103, +0.4 }, + { 5, 999.0, 0.0, 999.0, 999.0, 13.46, X_GT_MAX, 0, 0.0 }, + { 5, 0.0, 0.2, 0.0, 0.0, 13.46, OK, 103, +0.6 }, + { 5, 0.0, 1.1, 0.0, 0.0, 13.55, OK, 103, +1.5 }, + { 5, 999.0, 0.9, 999.0, 999.0, 13.55, X_GT_MAX, 0, 0.0 }, + { 5, 999.0, 999.0, 999.0, 0.0, 13.61, X_GT_MAX, 0, 0.0 }, + { 5, 999.0, 1.5, 999.0, 999.0, 13.61, X_GT_MAX, 0, 0.0 }, + { 5, 999.0, 1.7, 999.0, 0.2, 13.61, OK, 103, +2.1 }, }; #define N_TESTS ((int) (sizeof(test_data)/sizeof(test_data[0]))) @@ -212,7 +213,10 @@ int main(int argc, const char *const argv[]) { bool N_fail; int molecule_size; -double out_of_range_tolerance_min, out_of_range_tolerance_max; +fp boundary_off_centering_tolerance_min; +fp boundary_off_centering_tolerance_max; +fp boundary_extrapolation_tolerance_min; +fp boundary_extrapolation_tolerance_max; double x; switch (argc) @@ -222,44 +226,34 @@ case 1: N_fail = run_batch_tests(); if (N_fail == 0) then { - printf("*** all tests successful ***\n"); + printf("*** all %d tests successful ***\n", N_TESTS); return 0; } else { - printf("*** %d test(s) failed ***\n", N_fail); + printf("*** %d/%d test(s) failed ***\n", N_fail, N_TESTS); return 1; } -case 5: +case 7: if ( (sscanf(argv[1], "%d", &molecule_size) == 1) - && (sscanf(argv[2], "%lf", &out_of_range_tolerance_min) == 1) - && (sscanf(argv[3], "%lf", &out_of_range_tolerance_max) == 1) - && (sscanf(argv[4], "%lf", &x) == 1) ) + && (sscanf(argv[2], "%lf", &boundary_off_centering_tolerance_min) == 1) + && (sscanf(argv[3], "%lf", &boundary_off_centering_tolerance_max) == 1) + && (sscanf(argv[4], "%lf", &boundary_extrapolation_tolerance_min) == 1) + && (sscanf(argv[5], "%lf", &boundary_extrapolation_tolerance_max) == 1) + && (sscanf(argv[6], "%lf", &x) == 1) ) then { run_interactive_test(molecule_size, - out_of_range_tolerance_min, - out_of_range_tolerance_max, + boundary_off_centering_tolerance_min, + boundary_off_centering_tolerance_max, + boundary_extrapolation_tolerance_min, + boundary_extrapolation_tolerance_max, x); - return 0; + return 0; /*** NORMAL EXIT ***/ } - /* fall through */ + /* error ==> fall through */ default: - fprintf(stderr, - "usage:\n" - "# run a single test as specified:\n" - " %s molecule_size \\\n" - " %*s out_of_range_tolerance_min \\\n" - " %*s out_of_range_tolerance_max \\\n" - " %*s x\n" - "# run a preset set of tests:\n" - " %s\n" - , - argv[0], - (int) strlen(argv[0]), "", - (int) strlen(argv[0]), "", - (int) strlen(argv[0]), "", - argv[0]); + fprintf(stderr, help_msg); return 1; } } @@ -271,33 +265,39 @@ default: */ static void run_interactive_test(int molecule_size, - fp out_of_range_tolerance_min, - fp out_of_range_tolerance_max, + fp boundary_off_centering_tolerance_min, + fp boundary_off_centering_tolerance_max, + fp boundary_extrapolation_tolerance_min, + fp boundary_extrapolation_tolerance_max, fp x) { +int i_center; fp x_rel; -int m_min, m_max; printf("testing with molecule_size=%d\n", molecule_size); -printf(" out_of_range_tolerance_[min,max]=[%g,%g]\n", - (double) out_of_range_tolerance_min, - (double) out_of_range_tolerance_max); +printf(" boundary_off_centering_tolerance_[min,max]=[%g,%g]\n", + (double) boundary_off_centering_tolerance_min, + (double) boundary_off_centering_tolerance_max); +printf(" boundary_extrapolation_tolerance_[min,max]=[%g,%g]\n", + (double) boundary_extrapolation_tolerance_min, + (double) boundary_extrapolation_tolerance_max); printf(" x=%g\n", (double) x); { -const int i_center = LocalInterp_molecule_posn(grid_x0, grid_dx, - grid_i_min, grid_i_max, - molecule_size, - out_of_range_tolerance_min, - out_of_range_tolerance_max, - x, - &x_rel, - &m_min, &m_max); - -if ((i_center == INT_MIN) || (i_center == INT_MAX)) - then printf("i_center=%d\n", i_center); - else printf("i_center=%d x_rel=%g m_[min,max]=[%d,%d] i_[min,max]=[%d,%d]\n", - i_center, x_rel, m_min, m_max, i_center+m_min, i_center+m_max); +const int status + = LocalInterp_molecule_posn(grid_x0, grid_dx, + grid_i_min, grid_i_max, + molecule_size, + boundary_off_centering_tolerance_min, + boundary_off_centering_tolerance_max, + boundary_extrapolation_tolerance_min, + boundary_extrapolation_tolerance_max, + x, + &i_center, &x_rel); + +if (status < 0) + then printf("status=%d\n", status); + else printf("i_center=%d x_rel=%g\n", i_center, x_rel); } } @@ -319,41 +319,48 @@ int failure_count = 0; for (i = 0 ; i < N_TESTS ; ++i) { const struct args_results *p = & test_data[i]; - fp x_rel = 0.0; - int m_min = 0, m_max = 0; - int i_center = LocalInterp_molecule_posn(grid_x0, grid_dx, - grid_i_min, grid_i_max, - p->molecule_size, - p->out_of_range_tolerance_min, - p->out_of_range_tolerance_max, - p->x, - &x_rel, - &m_min, &m_max); - bool ok = ((i_center == INT_MIN) || (i_center == INT_MAX)) - ? (i_center == p->i_center) - : ( (i_center == p->i_center) - && fuzzy_EQ(x_rel, p->x_rel) - && (m_min == p->m_min) - && (m_max == p->m_max) ); - - int msglen = - printf("size=%d tol=[%g,%g] x=%g ==> ", - p->molecule_size, - (double) p->out_of_range_tolerance_min, - (double) p->out_of_range_tolerance_max, - (double) p->x); - printf("i_center=%d x_rel=%g m=[%d,%d]\n", - i_center, (double) x_rel, m_min, m_max); - - if (! ok) - then { + int i_center = 999; + fp x_rel = 999.999; + const int status = LocalInterp_molecule_posn + (grid_x0, grid_dx, + grid_i_min, grid_i_max, + p->molecule_size, + p->boundary_off_centering_tolerance_min, + p->boundary_off_centering_tolerance_max, + p->boundary_extrapolation_tolerance_min, + p->boundary_extrapolation_tolerance_max, + p->x, + &i_center, &x_rel); + bool ok = (status == 0) + ? ( (status == p->status) && (i_center == p->i_center) + && fuzzy_EQ(x_rel, p->x_rel) ) + : (status == p->status); + + printf("size=%d off_centering_tol=[%g,%g] extrapolation_tol=[%g,%g] x=%g", + p->molecule_size, + (double) p->boundary_off_centering_tolerance_min, + (double) p->boundary_off_centering_tolerance_max, + (double) p->boundary_extrapolation_tolerance_min, + (double) p->boundary_extrapolation_tolerance_max, + (double) p->x); + printf(": status=%d", status); + if (status == 0) + then printf(" i_center=%d x_rel=%g: ", + i_center, (double) x_rel); + + if (ok) + then printf(": ok\n"); + else { ++failure_count; - { - int error_msglen = printf("***** FAIL: "); - printf("%-*s", msglen-error_msglen, "expected"); - printf("i_center=%d x_rel=%g m_[min,max]=[%d,%d]\n", - p->i_center, (double) p->x_rel, p->m_min, p->m_max); - } + printf(": FAIL\n"); + if (p->status == 0) + then printf("*** expected status=%d i_center=%d x_rel=%g\n", + p->status, p->i_center, (double) p->x_rel); + else printf("*** expected status=%d\n", p->status); + if (status == 0) + then printf("*** got status=%d i_center=%d x_rel=%g\n", + status, i_center, (double) x_rel); + else printf("*** got status=%d\n", status); } } -- cgit v1.2.3