diff options
Diffstat (limited to 'src/GeneralizedPolynomial-Uniform/InterpLocalUniform.c')
| -rw-r--r-- | src/GeneralizedPolynomial-Uniform/InterpLocalUniform.c | 238 |
1 files changed, 190 insertions, 48 deletions
diff --git a/src/GeneralizedPolynomial-Uniform/InterpLocalUniform.c b/src/GeneralizedPolynomial-Uniform/InterpLocalUniform.c index d0b4c78..9e0bf9f 100644 --- a/src/GeneralizedPolynomial-Uniform/InterpLocalUniform.c +++ b/src/GeneralizedPolynomial-Uniform/InterpLocalUniform.c @@ -3,8 +3,9 @@ /* ** *****data structures and functions local to this file ***** ** - * LocalInterp_ILU_Lagrange - external entry point for Lagrange interpolator - * LocalInterp_ILU_Hermite - external entry point for Hermite interpolator + * LocalInterp_ILU_Lagrange_TP - driver for Lagrange interpolator (tensor prod) + * LocalInterp_ILU_Lagrange_MD - driver for Lagrange interpolator (max degree) + * LocalInterp_ILU_Hermite - driver point for Hermite interpolator ** ** LocalInterp_InterpLocalUniform - main driver routine ** @@ -56,7 +57,8 @@ #include "util_Table.h" #include "InterpLocalUniform.h" -#include "Lagrange/all_prototypes.h" +#include "Lagrange-tensor-product/all_prototypes.h" +#include "Lagrange-maximum-degree/all_prototypes.h" #include "Hermite/all_prototypes.h" /* the rcs ID and its dummy function to use it */ @@ -78,8 +80,11 @@ CCTK_FILEVERSION(CactusBase_LocalInterp_src_GeneralizedPolynomialUniform_InterpL /* this enum describes which interpolation operator we're using */ enum interp_operator { - interp_operator_Lagrange, /* "Lagrange polynomial interpolation" */ - interp_operator_Hermite, /* "Hermite polynomial interpolation" */ + interp_operator_Lagrange_TP, /* "Lagrange polynomial interpolation */ + /* (tensor product)" */ + interp_operator_Lagrange_MD, /* "Lagrange polynomial interpolation + /* (maximum degree)" */ + interp_operator_Hermite, /* "Hermite polynomial interpolation" */ N_INTERP_OPERATORS /* this must be the last entry in the enum */ }; @@ -185,6 +190,8 @@ typedef * 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.) + * + * In the comments on the initializers, "n.i." = "not implemented". */ static const p_interp_fn_t p_interp_fn_table[N_INTERP_OPERATORS] [MAX_N_DIMS+1] /* see above */ @@ -193,9 +200,69 @@ static const p_interp_fn_t p_interp_fn_table[N_INTERP_OPERATORS] [MAX_SMOOTHING+1] = { - /* in the comments on the initializers, "n.i." = "not implemented" */ { - /* interp_operator == interp_operator_Lagrange */ + /* interp_operator == interp_operator_Lagrange_TP */ + { + /* 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_U_LagrTP_1d_cube10, }, /* order=1 */ + { LocalInterp_U_LagrTP_1d_cube20, }, /* order=2 */ + { LocalInterp_U_LagrTP_1d_cube30, }, /* order=3 */ + { LocalInterp_U_LagrTP_1d_cube40, }, /* order=4 */ + { LocalInterp_U_LagrTP_1d_cube50, }, /* order=5 */ + { LocalInterp_U_LagrTP_1d_cube60, }, /* order=6 */ + }, + }, + + { + /* N_dims = 2 */ + { + /* molecule_family = molecule_family_cube */ + { NULL_INTERP_FN_PTR, }, /* order=0 (unused) */ + { LocalInterp_U_LagrTP_2d_cube10, }, /* order=1 */ + { LocalInterp_U_LagrTP_2d_cube20, }, /* order=2 */ + { LocalInterp_U_LagrTP_2d_cube30, }, /* order=3 */ + { LocalInterp_U_LagrTP_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_U_LagrTP_3d_cube10, }, /* order=1 */ + { LocalInterp_U_LagrTP_3d_cube20, }, /* order=2 */ + { LocalInterp_U_LagrTP_3d_cube30, }, /* order=3 */ + { LocalInterp_U_LagrTP_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_TP */ + }, + + { + /* interp_operator == interp_operator_Lagrange_MD */ { /* N_dims = 0 ==> unused */ { @@ -215,12 +282,12 @@ static const p_interp_fn_t p_interp_fn_table[N_INTERP_OPERATORS] { /* 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 */ + { LocalInterp_U_LagrMD_1d_cube10, }, /* order=1 */ + { LocalInterp_U_LagrMD_1d_cube20, }, /* order=2 */ + { LocalInterp_U_LagrMD_1d_cube30, }, /* order=3 */ + { LocalInterp_U_LagrMD_1d_cube40, }, /* order=4 */ + { LocalInterp_U_LagrMD_1d_cube50, }, /* order=5 */ + { LocalInterp_U_LagrMD_1d_cube60, }, /* order=6 */ }, }, @@ -229,10 +296,10 @@ static const p_interp_fn_t p_interp_fn_table[N_INTERP_OPERATORS] { /* 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 */ + { LocalInterp_U_LagrMD_2d_cube10, }, /* order=1 */ + { LocalInterp_U_LagrMD_2d_cube20, }, /* order=2 */ + { LocalInterp_U_LagrMD_2d_cube30, }, /* order=3 */ + { LocalInterp_U_LagrMD_2d_cube40, }, /* order=4 */ { NULL_INTERP_FN_PTR, }, /* order=5 (n.i.) */ { NULL_INTERP_FN_PTR, }, /* order=6 (n.i.) */ }, @@ -243,16 +310,16 @@ static const p_interp_fn_t p_interp_fn_table[N_INTERP_OPERATORS] { /* 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 */ + { LocalInterp_U_LagrMD_3d_cube10, }, /* order=1 */ + { LocalInterp_U_LagrMD_3d_cube20, }, /* order=2 */ + { LocalInterp_U_LagrMD_3d_cube30, }, /* order=3 */ + { LocalInterp_U_LagrMD_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 */ + /* end of interp_operator == interp_operator_Lagrange_MD */ }, { @@ -314,7 +381,7 @@ static const p_interp_fn_t p_interp_fn_table[N_INTERP_OPERATORS] }, /* end of interp_operator == interp_operator_Hermite */ - } + }, }; @@ -323,7 +390,7 @@ static const p_interp_fn_t p_interp_fn_table[N_INTERP_OPERATORS] /******************************************************************************/ /*@@ - @routine LocalInterp_ILU_Lagrange + @routine LocalInterp_ILU_Lagrange_TP @date 22 Oct 2001 @author Jonathan Thornburg <jthorn@aei.mpg.de> @desc @@ -332,10 +399,85 @@ static const p_interp_fn_t p_interp_fn_table[N_INTERP_OPERATORS] 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. + interpolation points. For multidimensional interpolation, it + uses the tensor-product basis for the interpolating function, + so the interpolant is guaranteed to be continuous and to pass + through the input data at the grid points (up to floating-point + roundoff errors). + + This interpolator can also + * do differentiation and/or smoothing in the process of the + interpolation + * 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_TP(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_TP, + "Lagrange polynomial interpolation (tensor product)", + LAGRANGE_BNDRY_OFF_CNTR_TOL_DEF, + LAGRANGE_BNDRY_EXTRAP_TOL_DEF, + 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_Lagrange_MD + @date 22 Oct 2001 + @author Jonathan Thornburg <jthorn@aei.mpg.de> + @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. For multidimensional interpolation, it + uses the maximum-degree basis for the interpolating function. + This means the interpolant is in general *not* continuous and + in general does *not* pass through the input data at the grid + points. + + This interpolator can also + * do differentiation and/or smoothing in the process of the + interpolation + * 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 @@ -343,27 +485,27 @@ static const p_interp_fn_t p_interp_fn_table[N_INTERP_OPERATORS] 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[]) +int LocalInterp_ILU_Lagrange_MD(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", +return LocalInterp_InterpLocalUniform(interp_operator_Lagrange_MD, + "Lagrange polynomial interpolation (maximum degree)", LAGRANGE_BNDRY_OFF_CNTR_TOL_DEF, LAGRANGE_BNDRY_EXTRAP_TOL_DEF, N_dims, |