aboutsummaryrefslogtreecommitdiff
path: root/src/GeneralizedPolynomial-Uniform/InterpLocalUniform.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/GeneralizedPolynomial-Uniform/InterpLocalUniform.c')
-rw-r--r--src/GeneralizedPolynomial-Uniform/InterpLocalUniform.c408
1 files changed, 328 insertions, 80 deletions
diff --git a/src/GeneralizedPolynomial-Uniform/InterpLocalUniform.c b/src/GeneralizedPolynomial-Uniform/InterpLocalUniform.c
index a137b06..f6982ac 100644
--- a/src/GeneralizedPolynomial-Uniform/InterpLocalUniform.c
+++ b/src/GeneralizedPolynomial-Uniform/InterpLocalUniform.c
@@ -17,9 +17,9 @@
Since this ultimately yields _some_ linear combination of the
data values, we precompute the coefficients for efficiency.
This is done with separate Maple-generated formulas for each
- combination of number of dimensions, molecule family, order,
- amount of Savitzky-Golay smoothing to be done, and
- differentiation operator.
+ combination of number of interpolation operator, number of
+ dimensions, molecule family, order, amount of Savitzky-Golay
+ smoothing to be done, and differentiation operator.
@enddesc
@version $Id$
@@*/
@@ -34,13 +34,175 @@
#include "util_Table.h"
#include "InterpLocalUniform.h"
-#include "all-prototypes.h"
+#include "Lagrange/all_prototypes.h"
/* the rcs ID and its dummy function to use it */
static const char *rcsid = "$Header$";
CCTK_FILEVERSION(CactusBase_LocalInterp_src_GeneralizedPolynomialUniform_InterpLocalUniform_c)
-/* prototypes for static functions (= local to this file) */
+/******************************************************************************/
+
+/*
+ * *****data structures and functions local to this file *****
+ */
+
+/*
+ * data structures local to this file
+ */
+
+/* this enum describes which interpolation operator we're using */
+enum interp_operator
+ {
+ interp_operator_Lagrange, /* "Lagrange polynomial interpolation" */
+ interp_operator_Hermite, /* "Hermite polynomial interpolation" */
+ N_INTERP_OPERATORS /* this must be the last entry in the enum */
+ };
+
+
+/*
+ * prototypes for functions local to this file
+ */
+static
+ int LocalInterp_InterpLocalUniform(enum interp_operator interp_operator,
+ const char* interp_operator_string,
+ int N_dims,
+ int param_table_handle,
+ /***** coordinate system *****/
+ const CCTK_REAL coord_origin[],
+ const CCTK_REAL coord_delta[],
+ /***** interpolation points *****/
+ int N_interp_points,
+ int interp_coords_type_code,
+ const void *const interp_coords[],
+ /***** input arrays *****/
+ int N_input_arrays,
+ const CCTK_INT input_array_dims[],
+ const CCTK_INT input_array_type_codes[],
+ const void *const input_arrays[],
+ /***** output arrays *****/
+ int N_output_arrays,
+ const CCTK_INT output_array_type_codes[],
+ void *const output_arrays[]);
+
+/******************************************************************************/
+/******************************************************************************/
+/******************************************************************************/
+
+/*@@
+ @routine LocalInterp_ILU_Lagrange
+ @date 22 Oct 2001
+ @author Jonathan Thornburg <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. It can also do differentiation and/or
+ smoothing in the process of the interpolation. It can also
+ return information about the Jacobian of the interpolated
+ values with respect to the input arrays.
+
+ This function does nothing except pass all its arguments
+ down to LocalInterp_InterpLocalUniform() with an extra
+ 2 arguments added to indicate the operator handle. See that
+ function for details on this function's arguments/results.
+ @enddesc
+ @@*/
+int LocalInterp_ILU_Lagrange(int N_dims,
+ int param_table_handle,
+ /***** coordinate system *****/
+ const CCTK_REAL coord_origin[],
+ const CCTK_REAL coord_delta[],
+ /***** interpolation points *****/
+ int N_interp_points,
+ int interp_coords_type_code,
+ const void *const interp_coords[],
+ /***** input arrays *****/
+ int N_input_arrays,
+ const CCTK_INT input_array_dims[],
+ const CCTK_INT input_array_type_codes[],
+ const void *const input_arrays[],
+ /***** output arrays *****/
+ int N_output_arrays,
+ const CCTK_INT output_array_type_codes[],
+ void *const output_arrays[])
+{
+return LocalInterp_InterpLocalUniform(interp_operator_Lagrange,
+ "Lagrange polynomial interpolation",
+ N_dims,
+ param_table_handle,
+ coord_origin, coord_delta,
+ N_interp_points,
+ interp_coords_type_code,
+ interp_coords,
+ N_input_arrays,
+ input_array_dims,
+ input_array_type_codes,
+ input_arrays,
+ N_output_arrays,
+ output_array_type_codes,
+ output_arrays);
+}
+
+/******************************************************************************/
+
+/*@@
+ @routine LocalInterp_ILU_Hermite
+ @date 22 Oct 2001
+ @author Jonathan Thornburg <jthorn@aei.mpg.de>
+ @desc
+ This function is a top-level driver for
+ LocalInterp::CCTK_InterpLocalUniform().
+
+ It does Hermite interpolation of a set of input arrays defined
+ on a uniform N-dimensional tensor-product grid, to a set of
+ interpolation points. It can also do differentiation and/or
+ smoothing in the process of the interpolation. It can also
+ return information about the Jacobian of the interpolated
+ values with respect to the input arrays.
+
+ This function does nothing except pass all its arguments
+ down to LocalInterp_InterpLocalUniform() with an extra
+ 2 arguments added to indicate the operator handle. See that
+ function for details on this function's arguments/results.
+ @enddesc
+ @@*/
+int LocalInterp_ILU_Hermite(int N_dims,
+ int param_table_handle,
+ /***** coordinate system *****/
+ const CCTK_REAL coord_origin[],
+ const CCTK_REAL coord_delta[],
+ /***** interpolation points *****/
+ int N_interp_points,
+ int interp_coords_type_code,
+ const void *const interp_coords[],
+ /***** input arrays *****/
+ int N_input_arrays,
+ const CCTK_INT input_array_dims[],
+ const CCTK_INT input_array_type_codes[],
+ const void *const input_arrays[],
+ /***** output arrays *****/
+ int N_output_arrays,
+ const CCTK_INT output_array_type_codes[],
+ void *const output_arrays[])
+{
+return LocalInterp_InterpLocalUniform(interp_operator_Hermite,
+ "Hermite polynomial interpolation",
+ N_dims,
+ param_table_handle,
+ coord_origin, coord_delta,
+ N_interp_points,
+ interp_coords_type_code,
+ interp_coords,
+ N_input_arrays,
+ input_array_dims,
+ input_array_type_codes,
+ input_arrays,
+ N_output_arrays,
+ output_array_type_codes,
+ output_arrays);
+}
/******************************************************************************/
/******************************************************************************/
@@ -51,7 +213,7 @@ CCTK_FILEVERSION(CactusBase_LocalInterp_src_GeneralizedPolynomialUniform_InterpL
@date 22 Oct 2001
@author Jonathan Thornburg <jthorn@aei.mpg.de>
@desc
- This function is the top-level driver for
+ This function is the main driver for
LocalInterp::CCTK_InterpLocalUniform().
It interpolates a set of input arrays defined on a uniform
@@ -63,7 +225,8 @@ CCTK_FILEVERSION(CactusBase_LocalInterp_src_GeneralizedPolynomialUniform_InterpL
This function is just a driver: it validates arguments,
extracts optional arguments from the parameter table,
- then decodes (N_dims, molecule_family, order, smoothing)
+ then decodes
+ (interp_operator, N_dims, molecule_family, order, smoothing)
and calls the appropriate subfunction to do the actual
interpolation, then finally stores some results back in
the parameter table.
@@ -71,6 +234,17 @@ CCTK_FILEVERSION(CactusBase_LocalInterp_src_GeneralizedPolynomialUniform_InterpL
***** misc arguments *****
+ @var interp_operator
+ @vdesc describes the interpolation operator
+ @vtype enum interp_operator interp_operator
+ @endvar
+
+ @var interp_operator_string
+ @vdesc gives the character-string name of the interpolation operator
+ (this is used only for printing error messages)
+ @vtype const char* interp_operator_string
+ @endvar
+
@var N_dims
@vdesc dimensionality of the interpolation
@vtype int N_dims (must be >= 1)
@@ -412,24 +586,27 @@ CCTK_FILEVERSION(CactusBase_LocalInterp_src_GeneralizedPolynomialUniform_InterpL
functions called by this function
@endreturndesc
@@*/
-int LocalInterp_InterpLocalUniform(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[])
+static
+ int LocalInterp_InterpLocalUniform(enum interp_operator interp_operator,
+ const char* interp_operator_string,
+ int N_dims,
+ int param_table_handle,
+ /***** coordinate system *****/
+ const CCTK_REAL coord_origin[],
+ const CCTK_REAL coord_delta[],
+ /***** interpolation points *****/
+ int N_interp_points,
+ int interp_coords_type_code,
+ const void *const interp_coords[],
+ /***** input arrays *****/
+ int N_input_arrays,
+ const CCTK_INT input_array_dims[],
+ const CCTK_INT input_array_type_codes[],
+ const void *const input_arrays[],
+ /***** output arrays *****/
+ int N_output_arrays,
+ const CCTK_INT output_array_type_codes[],
+ void *const output_arrays[])
{
/*
* Implementation Note:
@@ -1186,7 +1363,7 @@ if (status)
/******************************************************************************/
/*
- * decode (N_dims, molecule_family, order, smoothing)
+ * decode (interp_operator, N_dims, molecule_family, order, smoothing)
* and call the appropriate subfunction to do the actual interpolation
*/
{
@@ -1209,71 +1386,141 @@ typedef
* only use 1...X. (Having the array oversize like this slightly simplifies
* the array subscripting.)
*/
-static const p_interp_fn_t p_interp_fn_table[MAX_N_DIMS+1] /* see above */
+static const p_interp_fn_t p_interp_fn_table[N_INTERP_OPERATORS]
+ [MAX_N_DIMS+1] /* see above */
[N_MOLECULE_FAMILIES]
[MAX_ORDER+1] /* see above */
[MAX_SMOOTHING+1]
- = {
- /* in the comments on the initializers, "n.i." = "not implemented" */
+ = {
+
+ /* in the comments on the initializers, "n.i." = "not implemented" */
+ {
+ /* interp_operator == interp_operator_Lagrange */
+ {
+ /* N_dims = 0 ==> unused */
+ {
+ /* molecule_family = molecule_family_cube */
+ { NULL_INTERP_FN_PTR , }, /* order=0 (unused) */
+ { NULL_INTERP_FN_PTR , }, /* order=1 (unused) */
+ { NULL_INTERP_FN_PTR , }, /* order=2 (unused) */
+ { NULL_INTERP_FN_PTR , }, /* order=3 (unused) */
+ { NULL_INTERP_FN_PTR , }, /* order=4 (unused) */
+ { NULL_INTERP_FN_PTR , }, /* order=5 (unused) */
+ { NULL_INTERP_FN_PTR , }, /* order=6 (unused) */
+ },
+ },
+
+ {
+ /* N_dims = 1 */
+ {
+ /* molecule_family = molecule_family_cube */
+ { NULL_INTERP_FN_PTR , }, /* order=0 (unused) */
+ { LocalInterp_ULagrange_1dcube_10, }, /* order=1 */
+ { LocalInterp_ULagrange_1dcube_20, }, /* order=2 */
+ { LocalInterp_ULagrange_1dcube_30, }, /* order=3 */
+ { LocalInterp_ULagrange_1dcube_40, }, /* order=4 */
+ { LocalInterp_ULagrange_1dcube_50, }, /* order=5 */
+ { LocalInterp_ULagrange_1dcube_60, }, /* order=6 */
+ },
+ },
+
+ {
+ /* N_dims = 2 */
{
- /* 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) */
- },
+ /* molecule_family = molecule_family_cube */
+ { NULL_INTERP_FN_PTR , }, /* order=0 (unused) */
+ { LocalInterp_ULagrange_2dcube_10, }, /* order=1 */
+ { LocalInterp_ULagrange_2dcube_20, }, /* order=2 */
+ { LocalInterp_ULagrange_2dcube_30, }, /* order=3 */
+ { LocalInterp_ULagrange_2dcube_40, }, /* order=4 */
+ { NULL_INTERP_FN_PTR , }, /* order=5 (n.i.) */
+ { NULL_INTERP_FN_PTR , }, /* order=6 (n.i.) */
},
+ },
+ {
+ /* N_dims = 3 */
{
- /* N_dims = 1 */
- {
- /* molecule_family = molecule_family_cube */
- { NULL_INTERP_FN_PTR , }, /* order=0 (unused) */
- { LocalInterp_ILU_1d_cube_o1_s0, }, /* order=1 */
- { LocalInterp_ILU_1d_cube_o2_s0, }, /* order=2 */
- { LocalInterp_ILU_1d_cube_o3_s0, }, /* order=3 */
- { LocalInterp_ILU_1d_cube_o4_s0, }, /* order=4 */
- { LocalInterp_ILU_1d_cube_o5_s0, }, /* order=5 */
- { LocalInterp_ILU_1d_cube_o6_s0, }, /* order=6 */
- },
+ /* molecule_family = molecule_family_cube */
+ { NULL_INTERP_FN_PTR , }, /* order=0 (unused) */
+ { LocalInterp_ULagrange_3dcube_10, }, /* order=1 */
+ { LocalInterp_ULagrange_3dcube_20, }, /* order=2 */
+ { LocalInterp_ULagrange_3dcube_30, }, /* order=3 */
+ { LocalInterp_ULagrange_3dcube_40, }, /* order=4 */
+ { NULL_INTERP_FN_PTR , }, /* order=5 (n.i.) */
+ { NULL_INTERP_FN_PTR , }, /* order=6 (n.i.) */
},
+ },
+ /* end of interp_operator == interp_operator_Lagrange */
+ },
+
+ {
+ /* interp_operator == interp_operator_Hermite */
+ {
+ /* N_dims = 0 ==> unused */
{
- /* N_dims = 2 */
- {
- /* molecule_family = molecule_family_cube */
- { NULL_INTERP_FN_PTR , }, /* order=0 (unused) */
- { LocalInterp_ILU_2d_cube_o1_s0, }, /* order=1 */
- { LocalInterp_ILU_2d_cube_o2_s0, }, /* order=2 */
- { LocalInterp_ILU_2d_cube_o3_s0, }, /* order=3 */
- { LocalInterp_ILU_2d_cube_o4_s0, }, /* order=4 */
- { NULL_INTERP_FN_PTR , }, /* order=5 (n.i.) */
- { NULL_INTERP_FN_PTR , }, /* order=6 (n.i.) */
- },
+ /* 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 */
{
- /* N_dims = 3 */
- {
- /* molecule_family = molecule_family_cube */
- { NULL_INTERP_FN_PTR , }, /* order=0 (unused) */
- { LocalInterp_ILU_3d_cube_o1_s0, }, /* order=1 */
- { LocalInterp_ILU_3d_cube_o2_s0, }, /* order=2 */
- { LocalInterp_ILU_3d_cube_o3_s0, }, /* order=3 */
- { LocalInterp_ILU_3d_cube_o4_s0, }, /* order=4 */
- { NULL_INTERP_FN_PTR , }, /* order=5 (n.i.) */
- { NULL_INTERP_FN_PTR , }, /* order=6 (n.i.) */
- },
+ /* molecule_family = molecule_family_cube */
+ { NULL_INTERP_FN_PTR , }, /* order=0 (unused) */
+ { NULL_INTERP_FN_PTR , }, /* order=1 (n.i.) */
+ { NULL_INTERP_FN_PTR , }, /* order=2 (n.i.) */
+ { NULL_INTERP_FN_PTR , }, /* order=3 (n.i.) */
+ { NULL_INTERP_FN_PTR , }, /* order=4 (n.i.) */
+ { NULL_INTERP_FN_PTR , }, /* order=5 (n.i.) */
+ { NULL_INTERP_FN_PTR , }, /* order=6 (n.i.) */
},
- };
+ },
+
+ {
+ /* N_dims = 2 */
+ {
+ /* molecule_family = molecule_family_cube */
+ { NULL_INTERP_FN_PTR , }, /* order=0 (unused) */
+ { NULL_INTERP_FN_PTR , }, /* order=1 (n.i.) */
+ { NULL_INTERP_FN_PTR , }, /* order=2 (n.i.) */
+ { NULL_INTERP_FN_PTR , }, /* order=3 (n.i.) */
+ { NULL_INTERP_FN_PTR , }, /* order=4 (n.i.) */
+ { NULL_INTERP_FN_PTR , }, /* order=5 (n.i.) */
+ { NULL_INTERP_FN_PTR , }, /* order=6 (n.i.) */
+ },
+ },
+
+ {
+ /* N_dims = 3 */
+ {
+ /* molecule_family = molecule_family_cube */
+ { NULL_INTERP_FN_PTR , }, /* order=0 (unused) */
+ { NULL_INTERP_FN_PTR , }, /* order=1 (n.i.) */
+ { NULL_INTERP_FN_PTR , }, /* order=2 (n.i.) */
+ { NULL_INTERP_FN_PTR , }, /* order=3 (n.i.) */
+ { NULL_INTERP_FN_PTR , }, /* order=4 (n.i.) */
+ { NULL_INTERP_FN_PTR , }, /* order=5 (n.i.) */
+ { NULL_INTERP_FN_PTR , }, /* order=6 (n.i.) */
+ },
+ },
+
+ /* end of interp_operator == interp_operator_Hermite */
+ }
+
+ };
/* look up the subfunction to do the interpolation */
-p_interp_fn_t p_interp_fn = p_interp_fn_table[N_dims]
+p_interp_fn_t p_interp_fn = p_interp_fn_table[interp_operator]
+ [N_dims]
[molecule_family]
[order]
[smoothing];
@@ -1290,11 +1537,12 @@ if (p_interp_fn == NULL_INTERP_FN_PTR)
CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL,
__LINE__, __FILE__, CCTK_THORNSTRING,
"\n"
-" CCTK_InterpLocalUniform(): interpolation not implemented\n"
-" for the combination N_dims=%d,\n"
-" molecule_family=\"%s\",\n"
-" order=%d, smoothing=%d",
- N_dims, molecule_family_string, order, smoothing);
+" CCTK_InterpLocalUniform():\n"
+" interpolation not implemented for the combination\n"
+" interp_operator=\"%s\", N_dims=%d\n"
+" molecule_family=\"%s\", order=%d, smoothing=%d",
+ interp_operator_string, N_dims,
+ molecule_family_string, order, smoothing);
return UTIL_ERROR_BAD_INPUT; /*** ERROR RETURN ***/
}
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-03 00:07:19 +0000