diff options
Diffstat (limited to 'src/template.c')
-rw-r--r-- | src/template.c | 1676 |
1 files changed, 1676 insertions, 0 deletions
diff --git a/src/template.c b/src/template.c new file mode 100644 index 0000000..fe16818 --- /dev/null +++ b/src/template.c @@ -0,0 +1,1676 @@ +/*@@ + @file template.c + @date 23.Oct.2001 + @author Jonathan Thornburg <jthorn@aei.mpg.de> + @desc + This file is a template for generalized interpolation for + 1-, 2-, or 3-d molecules. It's used by defining various + C preprocessor macros, then #including this file. Each + such inclusion defines a single interpolation function + (plus other supporting functions local to this file), + which does generalized interpolation for a single molecule, + or more accurately for a single molecule for each possible + operation_codes[] value. + @enddesc + + @hdate 29.Aug.2002 + @hauthor Jonathan Thornburg <jthorn@aei.mpg.de> + @hdesc Restructure code and break it up into smaller functions + to reduce excessive compilation cpu/memory with big + 3-D molecules. + @endhdesc + + @hdate 28.Jan.2003 + @hauthor Jonathan Thornburg <jthorn@aei.mpg.de> + @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: + <math.h> + <limits.h> + <stdlib.h> + <string.h> + <stdio.h> + "util_ErrorCodes.h" + "cctk.h" + "InterpLocalUniform.h" + + The following preprocessor macros must be defined before + #including this file (note the macros which are file names + must all include the proper quotes for a #include in their + expansion, e.g. + #define DATA_VAR_DCL_FILE_NAME \ + "coeffs/2d.cube.size3/data-var.dcl.c" + @enddesc + + + @var FUNCTION_NAME + @vdesc The name of the interpolation function, e.g. + #define FUNCTION_NAME AEILocalInterp_U_LagTP_2cube_20 + @endvar + + @var N_DIMS + @vdesc The number of dimensions in which to interpolate, e.g. + #define N_DIMS 2 + The present implementation restricts this to 1, 2, + or 3, but this could easily be changed if needed. + Note that MAX_N_DIMS (defined in "InterpLocalUniform.c") + is a compile-time upper bound for N_DIMS, useful for + sizing arrays etc. + @endvar + + @var XYZ + @desc A comma-separated list of the (x,y,z) coordinates in the + interpolation, e.g. + #define XYZ x,y + @endvar + + @var FP_XYZ + @desc A comma-separated list of function-prototype declarations + of the (x,y,z) coordinates in the interpolation, e.g. + #define FP_XYZ fp x, fp y + @endvar + + @var STRIDE_IJK + @desc A comma-separated list of the (i,j,k) array strides for the + interpolation, e.g. + #define STRIDE_IJK stride_i, stride_j + @endvar + + @var JACOBIAN_MIJK_STRIDE + @desc A comma-separated list of the (i,j,k) strides for the + Jacobian, e.g. + #define JACOBIAN_MIJK_STRIDE \ + Jacobian_mi_stride, Jacobian_mj_stride + @endvar + + @var MOLECULE_MIN_M + @vdesc The minimum m coordinate in the molecule, e.g. + #define MOLECULE_MIN_M -1 + The present implementation takes this to be the same in each + dimension, but this could easily be changed if needed. + @endvar + + @var MOLECULE_MAX_M + @vdesc The maximum m coordinate in the molecule, e.g. + #define MOLECULE_MAX_M 1 + The present implementation takes this to be the same in each + dimension, but this could easily be changed if needed. + @endvar + + @var MOLECULE_SIZE + @vdesc The diameter of (number of points in) the molecules to be used, + e.g. + #define MOLECULE_SIZE 3 + The present implementation takes this to be the same in each + dimension, but this could easily be changed if needed. + @endvar + + @var HAVE_OP_{I,DX,DY,DXX,DXY,DYY,...} + @vdesc Each of these symbols should be defined or not defined + according as if the corresponding derivative operator is + to be supported or not supported by this function, e.g. + #define HAVE_OP_I + #define HAVE_OP_DX + #define HAVE_OP_DY + #define HAVE_OP_DXX + #define HAVE_OP_DXY + #define HAVE_OP_DYY + if we support all 1st and 2nd derivatives in 2-D. + @endvar + + @var DATA_STRUCT + @vdesc The name of a C struct containing all the data variables, + e.g. + #define DATA_STRUCT data_struct_2d_cube_size3 + @endvar + + @var COEFF_STRUCT + @vdesc The name of a C struct containing all the molecule coefficients, + e.g. + #define COEFFS_STRUCT coeffs_struct_2d_cube_size3 + @endvar + + @var LOAD_DATA_{REAL,REAL{4,8,16},COMPLEX,COMPLEX{8,16,32}} + @vdesc The name of a C function to load data from the input arrays + into the DATA_STRUCT structure. Typically these will be + chosem from among the functions defined in + "common/load-template.[ch]", e.g. + #define LOAD_DATA_REAL AEILocalInterp_load_2dcube3_r + #define LOAD_DATA_REAL4 AEILocalInterp_load_2dcube3_r4 + #define LOAD_DATA_REAL8 AEILocalInterp_load_2dcube3_r8 + #define LOAD_DATA_REAL16 AEILocalInterp_load_2dcube3_r16 + #define LOAD_DATA_COMPLEX AEILocalInterp_load_2dcube3_c + #define LOAD_DATA_COMPLEX8 AEILocalInterp_load_2dcube3_c8 + #define LOAD_DATA_COMPLEX16 AEILocalInterp_load_2dcube3_c16 + #define LOAD_DATA_COMPLEX32 AEILocalInterp_load_2dcube3_c32 + @endvar + + @var EVALUATE_MOLECULE + @vdesc The name of a C function to evaluate the dot product + of the DATA_STRUCT and COEFF_STRUCT structures. Typically + these will be chosem from among the functions defined in + "common/evaluate.[ch]", e.g. + #define EVALUATE_MOLECULE LocalInterp_eval_2d_cube6 + @endvar + + @var STORE_COEFFS + @vdesc The name of a C function to store a COEDFF_STRUCT structure + of the DATA_STRUCT and COEFF_STRUCT structures. Typically + this will be chosem from among the functions defined in + "common/store.[ch]", e.g. + #define STORE_COEFFS LocalInterp_store_2d_cube6 + @endvar + + @var COEFFS_{I,DX,DY,DXX,DXY,DYY,...}_COMPUTE_FILE_NAME + @vdesc Each of these macros should be the name of a file + (presumably machine-generated) containing a sequence + of C assignment statements to compute all the coefficients + in a corresponding coeffs_* structure as polynomials + in the variables (x,y,z), e.g. (for 2D size-3 molecules, + dx operator) + fp t36; + fp t42; + fp t41; + fp t40; + fp t39; + fp t37; + t36 = RATIONAL(1.0,3.0)*x; + t42 = RATIONAL(1.0,4.0)*y+t36; + t41 = t36+RATIONAL(-1.0,4.0)*y; + t40 = RATIONAL(1.0,6.0); + t39 = RATIONAL(-1.0,6.0); + t37 = RATIONAL(-2.0,3.0)*x; + coeffs_dx->coeff_m1_m1 = t39+t42; + coeffs_dx->coeff_0_m1 = t37; + coeffs_dx->coeff_p1_m1 = t40+t41; + coeffs_dx->coeff_m1_0 = t36+t39; + coeffs_dx->coeff_0_0 = t37; + coeffs_dx->coeff_p1_0 = t36+t40; + coeffs_dx->coeff_m1_p1 = t39+t41; + coeffs_dx->coeff_0_p1 = t37; + coeffs_dx->coeff_p1_p1 = t40+t42; + + As illustrated, the code may use the macro RATIONAL + (defined later in this file) to represent rational-number + coefficients, and it may also declare temporary variables + as needed. + + Note that this is the only included code which depends + on the actual interpolation scheme used; all the rest + just depends on the interpolation dimension and molecule + family and size. + @endvar + + @version $Header$ + @@*/ + +/******************************************************************************/ + +/* + * ***** prototypes for functions local to this file ***** + */ + +#ifdef HAVE_OP_I + static + void compute_coeffs_I(FP_XYZ, struct COEFFS_STRUCT *coeffs_I); +#endif +#ifdef HAVE_OP_DX + static + void compute_coeffs_dx(FP_XYZ, struct COEFFS_STRUCT *coeffs_dx); +#endif +#ifdef HAVE_OP_DY + static + void compute_coeffs_dy(FP_XYZ, struct COEFFS_STRUCT *coeffs_dy); +#endif +#ifdef HAVE_OP_DZ + static + void compute_coeffs_dz(FP_XYZ, struct COEFFS_STRUCT *coeffs_dz); +#endif +#ifdef HAVE_OP_DXX + static + void compute_coeffs_dxx(FP_XYZ, struct COEFFS_STRUCT *coeffs_dxx); +#endif +#ifdef HAVE_OP_DXY + static + void compute_coeffs_dxy(FP_XYZ, struct COEFFS_STRUCT *coeffs_dxy); +#endif +#ifdef HAVE_OP_DXZ + static + void compute_coeffs_dxz(FP_XYZ, struct COEFFS_STRUCT *coeffs_dxz); +#endif +#ifdef HAVE_OP_DYY + static + void compute_coeffs_dyy(FP_XYZ, struct COEFFS_STRUCT *coeffs_dyy); +#endif +#ifdef HAVE_OP_DYZ + static + void compute_coeffs_dyz(FP_XYZ, struct COEFFS_STRUCT *coeffs_dyz); +#endif +#ifdef HAVE_OP_DZZ + static + void compute_coeffs_dzz(FP_XYZ, struct COEFFS_STRUCT *coeffs_dzz); +#endif + +/******************************************************************************/ + +/*@@ + @routine FUNCTION_NAME + @date 23.Oct.2001 + @author Jonathan Thornburg <jthorn@aei.mpg.de> + @desc + This function does generalized interpolation of one or more + 2d arrays to arbitrary points. For details, see the header + comments for InterpLocalUniform() (in "InterpLocalUniform.c" + in this same directory). + + This function's arguments are mostly all a subset of those of + InterpLocalUniform() ; the difference is that this function + takes all its arguments explicitly, whereas InputLocalArrays() + takes some of them indirectly via a key/value parameter table. + Here we document only the "new" explicit arguments, and these + 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 + point in the pointed-to structure. + @vtype struct error_flags* error_flags; + @vio out + + @var molecule_structure_flags + @vdesc If this pointer is non-NULL, we store flags describing + the interpolation molecule's structure in the pointed-to + structure. + @vtype struct molecule_structure_flags* molecule_structure_flags; + @vio out + + @var molecule_min_max_m_info + @vdesc If this pointer is non-NULL, we store the interpolation + molecule's min/max m in the pointed-to structure. + @vtype struct molecule_min_max_m_info* molecule_min_max_m_info; + @vio out + + @var molecule_positions + @vdesc If this pointer is non-NULL, we store the interpolation + molecule's positions in the (caller-supplied) arrays + pointed to by this pointer. + @vtype CCTK_INT* const molecule_positions[]; + @vio out + + @var Jacobian_info + @vdesc If this pointer is non-NULL, we store the Jacobian of + the interpolation in the arrays (and in the manner) + pointed to by this structure. + @vtype struct Jacobian_info* Jacobian_info; + @vio out + + @returntype int + @returndesc This function's return results are a subset of those of + InterpLocalUniform(): + 0 success + UTIL_ERROR_BAD_INPUT one of the input arguments is invalid + CCTK_ERROR_INTERP_POINT_OUTSIDE + interpolation point is out of range + @endreturndesc + + @@*/ +int FUNCTION_NAME(/***** 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[], + 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[], + const CCTK_INT input_array_strides[], + const CCTK_INT input_array_min_subscripts[], + const CCTK_INT input_array_max_subscripts[], + 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[], + /***** operation info */ + const CCTK_INT operand_indices[], + const CCTK_INT operation_codes[], + /***** other return results *****/ + struct error_flags* error_flags, + struct molecule_structure_flags* molecule_structure_flags, + struct molecule_min_max_m_info* molecule_min_max_m_info, + CCTK_INT* const molecule_positions[], + struct Jacobian_info* Jacobian_info) +{ +/* + * ***** Naming conventions ***** + * + * in, out = 0-origin indices each selecting an input/output array + * pt = 0-origin index selecting an interpolation point + * part = 0-origin index selecting real/imaginary part of a complex number + */ + +/* + * ***** Implementation notes: ***** + * + * The basic outline of this function is as follows: + * + * store molecule structure flags + * if (querying molecule min/max m) + * then store molecule min/max m + * compute "which derivatives are wanted" flags + * precompute 1/dx factors + * set up input array strides + * set up Jacobian strides (or dummy values if no Jacobian info) + * for (int pt = 0 ; pt < N_interp_point ; ++pt) + * { + * struct DATA_STRUCT data; + * struct COEFFS_STRUCT coeffs_I, coeffs_dx, ..., coeffs_dzz; + * for (int axis = 0 ; axis < N_dims ; ++axis) + * { + * switch (interp_coords_type_codes) + * { + * case CCTK_VARIABLE_REAL: + * load interp coords for this axis from + * const CCTK_REAL *const interp_coords[] + * break; + * ... + * } + * } + * compute position of interpolation molecules + * with respect to the grid into the XYZ variables + * if (querying molecule positions) + * then store this molecule position + * compute molecule center 1-d position in input arrays + * + * /# compute interpolation coefficients at this point #/ + * if (want_I) + * then compute_coeffs_I(XYZ, &coeffs_I); + * if (want_dx) + * then compute_coeffs_dx(XYZ, &coeffs_dx); + * ... + * if (want_dzz) + * then compute_coeffs_dzz(XYZ, &coeffs_dx); + * + * for (int out = 0 ; out < N_output_arrays ; ++out) + * { + * const int in = operand_indices[out]; + * ***decode*** the output array datatype + * to determine whether it's real or complex + * const int N_output_parts = output array is complex ? 2 : 1; + * for (int part = 0 ; part < N_output_parts ; ++part) + * { + * if ( (input_arrays[in] != NULL) + * && ( (input_arrays[in] != value at last load) + * || (part != value at last load) ) ) + * then { + * save input_arrays[in] and part for + * "previous value" test above + * ***decode*** the input array datatype + * to determine whether it's real or complex + * const int N_input_parts + * = input array is complex ? 2 : 1; + * if (N_input_parts != N_output_parts) + * then error(...) + * switch (input_array_type_codes[in]) + * { + * case CCTK_VARIABLE_REAL: + * LOAD_DATA_REAL(input_array_ptr, + * STRIDE_IJK, + * &data); + * break; + * ... + * case CCTK_VARIABLE_COMPLEX: + * LOAD_DATA_COMPLEX(input_array_ptr, + * STRIDE_IJK, part, + * &data); + * break; + * ... + * } + * } + * + * if (output_arrays[out] != NULL) + * then { + * /# interpolate at this point #/ + * fp result; + * switch (operation_codes[out]) + * { + * case 0: + * result = EVALUATE_MOLECULE(&coeffs_I, + * &data); + * break; + * case 1: + * result = inverse_dx + * * EVALUATE_MOLECULE(&coeffs_dx, + * &data); + * break; + * ... + * case 33: + * result = inverse_dz * inverse_dz + * * EVALUATE_MOLECULE(&coeffs_dzz, + * &data); + * break; + * } + * + * /# store result in output array #/ + * switch (output_array_type_codes[out]) + * { + * case CCTK_VARIABLE_REAL: + * store result in + * CCTK_REAL *const output_arrays[] + * break; + * ... + * case CCTK_VARIABLE_COMPLEX: + * store result in + * CCTK_COMPLEX *const output_arrays[][part] + * break; + * ... + * } + * } + * if (querying Jacobian && (Jacobian_pointer[out] != NULL)) + * then { + * /# store Jacobian coefficients at this point #/ + * CCTK_REAL *const Jacobian_ptr + * = Jacobian_info->Jacobian_pointer[out] + * + Jacobian_info->Jacobian_offset[out]; + * switch (operation_code) + * { + * case 0: + * STORE_COEFFS(1.0, &coeffs_I, + * Jacobian_ptr, JACOBIAN_MIJK_STRIDE, + * pt, part, + * break; + * case 1: + * STORE_COEFFS(inverse_dx, &coeffs_dx, + * Jacobian_ptr, JACOBIAN_MIJK_STRIDE, + * pt, part, + * break; + * ... + * case 33: + * STORE_COEFFS(inverse_dz*inverse_dz, &coeffs_dzz, + * Jacobian_ptr, JACOBIAN_MIJK_STRIDE, + * pt, part, + * break; + * } + * } + * } + * } + * } + * + * For complex datatypes we pointer-alias the N-dimensional input array + * to a 2-dimensional array where the 1st axis is the 1-D subscript + * corresponding to the N input axes, and the 2nd axes has 2 elements + * subscripted by [part] for the (real,imaginary) components of the + * complex values. + * + * We do all floating-point computations in type "fp" (typically a typedef + * for CCTK_REAL), so arrays of higher precision than this will incur extra + * rounding errors. In practice these should be negligible compared to the + * interpolation errors. + */ + + +/* layout of axes in N_dims-element arrays */ +#define X_AXIS 0 +#define Y_AXIS 1 +#define Z_AXIS 2 +#if (N_DIMS > 3) + #error "N_DIMS may not be > 3!" +#endif + +/* layout of axes and min/max ends in out_of_range_tolerance[] array */ +#define X_AXIS_MIN 0 +#define X_AXIS_MAX 1 +#define Y_AXIS_MIN 2 +#define Y_AXIS_MAX 3 +#define Z_AXIS_MIN 4 +#define Z_AXIS_MAX 5 +#if (N_DIMS > 3) + #error "N_DIMS may not be > 3!" +#endif + +/* basic sanity check on molecule size */ +#define MOLECULE_M_COUNT (MOLECULE_MAX_M - MOLECULE_MIN_M + 1) +#if (MOLECULE_SIZE != MOLECULE_M_COUNT) + #error "MOLECULE_SIZE inconsistent with MOLECULE_{MIN,MAX}_M!" +#endif + + +/* macros used by machine-generated interpolation coefficient expressions */ +/* + * FIXME: right now this is used as (eg) RATIONAL(1.0,2.0); + * it might be cleaner if it were RATIONAL(1,2) with the + * preprocessor ## operator used to glue on the .0 + * (I _think_ this is portable, but is it really?) + */ +#define RATIONAL(num,den) (num/den) + + +/* + * store molecule structure flags, molecule min/max m (if requested) + */ +if (molecule_structure_flags != NULL) + then { + molecule_structure_flags->MSS_is_fn_of_interp_coords = 0; + molecule_structure_flags->MSS_is_fn_of_which_operation = 0; + molecule_structure_flags->MSS_is_fn_of_input_array_values = 0; + molecule_structure_flags->Jacobian_is_fn_of_input_array_values = 0; + } +if (molecule_min_max_m_info != NULL) + then { + int axis; + for (axis = 0 ; axis < N_DIMS ; ++axis) + { + molecule_min_max_m_info->molecule_min_m[axis] = MOLECULE_MIN_M; + molecule_min_max_m_info->molecule_max_m[axis] = MOLECULE_MAX_M; + } + } + + +/* + * compute "which derivatives are wanted" flags + */ + { +#ifdef HAVE_OP_I + bool want_I = false; +#endif +#ifdef HAVE_OP_DX + bool want_dx = false; +#endif +#ifdef HAVE_OP_DY + bool want_dy = false; +#endif +#ifdef HAVE_OP_DZ + bool want_dz = false; +#endif +#ifdef HAVE_OP_DXX + bool want_dxx = false; +#endif +#ifdef HAVE_OP_DXY + bool want_dxy = false; +#endif +#ifdef HAVE_OP_DXZ + bool want_dxz = false; +#endif +#ifdef HAVE_OP_DYY + bool want_dyy = false; +#endif +#ifdef HAVE_OP_DYZ + bool want_dyz = false; +#endif +#ifdef HAVE_OP_DZZ + bool want_dzz = false; +#endif + + { +int out; + for (out = 0 ; out < N_output_arrays ; ++out) + { + switch (operation_codes[out]) + { + #ifdef HAVE_OP_I + case 0: want_I = true; break; + #endif + #ifdef HAVE_OP_DX + case 1: want_dx = true; break; + #endif + #ifdef HAVE_OP_DY + case 2: want_dy = true; break; + #endif + #ifdef HAVE_OP_DZ + case 3: want_dz = true; break; + #endif + #ifdef HAVE_OP_DXX + case 11: want_dxx = true; break; + #endif + #ifdef HAVE_OP_DXY + case 12: + case 21: want_dxy = true; break; + #endif + #ifdef HAVE_OP_DXZ + case 13: + case 31: want_dxz = true; break; + #endif + #ifdef HAVE_OP_DYY + case 22: want_dyy = true; break; + #endif + #ifdef HAVE_OP_DYZ + case 23: + case 32: want_dyz = true; break; + #endif + #ifdef HAVE_OP_DZZ + case 33: want_dzz = true; break; + #endif + default: + CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, +"\n" +" CCTK_InterpLocalUniform(): generalized interpolation operation_code %d\n" +" is not supported!\n" +" (0-origin) output #out=%d" + , + (int) operation_codes[out], + out); /*NOTREACHED*/ + return UTIL_ERROR_BAD_INPUT; /*** ERROR RETURN ***/ + } + } + } + + +/* + * save origin/delta variables, precompute 1/delta factors + * ... in theory we could compute only those factors we're going to use, + * but it's not worth the trouble, so we just compute them all + */ + { +#if N_DIMS >= 1 + #if defined(HAVE_OP_DX) \ + || defined(HAVE_OP_DXY) || defined(HAVE_OP_DXZ) \ + || defined(HAVE_OP_DXX) + const fp inverse_dx = 1.0 / coord_delta[X_AXIS]; + #endif +#endif +#if N_DIMS >= 2 + #if defined(HAVE_OP_DY) \ + || defined(HAVE_OP_DXY) || defined(HAVE_OP_DYZ) \ + || defined(HAVE_OP_DYY) + const fp inverse_dy = 1.0 / coord_delta[Y_AXIS]; + #endif +#endif +#if N_DIMS >= 3 + #if defined(HAVE_OP_DZ) \ + || defined(HAVE_OP_DXZ) || defined(HAVE_OP_DYZ) \ + || defined(HAVE_OP_DZZ) + const fp inverse_dz = 1.0 / coord_delta[Z_AXIS]; + #endif +#endif +#if (N_DIMS > 3) + #error "N_DIMS may not be > 3!" +#endif + + +/* + * set up input array strides + */ +#if (N_DIMS >= 1) + const int stride_i = input_array_strides[X_AXIS]; +#endif +#if (N_DIMS >= 2) + const int stride_j = input_array_strides[Y_AXIS]; +#endif +#if (N_DIMS >= 3) + const int stride_k = input_array_strides[Z_AXIS]; +#endif +#if (N_DIMS > 3) + #error "N_DIMS may not be > 3!" +#endif + + +/* + * set up Jacobian m strides, or dummy values if no Jacobian info + */ +#if N_DIMS >= 1 + const int Jacobian_mi_stride = (Jacobian_info == NULL) + ? 0 + : Jacobian_info->Jacobian_m_strides[X_AXIS]; +#endif +#if (N_DIMS >= 2) + const int Jacobian_mj_stride = (Jacobian_info == NULL) + ? 0 + : Jacobian_info->Jacobian_m_strides[Y_AXIS]; +#endif +#if (N_DIMS >= 3) + const int Jacobian_mk_stride = (Jacobian_info == NULL) + ? 0 + : Jacobian_info->Jacobian_m_strides[Z_AXIS]; +#endif +#if (N_DIMS > 3) + #error "N_DIMS may not be > 3!" +#endif + +/* + * interpolate at each point + */ +int pt; + for (pt = 0 ; pt < N_interp_points ; ++pt) + { + /* this struct holds a molecule-sized piece of a single */ + /* real input array, or of a real/complex component of */ + /* a complex input array */ + struct DATA_STRUCT data; + + /* each of these structs holds a molecule-sized set of */ + /* interpolation coefficients -- recall we are representing */ + /* the interpolant as a linear combination of the data values */ + #ifdef HAVE_OP_I + struct COEFFS_STRUCT coeffs_I; + #endif + #ifdef HAVE_OP_DX + struct COEFFS_STRUCT coeffs_dx; + #endif + #ifdef HAVE_OP_DY + struct COEFFS_STRUCT coeffs_dy; + #endif + #ifdef HAVE_OP_DZ + struct COEFFS_STRUCT coeffs_dz; + #endif + #ifdef HAVE_OP_DXX + struct COEFFS_STRUCT coeffs_dxx; + #endif + #ifdef HAVE_OP_DXY + struct COEFFS_STRUCT coeffs_dxy; + #endif + #ifdef HAVE_OP_DXZ + struct COEFFS_STRUCT coeffs_dxz; + #endif + #ifdef HAVE_OP_DYY + struct COEFFS_STRUCT coeffs_dyy; + #endif + #ifdef HAVE_OP_DYZ + struct COEFFS_STRUCT coeffs_dyz; + #endif + #ifdef HAVE_OP_DZZ + struct COEFFS_STRUCT coeffs_dzz; + #endif + + + /* + * load the interpolation point coordinates + * from the interp_coords[] arrays + * FIXME: Maybe it would be better (faster) to do this + * with N_DIMS open-coded calls on a function? + * But then we'd have to have a sentinal value + * return for the unknown-type-code error case. + * Yuk! :( :( + */ + fp interp_coords_fp[N_DIMS]; + { + int axis; + for (axis = 0 ; axis < N_DIMS ; ++axis) + { + /* pointer to array of interp coords for this axis */ + const void *const interp_coords_ptr = interp_coords[axis]; + + switch (interp_coords_type_code) + { + case CCTK_VARIABLE_REAL: + { + const CCTK_REAL *const interp_coords_ptr_real + = (const CCTK_REAL *) interp_coords_ptr; + interp_coords_fp[axis] = interp_coords_ptr_real[pt]; + break; + } + + #ifdef HAVE_CCTK_REAL4 + case CCTK_VARIABLE_REAL4: + { + const CCTK_REAL4 *const interp_coords_ptr_real4 + = (const CCTK_REAL4 *) interp_coords_ptr; + interp_coords_fp[axis] = interp_coords_ptr_real4[pt]; + break; + } + #endif + + #ifdef HAVE_CCTK_REAL8 + case CCTK_VARIABLE_REAL8: + { + const CCTK_REAL8 *const interp_coords_ptr_real8 + = (const CCTK_REAL8 *) interp_coords_ptr; + interp_coords_fp[axis] = interp_coords_ptr_real8[pt]; + break; + } + #endif + + #ifdef HAVE_CCTK_REAL16 + case CCTK_VARIABLE_REAL16: + { + /* FIXME: maybe we should warn (once per cactus run) */ + /* that we may be doing arithmetic in lower */ + /* precision than the interp coords? */ + const CCTK_REAL16 *const interp_coords_ptr_real16 + = (const CCTK_REAL16 *) interp_coords_ptr; + interp_coords_fp[axis] + = interp_coords_ptr_real16[pt]; + break; + } + #endif + + default: + CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, +"\n" +" CCTK_InterpLocalUniform(): interpolation coordinates datatype %d\n" +" is not supported!" + , + interp_coords_type_code); + /*NOTREACHED*/ + return UTIL_ERROR_BAD_INPUT; /*** ERROR RETURN ***/ + /* end of switch (interp_coords_type_code) */ + } + + /* end of for (axis = ...) loop */ + } + } + +#ifdef PICKLE + #if (N_DIMS == 1) +printf("pt=%d interp coord %.16g\n", + pt, (double) interp_coords_fp[0]); + #endif + #if (N_DIMS == 3) +printf("pt=%d interp coords %.16g %.16g %.16g\n", + pt, (double) interp_coords_fp[0], + (double) interp_coords_fp[1], + (double) interp_coords_fp[2]); + #endif +#endif + + + /* + * compute position of interpolation molecules with respect to + * the grid, i.e. compute (x,y,z) coordinates of interpolation + * point relative to molecule center, in units of the grid spacing + * + * 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 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 + */ + { + 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 = AEILocalInterp_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_direction = -1; + } + else { + error_flags->error_ibndry = ibndry_max; + error_flags->error_direction = +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 + +#ifdef PICKLE + #if (N_DIMS == 1) +printf("interp center_i = %d\n", center_i); +printf("interp grid-relative x = %.16g\n", (double) x); + #endif + #if (N_DIMS == 3) +printf("interp center_ijk = %d %d %d\n", center_i, center_j, center_k); +printf("interp grid-relative xyz = %.16g %.16g %.16g\n", + (double) x, (double) y, (double) z); + #endif +#endif + + /* + * compute 1-d position of molecule center in input arrays + */ + { + #if (N_DIMS == 1) + const int molecule_center_posn = stride_i*center_i; + #elif (N_DIMS == 2) + const int molecule_center_posn = stride_i*center_i + + stride_j*center_j; + #elif (N_DIMS == 3) + const int molecule_center_posn = stride_i*center_i + + stride_j*center_j + + stride_k*center_k; + #else + #error "N_DIMS may not be > 3!" + #endif + + + /* + * molecule position queries + */ + if (molecule_positions != NULL) + then { + #if (N_DIMS >= 1) + molecule_positions[X_AXIS][pt] = center_i; + #endif + #if (N_DIMS >= 2) + molecule_positions[Y_AXIS][pt] = center_j; + #endif + #if (N_DIMS >= 3) + molecule_positions[Z_AXIS][pt] = center_k; + #endif + } + + + /* + * compute the coefficients at this point for whichever + * operation_codes[] values (derivatatives) are wanted + * (these are polynomials in the variables (x,y,z)) + */ + #ifdef HAVE_OP_I + if (want_I) + then compute_coeffs_I(XYZ, &coeffs_I); + #endif + #ifdef HAVE_OP_DX + if (want_dx) + then compute_coeffs_dx(XYZ, &coeffs_dx); + #endif + #ifdef HAVE_OP_DY + if (want_dy) + then compute_coeffs_dy(XYZ, &coeffs_dy); + #endif + #ifdef HAVE_OP_DZ + if (want_dz) + then compute_coeffs_dz(XYZ, &coeffs_dz); + #endif + #ifdef HAVE_OP_DXX + if (want_dxx) + then compute_coeffs_dxx(XYZ, &coeffs_dxx); + #endif + #ifdef HAVE_OP_DXY + if (want_dxy) + then compute_coeffs_dxy(XYZ, &coeffs_dxy); + #endif + #ifdef HAVE_OP_DXZ + if (want_dxz) + then compute_coeffs_dxz(XYZ, &coeffs_dxz); + #endif + #ifdef HAVE_OP_DYY + if (want_dyy) + then compute_coeffs_dyy(XYZ, &coeffs_dyy); + #endif + #ifdef HAVE_OP_DYZ + if (want_dyz) + then compute_coeffs_dyz(XYZ, &coeffs_dyz); + #endif + #ifdef HAVE_OP_DZZ + if (want_dzz) + then compute_coeffs_dzz(XYZ, &coeffs_dzz); + #endif + + + /* + * compute each output array at this point + */ + { + int out; + + /* + * next 2 initializers must be invalid values to make sure we + * execute the ***load*** the first time in the test at the + * top of the part loop below + */ + const void* input_array_ptr__last_load = NULL; + int part__last_load = -1; + + for (out = 0 ; out < N_output_arrays ; ++out) + { + const int in = operand_indices[out]; + const void* const input_array_ptr = input_arrays[in]; + + /* + * ***decode*** the output array datatype + * to determine whether it's real or complex, + */ + const int N_output_parts + = AEILocalInterp_decode_N_parts(output_array_type_codes[out]); + if (! ((N_output_parts == 1) || (N_output_parts == 2))) + then { +CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, +"\n" +" CCTK_InterpLocalUniform():\n" +" output array doesn't seem to be a real or complex number,\n" +" or more precisely, output array has number of \"real parts\"\n" +" (1=real, 2=complex) which isn't 1 or 2!\n" +" 0-origin output #out=%d\n" +" datatype code=%d\n" +" (datatype codes are defined by the Cactus flesh,\n" +" see src/include/cctk_Constants.h)\n" +" ==> N_parts=%d" + , + out, (int) output_array_type_codes[out], N_output_parts); + /*NOTREACHED*/ + return UTIL_ERROR_BAD_INPUT; /*** ERROR RETURN ***/ + } + + { + int part; + for (part = 0 ; part < N_output_parts ; ++part) + { + if ( (input_array_ptr != NULL) + && + ((input_array_ptr != input_array_ptr__last_load) + || (part != part__last_load)) ) + then { + /* remember when we did the following load */ + input_array_ptr__last_load = input_array_ptr; + part__last_load = part; + + /* + * ***decode*** the input array datatype + * to determine whether it's real or complex, + */ + { + const int N_input_parts + = AEILocalInterp_decode_N_parts(input_array_type_codes[in]); + if (N_input_parts != N_output_parts) + then { +CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, +"\n" +" CCTK_InterpLocalUniform():\n" +" data types are incompatible between input and output arrays, or\n" +" more precisely, number of \"real parts\" (1=real, 2=complex) differ!\n" +" 0-origin input #in =%d datatype code=%d N_parts=%d\n" +" 0-origin output #out=%d datatype code=%d N_parts=%d\n" +" (datatype codes are defined by the Cactus flesh,\n" +" see src/include/cctk_Constants.h)" + , + in, (int) input_array_type_codes[in], N_input_parts, + out, (int) output_array_type_codes[out], N_output_parts); + /*NOTREACHED*/ + return UTIL_ERROR_BAD_INPUT; + /*** ERROR RETURN ***/ + } + + /* + * load the molecule-sized piece of + * input_arrays[in][part] at this molecule + * position, into the data struct + */ + { + const int input_posn = molecule_center_posn + + input_array_offsets[in]; + switch (input_array_type_codes[in]) + { +case CCTK_VARIABLE_REAL: + { + const CCTK_REAL *const input_array_ptr_real + = ((const CCTK_REAL *) input_array_ptr) + input_posn; + LOAD_DATA_REAL(input_array_ptr_real, + STRIDE_IJK, + &data); + break; + } + +#ifdef HAVE_CCTK_REAL4 +case CCTK_VARIABLE_REAL4: + { + const CCTK_REAL4 *const input_array_ptr_real4 + = ((const CCTK_REAL4 *) input_array_ptr) + input_posn; + LOAD_DATA_REAL4(input_array_ptr_real4, + STRIDE_IJK, + &data); + break; + } +#endif + +#ifdef HAVE_CCTK_REAL8 +case CCTK_VARIABLE_REAL8: + { + const CCTK_REAL8 *const input_array_ptr_real8 + = ((const CCTK_REAL8 *) input_array_ptr) + input_posn; + LOAD_DATA_REAL8(input_array_ptr_real8, + STRIDE_IJK, + &data); + break; + } +#endif + +#ifdef HAVE_CCTK_REAL16 +case CCTK_VARIABLE_REAL16: + { + /* FIXME: maybe we should warn (once per cactus run) that we may be */ + /* doing arithmetic in lower precision than the data type? */ + const CCTK_REAL16 *const input_array_ptr_real16 + = ((const CCTK_REAL16 *) input_array_ptr) + input_posn; + LOAD_DATA_REAL16(input_array_ptr_real16, + STRIDE_IJK, + &data); + break; + } +#endif + +case CCTK_VARIABLE_COMPLEX: + { + const CCTK_REAL (*const input_array_ptr_complex)[COMPLEX_N_PARTS] + = ((const CCTK_REAL (*)[COMPLEX_N_PARTS]) input_array_ptr) + + input_posn; + LOAD_DATA_COMPLEX(input_array_ptr_complex, + STRIDE_IJK, part, + &data); + break; + } + +#ifdef HAVE_CCTK_COMPLEX8 +case CCTK_VARIABLE_COMPLEX8: + { + const CCTK_REAL4 (*const input_array_ptr_complex8)[COMPLEX_N_PARTS] + = ((const CCTK_REAL4 (*)[COMPLEX_N_PARTS]) input_array_ptr) + + input_posn; + LOAD_DATA_COMPLEX8(input_array_ptr_complex8, + STRIDE_IJK, part, + &data); + break; + } +#endif + +#ifdef HAVE_CCTK_COMPLEX16 +case CCTK_VARIABLE_COMPLEX16: + { + const CCTK_REAL8 (*const input_array_ptr_complex16)[COMPLEX_N_PARTS] + = ((const CCTK_REAL8 (*)[COMPLEX_N_PARTS]) input_array_ptr) + + input_posn; + LOAD_DATA_COMPLEX16(input_array_ptr_complex16, + STRIDE_IJK, part, + &data); + break; + } +#endif + + +#ifdef HAVE_CCTK_COMPLEX32 +case CCTK_VARIABLE_COMPLEX32: + { + /* FIXME: maybe we should warn (once per cactus run) that we may be */ + /* doing arithmetic in lower precision than the data type? */ + const CCTK_REAL16 (*const input_array_ptr_complex32)[COMPLEX_N_PARTS] + = ((const CCTK_REAL16 (*)[COMPLEX_N_PARTS]) input_array_ptr) + + input_posn; + LOAD_DATA_COMPLEX32(input_array_ptr_complex32, + STRIDE_IJK, part, + &data); + break; + } +#endif + +default: +CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, +"\n" +" CCTK_InterpLocalUniform(): input datatype %d not supported!\n" +" (0-origin) input #in=%d" + , + (int) input_array_type_codes[in], + in); /*NOTREACHED*/ +return UTIL_ERROR_BAD_INPUT; /*** ERROR RETURN ***/ + /* end of switch (input_array_type_codes[in]) */ + } + } + + } + /* end of load input array values */ + } + + + if (output_arrays[out] != NULL) + then { + /* + * interpolate at this point + */ + fp result; + + switch (operation_codes[out]) + { + #ifdef HAVE_OP_I + case 0: + result = EVALUATE_MOLECULE(&coeffs_I, + &data); + break; + #endif + #ifdef HAVE_OP_DX + case 1: + result = inverse_dx + * EVALUATE_MOLECULE(&coeffs_dx, + &data); + break; + #endif + #ifdef HAVE_OP_DY + case 2: + result = inverse_dy + * EVALUATE_MOLECULE(&coeffs_dy, + &data); + break; + #endif + #ifdef HAVE_OP_DZ + case 3: + result = inverse_dz + * EVALUATE_MOLECULE(&coeffs_dz, + &data); + break; + #endif + #ifdef HAVE_OP_DXX + case 11: + result = inverse_dx * inverse_dx + * EVALUATE_MOLECULE(&coeffs_dxx, + &data); + break; + #endif + #ifdef HAVE_OP_DXY + case 12: + case 21: + result = inverse_dx * inverse_dy + * EVALUATE_MOLECULE(&coeffs_dxy, + &data); + break; + #endif + #ifdef HAVE_OP_DXZ + case 13: + case 31: + result = inverse_dx * inverse_dz + * EVALUATE_MOLECULE(&coeffs_dxz, + &data); + break; + #endif + #ifdef HAVE_OP_DYY + case 22: + result = inverse_dy * inverse_dy + * EVALUATE_MOLECULE(&coeffs_dyy, + &data); + break; + #endif + #ifdef HAVE_OP_DYZ + case 23: + case 32: + result = inverse_dy * inverse_dz + * EVALUATE_MOLECULE(&coeffs_dyz, + &data); + break; + #endif + #ifdef HAVE_OP_DZZ + case 33: + result = inverse_dz * inverse_dz + * EVALUATE_MOLECULE(&coeffs_dzz, + &data); + break; + #endif + default: +CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, +"\n" +" CCTK_InterpLocalUniform(): generalized interpolation operation_code %d\n" +" is not supported!\n" +" (0-origin) output #out=%d" + , + (int) operation_codes[out], + out); + return UTIL_ERROR_BAD_INPUT; + /*** ERROR RETURN ***/ + /* end of switch (operation_codes[out]) */ + } + + + /* + * store result in output array + */ + switch (output_array_type_codes[out]) + { + +case CCTK_VARIABLE_REAL: + { + CCTK_REAL *const output_array_ptr_real + = (CCTK_REAL *) output_arrays[out]; + output_array_ptr_real[pt] = (CCTK_REAL) result; + break; + } + +#ifdef HAVE_CCTK_REAL4 +case CCTK_VARIABLE_REAL4: + { + CCTK_REAL4 *const output_array_ptr_real4 + = (CCTK_REAL4 *) output_arrays[out]; + output_array_ptr_real4[pt] = (CCTK_REAL4) result; + break; + } +#endif + +#ifdef HAVE_CCTK_REAL8 +case CCTK_VARIABLE_REAL8: + { + CCTK_REAL8 *const output_array_ptr_real8 + = (CCTK_REAL8 *) output_arrays[out]; + output_array_ptr_real8[pt] = (CCTK_REAL8) result; + break; + } +#endif + +#ifdef HAVE_CCTK_REAL16 +case CCTK_VARIABLE_REAL16: + { + CCTK_REAL16 *const output_array_ptr_real16 + = (CCTK_REAL16 *) output_arrays[out]; + output_array_ptr_real16[pt] = (CCTK_REAL16) result; + break; + } +#endif + +case CCTK_VARIABLE_COMPLEX: + { + CCTK_REAL (*const output_array_ptr_complex)[COMPLEX_N_PARTS] + = (CCTK_REAL (*)[COMPLEX_N_PARTS]) output_arrays[out]; + output_array_ptr_complex[pt][part] = (CCTK_REAL) result; + break; + } + +#ifdef HAVE_CCTK_COMPLEX8 +case CCTK_VARIABLE_COMPLEX8: + { + CCTK_REAL4 (*const output_array_ptr_complex8)[COMPLEX_N_PARTS] + = (CCTK_REAL4 (*)[COMPLEX_N_PARTS]) output_arrays[out]; + output_array_ptr_complex8[pt][part] = (CCTK_REAL4) result; + break; + } +#endif /* HAVE_CCTK_COMPLEX8 */ + +#ifdef HAVE_CCTK_COMPLEX16 +case CCTK_VARIABLE_COMPLEX16: + { + CCTK_REAL8 (*const output_array_ptr_complex16)[COMPLEX_N_PARTS] + = (CCTK_REAL8 (*)[COMPLEX_N_PARTS]) output_arrays[out]; + output_array_ptr_complex16[pt][part] = (CCTK_REAL8) result; + break; + } +#endif /* HAVE_CCTK_COMPLEX16 */ + +#ifdef HAVE_CCTK_COMPLEX32 +case CCTK_VARIABLE_COMPLEX32: + { + CCTK_REAL16 (*const output_array_ptr_complex32)[COMPLEX_N_PARTS] + = (CCTK_REAL16 (*)[COMPLEX_N_PARTS]) output_arrays[out]; + output_array_ptr_complex32[pt][part] = (CCTK_REAL16) result; + break; + } +#endif /* HAVE_CCTK_COMPLEX32 */ + +default: + CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, +"\n" +" CCTK_InterpLocalUniform(): output datatype %d not supported!\n" +" (0-origin) output #out=%d" + , + (int) output_array_type_codes[out], + out); + return UTIL_ERROR_BAD_INPUT; /*** ERROR RETURN ***/ + /* end of switch (output type code) */ + } + + /* end of if (output_arrays[out] != NULL) */ + } + + + /* + * handle querying the Jacobian + */ + if ( (Jacobian_info != NULL) + && (Jacobian_info->Jacobian_pointer[out] != NULL)) + then { + CCTK_REAL *const Jacobian_ptr + = Jacobian_info->Jacobian_pointer[out] + + Jacobian_info->Jacobian_offset[out] + + Jacobian_info->Jacobian_interp_point_stride*pt + + Jacobian_info->Jacobian_part_stride*part; + + switch (operation_codes[out]) + { + #ifdef HAVE_OP_I + case 0: + STORE_COEFFS(1.0, &coeffs_I, + Jacobian_ptr, + JACOBIAN_MIJK_STRIDE); + break; + #endif + #ifdef HAVE_OP_DX + case 1: + STORE_COEFFS(inverse_dx, &coeffs_dx, + Jacobian_ptr, + JACOBIAN_MIJK_STRIDE); + break; + #endif + #ifdef HAVE_OP_DY + case 2: + STORE_COEFFS(inverse_dy, &coeffs_dy, + Jacobian_ptr, + JACOBIAN_MIJK_STRIDE); + break; + #endif + #ifdef HAVE_OP_DZ + case 3: + STORE_COEFFS(inverse_dz, &coeffs_dz, + Jacobian_ptr, + JACOBIAN_MIJK_STRIDE); + break; + #endif + #ifdef HAVE_OP_DXX + case 11: + STORE_COEFFS(inverse_dx*inverse_dx, &coeffs_dxx, + Jacobian_ptr, + JACOBIAN_MIJK_STRIDE); + break; + #endif + #ifdef HAVE_OP_DXY + case 12: + case 21: + STORE_COEFFS(inverse_dx*inverse_dy, &coeffs_dxy, + Jacobian_ptr, + JACOBIAN_MIJK_STRIDE); + break; + #endif + #ifdef HAVE_OP_DXZ + case 13: + case 31: + STORE_COEFFS(inverse_dx*inverse_dz, &coeffs_dxz, + Jacobian_ptr, + JACOBIAN_MIJK_STRIDE); + break; + #endif + #ifdef HAVE_OP_DYY + case 22: + STORE_COEFFS(inverse_dy*inverse_dy, &coeffs_dyy, + Jacobian_ptr, + JACOBIAN_MIJK_STRIDE); + break; + #endif + #ifdef HAVE_OP_DYZ + case 23: + case 32: + STORE_COEFFS(inverse_dy*inverse_dz, &coeffs_dyz, + Jacobian_ptr, + JACOBIAN_MIJK_STRIDE); + break; + #endif + #ifdef HAVE_OP_DZZ + case 33: + STORE_COEFFS(inverse_dz*inverse_dz, &coeffs_dzz, + Jacobian_ptr, + JACOBIAN_MIJK_STRIDE); + break; + #endif + default: +CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, +"\n" +" CCTK_InterpLocalUniform(): generalized interpolation operation_code %d\n" +" is not supported!\n" +" (0-origin) output #out=%d" + , + (int) operation_codes[out], + out); + return UTIL_ERROR_BAD_INPUT; + /*** ERROR RETURN ***/ + /* end of switch(operation_codes[out])*/ + } + /* end of Jacobian-query code */ + } + + /* end of for (part = ...) loop */ + } + } + /* end of for (out = ...) loop */ + } + } + } + } + } + + /* end of for (pt = ...) loop */ + } + +return 0; /*** NORMAL RETURN ***/ + } + } +} + +/******************************************************************************/ + +/*@@ + @routine compute_coeffs_I + @routine compute_coeffs_dx + @routine compute_coeffs_dy + @routine compute_coeffs_dz + @routine compute_coeffs_dxx + @routine compute_coeffs_dxy + @routine compute_coeffs_dxz + @routine compute_coeffs_dyy + @routine compute_coeffs_dyz + @routine compute_coeffs_dzz + @date 29.Aug.2002 + @author Jonathan Thornburg <jthorn@aei.mpg.de> + @desc + Each of these functions computes the corresponding set of + interpolation coefficients at the point given by the XYZ + variables, using the machine-generated experessions in the + files + COEFFS_{I,DX,DY,DXX,DXY,DYY,...}_COMPUTE_FILE_NAME + + These functions *must* be static, because this whole file + is #included (and hence compiled) multiple times. + @enddesc + @@*/ + +#ifdef HAVE_OP_I + static + void compute_coeffs_I(FP_XYZ, struct COEFFS_STRUCT *coeffs_I) + { + #include COEFFS_I_COMPUTE_FILE_NAME + } +#endif + +#ifdef HAVE_OP_DX + static + void compute_coeffs_dx(FP_XYZ, struct COEFFS_STRUCT *coeffs_dx) + { + #include COEFFS_DX_COMPUTE_FILE_NAME + } +#endif + +#ifdef HAVE_OP_DY + static + void compute_coeffs_dy(FP_XYZ, struct COEFFS_STRUCT *coeffs_dy) + { + #include COEFFS_DY_COMPUTE_FILE_NAME + } +#endif + +#ifdef HAVE_OP_DZ + static + void compute_coeffs_dz(FP_XYZ, struct COEFFS_STRUCT *coeffs_dz) + { + #include COEFFS_DZ_COMPUTE_FILE_NAME + } +#endif + +#ifdef HAVE_OP_DXX + static + void compute_coeffs_dxx(FP_XYZ, struct COEFFS_STRUCT *coeffs_dxx) + { + #include COEFFS_DXX_COMPUTE_FILE_NAME + } +#endif + +#ifdef HAVE_OP_DXY + static + void compute_coeffs_dxy(FP_XYZ, struct COEFFS_STRUCT *coeffs_dxy) + { + #include COEFFS_DXY_COMPUTE_FILE_NAME + } +#endif + +#ifdef HAVE_OP_DXZ + static + void compute_coeffs_dxz(FP_XYZ, struct COEFFS_STRUCT *coeffs_dxz) + { + #include COEFFS_DXZ_COMPUTE_FILE_NAME + } +#endif + +#ifdef HAVE_OP_DYY + static + void compute_coeffs_dyy(FP_XYZ, struct COEFFS_STRUCT *coeffs_dyy) + { + #include COEFFS_DYY_COMPUTE_FILE_NAME + } +#endif + +#ifdef HAVE_OP_DYZ + static + void compute_coeffs_dyz(FP_XYZ, struct COEFFS_STRUCT *coeffs_dyz) + { + #include COEFFS_DYZ_COMPUTE_FILE_NAME + } +#endif + +#ifdef HAVE_OP_DZZ + static + void compute_coeffs_dzz(FP_XYZ, struct COEFFS_STRUCT *coeffs_dzz) + { + #include COEFFS_DZZ_COMPUTE_FILE_NAME + } +#endif |