aboutsummaryrefslogtreecommitdiff
path: root/src/GeneralizedPolynomial-Uniform/template.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/GeneralizedPolynomial-Uniform/template.c')
-rw-r--r--src/GeneralizedPolynomial-Uniform/template.c1675
1 files changed, 0 insertions, 1675 deletions
diff --git a/src/GeneralizedPolynomial-Uniform/template.c b/src/GeneralizedPolynomial-Uniform/template.c
deleted file mode 100644
index 454f733..0000000
--- a/src/GeneralizedPolynomial-Uniform/template.c
+++ /dev/null
@@ -1,1675 +0,0 @@
-/*@@
- @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 LocalInterp_ILA_2d_cube_ord2_sm0
- @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 FETCH_DATA_{REAL,REAL{4,8,16},COMPLEX,COMPLEX{8,16,32}}
- @vdesc The name of a C function to fetch data from the input arrays
- into the DATA_STRUCT structure. Typically these will be
- chosem from among the functions defined in
- "common/fetch-template.[ch]", e.g.
- #define FETCH_DATA_REAL LocalInterp_fetch_2d_cube3_r
- #define FETCH_DATA_REAL4 LocalInterp_fetch_2d_cube3_r4
- #define FETCH_DATA_REAL8 LocalInterp_fetch_2d_cube3_r8
- #define FETCH_DATA_REAL16 LocalInterp_fetch_2d_cube3_r16
- #define FETCH_DATA_COMPLEX LocalInterp_fetch_2d_cube3_c
- #define FETCH_DATA_COMPLEX8 LocalInterp_fetch_2d_cube3_c8
- #define FETCH_DATA_COMPLEX16 LocalInterp_fetch_2d_cube3_c16
- #define FETCH_DATA_COMPLEX32 LocalInterp_fetch_2d_cube3_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:
- * fetch 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 fetch)
- * || (part != value at last fetch) ) )
- * 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:
- * FETCH_DATA_REAL(input_array_ptr,
- * STRIDE_IJK,
- * &data);
- * break;
- * ...
- * case CCTK_VARIABLE_COMPLEX:
- * FETCH_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
-
-
- /*
- * fetch 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 = 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],
- &center_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 ***fetch*** the first time in the test at the
- * top of the part loop below
- */
- const void* input_array_ptr__last_fetch = NULL;
- int part__last_fetch = -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
- = LocalInterp_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_fetch)
- || (part != part__last_fetch)) )
- then {
- /* remember when we did the following fetch */
- input_array_ptr__last_fetch = input_array_ptr;
- part__last_fetch = part;
-
- /*
- * ***decode*** the input array datatype
- * to determine whether it's real or complex,
- */
- {
- const int N_input_parts
- = LocalInterp_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 ***/
- }
-
- /*
- * fetch 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;
- FETCH_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;
- FETCH_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;
- FETCH_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;
- FETCH_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;
- FETCH_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;
- FETCH_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;
- FETCH_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;
- FETCH_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 fetch 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
generated by cgit v1.2.3 (git 2.39.1) at 2024-05-21 06:50:13 +0000