/*@@ @file template.c @date 23 Oct 2001 @author Jonathan Thornburg @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, which does generalized interpolation for a single combination of (N_dims, molecule_family, order, smoothing). The following header files must be #included before #including this file: "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 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,...} 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. @endvar @var DATA_VAR_DCL_FILE_NAME @vdesc The name of a file (presumably machine-generated) containing a sequence of one or more C declarations for a molecule-sized set of "data variables", eg (for 2D size-3 molecules) fp data_m1_m1; fp data_0_m1; fp data_p1_m1; fp data_m1_0; fp data_0_0; fp data_p1_0; fp data_m1_p1; fp data_0_p1; fp data_p1_p1; @endvar @var COEFF_{I,DX,DY,DXX,DXY,DYY,...}_DCL_FILE_NAME @vdesc Each of these macros should be the name of a file (presumably machine-generated) containing a sequence of one or more C declarations for a molecule-sized set of coefficients for the corresponding derivative operator, eg (for dx with 2D size-3 molecules) fp coeff_dx_m1_m1; fp coeff_dx_0_m1; fp coeff_dx_p1_m1; fp coeff_dx_m1_0; fp coeff_dx_0_0; fp coeff_dx_p1_0; fp coeff_dx_m1_p1; fp coeff_dx_0_p1; fp coeff_dx_p1_p1; @endvar @var COEFF_{I,DX,DY,DXX,DXY,DYY,...}_VAR_STORE_FILE_NAME @vdesc The name of a file (presumably machine-generated) containing a sequence of C assignment statements to assign the interpolation coefficients to the corresponding COEFF(...) expressions, eg (for 2D size-3 molecules) COEFF(-1,-1) = coeff_dx_m1_m1; COEFF(0,-1) = coeff_dx_0_m1; COEFF(1,-1) = coeff_dx_p1_m1; COEFF(-1,0) = coeff_dx_m1_0; COEFF(0,0) = coeff_dx_0_0; COEFF(1,0) = coeff_dx_p1_0; COEFF(-1,1) = coeff_dx_m1_p1; COEFF(0,1) = coeff_dx_0_p1; COEFF(1,1) = coeff_dx_p1_p1; @endvar @var DATA_VAR_ASSIGN_FILE_NAME @vdesc The name of a file (presumably machine-generated) containing a sequence of C assignment statements to assign DATA(...) to the corresponding data variables for each molecule component, eg (for 2D size-3 molecules) data_m1_m1 = DATA(-1,-1); data_0_m1 = DATA(0,-1); data_p1_m1 = DATA(1,-1); data_m1_0 = DATA(-1,0); data_0_0 = DATA(0,0); data_p1_0 = DATA(1,0); data_m1_p1 = DATA(-1,+1); data_0_p1 = DATA(0,+1); data_p1_p1 = DATA(1,+1); @endvar @var INTERP_{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 C assignment statement (or a sequence of such statements) to compute the variable result as the molecule-sized linear combination of the data variables (cf DATA_VAR_{DCL,ASSIGN}_FILE_NAME) and the interpolation coefficients (cf COEFF_{I,DX,DY,DXX,DXY,DYY,...}_COMPUTE_FILE_NAME), eg (for 2D size-3 molecules, dx operator) result = coeff_dx_m1_m1*data_m1_m1 + coeff_dx_0_m1*data_0_m1 + coeff_dx_p1_m1*data_p1_m1 + coeff_dx_p2_m1*data_p2_m1 + coeff_dx_m1_0*data_m1_0 + coeff_dx_0_0*data_0_0 + coeff_dx_p1_0*data_p1_0 + coeff_dx_p2_0*data_p2_0 + coeff_dx_m1_p1*data_m1_p1 + coeff_dx_0_p1*data_0_p1 + coeff_dx_p1_p1*data_p1_p1 + coeff_dx_p2_p1*data_p2_p1 + coeff_dx_m1_p2*data_m1_p2 + coeff_dx_0_p2*data_0_p2 + coeff_dx_p1_p2*data_p1_p2 + coeff_dx_p2_p2*data_p2_p2; Note that this may *NOT* include any variable declarations; it should be valid to appear in the middle of a sequence of C statements. @endvar @var COEFF_{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 (molecule-sized set of) interpolation coefficients for the specified derivative operator, as polynomials in the variables (x,y,z), eg (for 2D size-3 molecules, dx operator) fp t35, t42, t41, t40, t39, t36; t35 = RATIONAL(1.0,3.0)*x; t42 = t35+RATIONAL(-1.0,4.0)*y; t41 = t35+RATIONAL(1.0,4.0)*y; t40 = RATIONAL(-1.0,6.0); t39 = RATIONAL(1.0,6.0); t36 = RATIONAL(-2.0,3.0)*x; coeff_dx_m1_m1 = t40+t41; coeff_dx_0_m1 = t36; coeff_dx_p1_m1 = t39+t42; coeff_dx_m1_0 = t40+t35; coeff_dx_0_0 = t36; coeff_dx_p1_0 = t35+t39; coeff_dx_m1_p1 = t40+t42; coeff_dx_0_p1 = t36; coeff_dx_p1_p1 = t39+t41; 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 these are the only coefficients which depend on the actual interpolation scheme used; all the others just depend on the interpolation dimension and molecule family and size. @endvar @version $Header$ @@*/ /******************************************************************************/ /*@@ @routine FUNCTION_NAME @date 23 Oct 2001 @author Jonathan Thornburg @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. @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_X_RANGE interpolation point is out of range @endreturndesc @@*/ int FUNCTION_NAME(/***** coordinate system *****/ const CCTK_REAL coord_system_origin[], const CCTK_REAL grid_spacing[], /***** interpolation points *****/ int N_interp_points, int interp_coords_type_code, const void* const interp_coords[], const CCTK_REAL out_of_range_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) { /* * 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 * for (int pt = 0 ; pt < N_interp_point ; ++pt) * { * declare all the coefficients * declare all the data-values variables * ***fetch*** interpolation point coordinates * compute position of interpolation molecules * with respect to the grid * if (querying molecule positions) * then store this molecule position * compute coefficients for all derivatives which are wanted * for (int out = 0 ; out < N_output_arrays ; ++out) * { * const int in = operand_indices[out]; * ***decode*** the input/output array datatypes * to determine whether they're real or complex * (they must both be the same in this regard), then * const int N_parts = data is complex ? 2 : 1; * for (int part = 0 ; part < N_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 * ***fetch*** the molecule-sized piece * of input_arrays[in][part] * at this molecule position * into local "data" variables * } * if (output_arrays[out] != NULL) * then { * fp result; * switch (operation_code) * { * case 0: * result = compute interpolant * as a linear combination * of data variables * and op=0 coeffs * break; * case 1: * result = compute interpolant * as a linear combination * of data variables * and op=1 coeffs * result *= inverse_dx; * break; * case ... * } * ***store*** result in output array * } * if (querying Jacobian * && (Jacobian_pointer[out] != NULL)) * then { * fp factor; * switch (operation_code) * { * case 0: * store op=0 Jacobian values * break; * case 1: * factor = inverse_dx; * store factor * * op=1 Jacobian values * break; * case ... * } * } * } * } * } * * Here "***decode***", "***fetch***", and "***store***" are all * actually switches on the various array datatypes. For complex * datatypes "***fetch***" and "***store***" 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. * * At present 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 "truncation" interpolation errors. */ /* * 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 */ /* number of real "parts" in a complex number */ #define COMPLEX_N_PARTS 2 /* 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 /* input array size, strides, and subscripting computation */ #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 == 1) #define SUB1(i) (i*stride_i) #elif (N_DIMS == 2) #define SUB2(i,j) (i*stride_i + j*stride_j) #elif (N_DIMS == 3) #define SUB3(i,j,k) (i*stride_i + j*stride_j + k*stride_k) #else #error "N_DIMS must be 1, 2, or 3!" #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) /* * Jacobian structure info */ 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 flags specifying which derivatives are wanted */ { #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 const fp origin_x = coord_system_origin[X_AXIS]; const fp delta_x = grid_spacing[X_AXIS]; #if defined(HAVE_OP_DX) || defined(HAVE_OP_DXY) || defined(HAVE_OP_DXZ) const fp inverse_delta_x = 1.0 / delta_x; #endif #ifdef HAVE_OP_DXX const fp inverse_delta_x2 = 1.0 / (delta_x*delta_x); #endif #endif #if N_DIMS >= 2 const fp origin_y = coord_system_origin[Y_AXIS]; const fp delta_y = grid_spacing[Y_AXIS]; #if defined(HAVE_OP_DY) || defined(HAVE_OP_DXY) || defined(HAVE_OP_DYZ) const fp inverse_delta_y = 1.0 / delta_y; #endif #ifdef HAVE_OP_DYY const fp inverse_delta_y2 = 1.0 / (delta_y*delta_y); #endif #endif #if N_DIMS >= 3 const fp origin_z = coord_system_origin[Z_AXIS]; const fp delta_z = grid_spacing[Z_AXIS]; #if defined(HAVE_OP_DZ) || defined(HAVE_OP_DXZ) || defined(HAVE_OP_DYZ) const fp inverse_delta_z = 1.0 / delta_z; #endif #ifdef HAVE_OP_DZZ const fp inverse_delta_z2 = 1.0 / (delta_z*delta_z); #endif #endif /* * interpolate at each point */ int pt; for (pt = 0 ; pt < N_interp_points ; ++pt) { /* declare all the data-values variables */ #include DATA_VAR_DCL_FILE_NAME /* declare all the interpolation coefficients */ #ifdef HAVE_OP_I #include COEFF_I_DCL_FILE_NAME #endif #ifdef HAVE_OP_DX #include COEFF_DX_DCL_FILE_NAME #endif #ifdef HAVE_OP_DY #include COEFF_DY_DCL_FILE_NAME #endif #ifdef HAVE_OP_DZ #include COEFF_DZ_DCL_FILE_NAME #endif #ifdef HAVE_OP_DXX #include COEFF_DXX_DCL_FILE_NAME #endif #ifdef HAVE_OP_DXY #include COEFF_DXY_DCL_FILE_NAME #endif #ifdef HAVE_OP_DXZ #include COEFF_DXZ_DCL_FILE_NAME #endif #ifdef HAVE_OP_DYY #include COEFF_DYY_DCL_FILE_NAME #endif #ifdef HAVE_OP_DYZ #include COEFF_DYZ_DCL_FILE_NAME #endif #ifdef HAVE_OP_DZZ #include COEFF_DZZ_DCL_FILE_NAME #endif /* * ***fetch*** interpolation point coordinates */ 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 */ } /* * 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 (likely) * better performance: the temp variables have their addresses * taken and so may not be register-allocated, whereas the * final (x,y,z) are "clean" and thus more likely to be * register-allocated */ { #if (N_DIMS >= 1) fp x_temp; const int center_i = LocalInterp_molecule_posn(origin_x, delta_x, input_array_min_subscripts[X_AXIS], input_array_max_subscripts[X_AXIS], MOLECULE_SIZE, out_of_range_tolerance[X_AXIS_MIN], out_of_range_tolerance[X_AXIS_MAX], interp_coords_fp[X_AXIS], &x_temp, (int *) NULL, (int *) NULL); const fp x = x_temp; #endif #if (N_DIMS >= 2) fp y_temp; const int center_j = LocalInterp_molecule_posn(origin_y, delta_y, input_array_min_subscripts[Y_AXIS], input_array_max_subscripts[Y_AXIS], MOLECULE_SIZE, out_of_range_tolerance[Y_AXIS_MIN], out_of_range_tolerance[Y_AXIS_MAX], interp_coords_fp[Y_AXIS], &y_temp, (int *) NULL, (int *) NULL); const fp y = y_temp; #endif #if (N_DIMS >= 3) fp z_temp; const int center_k = LocalInterp_molecule_posn(origin_z, delta_z, input_array_min_subscripts[Z_AXIS], input_array_max_subscripts[Z_AXIS], MOLECULE_SIZE, out_of_range_tolerance[Z_AXIS_MIN], out_of_range_tolerance[Z_AXIS_MAX], interp_coords_fp[Z_AXIS], &z_temp, (int *) NULL, (int *) NULL); const fp z = z_temp; #endif /* is the interpolation point out-of-range? */ #if (N_DIMS >= 1) if ((center_i == INT_MIN) || (center_i == INT_MAX)) then { if (error_flags != NULL) then { error_flags->error_pt = pt; error_flags->error_axis = X_AXIS; error_flags->error_end = (center_i > 0) ? +1 : -1; } return CCTK_ERROR_INTERP_POINT_X_RANGE; /*** ERROR RETURN ***/ } #endif #if (N_DIMS >= 2) if ((center_j == INT_MIN) || (center_j == INT_MAX)) then { if (error_flags != NULL) then { error_flags->error_pt = pt; error_flags->error_axis = Y_AXIS; error_flags->error_end = (center_j > 0) ? +1 : -1; } return CCTK_ERROR_INTERP_POINT_X_RANGE; /*** ERROR RETURN ***/ } #endif #if (N_DIMS >= 3) if ((center_k == INT_MIN) || (center_k == INT_MAX)) then { if (error_flags != NULL) then { error_flags->error_pt = pt; error_flags->error_axis = Z_AXIS; error_flags->error_end = (center_k > 0) ? +1 : -1; } return CCTK_ERROR_INTERP_POINT_X_RANGE; /*** ERROR RETURN ***/ } #endif /* 1D subscript of molecule center in input arrays */ /* (this is used in DATA data-fetching macros below) */ { #if (N_DIMS == 1) const int center_sub = SUB1(center_i); #elif (N_DIMS == 2) const int center_sub = SUB2(center_i, center_j); #elif (N_DIMS == 3) const int center_sub = SUB3(center_i, center_j, center_k); #else #error "N_DIMS must be 1, 2, or 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 for whichever derivatives are wanted * using machine-generated coefficient expressions * (these are polynomials in the variables (x,y,z)) */ #ifdef HAVE_OP_I if (want_I) then { #include COEFF_I_COMPUTE_FILE_NAME } #endif #ifdef HAVE_OP_DX if (want_dx) then { #include COEFF_DX_COMPUTE_FILE_NAME } #endif #ifdef HAVE_OP_DY if (want_dy) then { #include COEFF_DY_COMPUTE_FILE_NAME } #endif #ifdef HAVE_OP_DZ if (want_dz) then { #include COEFF_DZ_COMPUTE_FILE_NAME } #endif #ifdef HAVE_OP_DXX if (want_dxx) then { #include COEFF_DXX_COMPUTE_FILE_NAME } #endif #ifdef HAVE_OP_DXY if (want_dxy) then { #include COEFF_DXY_COMPUTE_FILE_NAME } #endif #ifdef HAVE_OP_DXZ if (want_dxz) then { #include COEFF_DXZ_COMPUTE_FILE_NAME } #endif #ifdef HAVE_OP_DYY if (want_dyy) then { #include COEFF_DYY_COMPUTE_FILE_NAME } #endif #ifdef HAVE_OP_DYZ if (want_dyz) then { #include COEFF_DYZ_COMPUTE_FILE_NAME } #endif #ifdef HAVE_OP_DZZ if (want_dzz) then { #include COEFF_DZZ_COMPUTE_FILE_NAME } #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 int input_center_sub = input_array_offsets[in] + center_sub; /* * ***decode*** the input/output array datatypes * to determine whether they're real or complex, * and verify that they're both the same in this regard * ==> define * const int N_parts = data is complex ? 2 : 1; */ const int N_input_parts = LocalInterp_decode_N_parts(input_array_type_codes[in]); const int N_output_parts = LocalInterp_decode_N_parts(output_array_type_codes[out]); if (N_input_parts != N_output_parts) then { CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, "\n" " CCTK_InterpLocalUniform(): can't do real input --> complex output or\n" " complex input --> real output interpolation!\n" " (0-origin) input #in=%d datatype = %d\n" " (0-origin) output #out=%d datatype = %d\n" , in, (int) input_array_type_codes[in], out, (int) output_array_type_codes[out]); /*NOTREACHED*/ return UTIL_ERROR_BAD_INPUT; /*** ERROR RETURN ***/ } { const int N_parts = N_input_parts; const void* const input_array_ptr = input_arrays[in]; int part; for (part = 0 ; part < N_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; /* * ***fetch*** the molecule-sized piece * of input_arrays[in][part] * at this molecule position * into local "data" variables */ switch (input_array_type_codes[in]) { case CCTK_VARIABLE_REAL: { const CCTK_REAL *const input_array_ptr_real = (const CCTK_REAL *) input_array_ptr; #undef DATA #if (N_DIMS == 1) #define DATA(mi) \ input_array_ptr_real[input_center_sub + SUB1(mi)] #elif (N_DIMS == 2) #define DATA(mi,mj) \ input_array_ptr_real[input_center_sub + SUB2(mi,mj)] #elif (N_DIMS == 3) #define DATA(mi,mj,mk) \ input_array_ptr_real[input_center_sub + SUB3(mi,mj,mk)] #else #error "N_DIMS must be 1, 2, or 3!" #endif #include DATA_VAR_ASSIGN_FILE_NAME break; } #ifdef HAVE_CCTK_REAL4 case CCTK_VARIABLE_REAL4: { const CCTK_REAL4 *const input_array_ptr_real4 = (const CCTK_REAL4 *) input_array_ptr; #undef DATA #if (N_DIMS == 1) #define DATA(mi) \ input_array_ptr_real4[input_center_sub + SUB1(mi)] #elif (N_DIMS == 2) #define DATA(mi,mj) \ input_array_ptr_real4[input_center_sub + SUB2(mi,mj)] #elif (N_DIMS == 3) #define DATA(mi,mj,mk) \ input_array_ptr_real4[input_center_sub + SUB3(mi,mj,mk)] #else #error "N_DIMS must be 1, 2, or 3!" #endif #include DATA_VAR_ASSIGN_FILE_NAME break; } #endif /* HAVE_CCTK_REAL4 */ #ifdef HAVE_CCTK_REAL8 case CCTK_VARIABLE_REAL8: { const CCTK_REAL8 *const input_array_ptr_real8 = (const CCTK_REAL8 *) input_array_ptr; #undef DATA #if (N_DIMS == 1) #define DATA(mi) \ input_array_ptr_real8[input_center_sub + SUB1(mi)] #elif (N_DIMS == 2) #define DATA(mi,mj) \ input_array_ptr_real8[input_center_sub + SUB2(mi,mj)] #elif (N_DIMS == 3) #define DATA(mi,mj,mk) \ input_array_ptr_real8[input_center_sub + SUB3(mi,mj,mk)] #else #error "N_DIMS must be 1, 2, or 3!" #endif #include DATA_VAR_ASSIGN_FILE_NAME break; } #endif /* HAVE_CCTK_REAL8 */ #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; #undef DATA #if (N_DIMS == 1) #define DATA(mi) \ input_array_ptr_real16[input_center_sub + SUB1(mi)] #elif (N_DIMS == 2) #define DATA(mi,mj) \ input_array_ptr_real16[input_center_sub + SUB2(mi,mj)] #elif (N_DIMS == 3) #define DATA(mi,mj,mk) \ input_array_ptr_real16[input_center_sub + SUB3(mi,mj,mk)] #else #error "N_DIMS must be 1, 2, or 3!" #endif #include DATA_VAR_ASSIGN_FILE_NAME break; } #endif /* HAVE_CCTK_REAL16 */ case CCTK_VARIABLE_COMPLEX: { const CCTK_REAL (*const input_array_ptr_complex)[COMPLEX_N_PARTS] = (const CCTK_REAL (*)[COMPLEX_N_PARTS]) input_array_ptr; #undef DATA #if (N_DIMS == 1) #define DATA(mi) \ input_array_ptr_complex[input_center_sub + SUB1(mi)][part] #elif (N_DIMS == 2) #define DATA(mi,mj) \ input_array_ptr_complex[input_center_sub + SUB2(mi,mj)][part] #elif (N_DIMS == 3) #define DATA(mi,mj,mk) \ input_array_ptr_complex[input_center_sub + SUB3(mi,mj,mk)][part] #else #error "N_DIMS must be 1, 2, or 3!" #endif #include DATA_VAR_ASSIGN_FILE_NAME 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; #undef DATA #if (N_DIMS == 1) #define DATA(mi) \ input_array_ptr_complex8[input_center_sub + SUB1(mi)][part] #elif (N_DIMS == 2) #define DATA(mi,mj) \ input_array_ptr_complex8[input_center_sub + SUB2(mi,mj)][part] #elif (N_DIMS == 3) #define DATA(mi,mj,mk) \ input_array_ptr_complex8[input_center_sub + SUB3(mi,mj,mk)][part] #else #error "N_DIMS must be 1, 2, or 3!" #endif #include DATA_VAR_ASSIGN_FILE_NAME break; } #endif /* HAVE_CCTK_COMPLEX8 */ #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; #undef DATA #if (N_DIMS == 1) #define DATA(mi) \ input_array_ptr_complex16[input_center_sub + SUB1(mi)][part] #elif (N_DIMS == 2) #define DATA(mi,mj) \ input_array_ptr_complex16[input_center_sub + SUB2(mi,mj)][part] #elif (N_DIMS == 3) #define DATA(mi,mj,mk) \ input_array_ptr_complex16[input_center_sub + SUB3(mi,mj,mk)][part] #else #error "N_DIMS must be 1, 2, or 3!" #endif #include DATA_VAR_ASSIGN_FILE_NAME break; } #endif /* HAVE_CCTK_COMPLEX16 */ #ifdef HAVE_CCTK_COMPLEX32 case CCTK_VARIABLE_COMPLEX32: { const CCTK_REAL16 (*const input_array_ptr_complex32)[COMPLEX_N_PARTS] = (const CCTK_REAL16 (*)[COMPLEX_N_PARTS]) input_array_ptr; #undef DATA #if (N_DIMS == 1) #define DATA(mi) \ input_array_ptr_complex32[input_center_sub + SUB1(mi)][part] #elif (N_DIMS == 2) #define DATA(mi,mj) \ input_array_ptr_complex32[input_center_sub + SUB2(mi,mj)][part] #elif (N_DIMS == 3) #define DATA(mi,mj,mk) \ input_array_ptr_complex32[input_center_sub + SUB3(mi,mj,mk)][part] #else #error "N_DIMS must be 1, 2, or 3!" #endif #include DATA_VAR_ASSIGN_FILE_NAME break; } #endif /* HAVE_CCTK_COMPLEX32 */ 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 { /* * compute the interpolant itself * as a linear combination of the data variables */ fp result; switch (operation_codes[out]) { #ifdef HAVE_OP_I case 0: #include INTERP_I_COMPUTE_FILE_NAME break; #endif #ifdef HAVE_OP_DX case 1: #include INTERP_DX_COMPUTE_FILE_NAME result *= inverse_delta_x; break; #endif #ifdef HAVE_OP_DY case 2: #include INTERP_DY_COMPUTE_FILE_NAME result *= inverse_delta_y; break; #endif #ifdef HAVE_OP_DZ case 3: #include INTERP_DZ_COMPUTE_FILE_NAME result *= inverse_delta_z; break; #endif #ifdef HAVE_OP_DXX case 11: #include INTERP_DXX_COMPUTE_FILE_NAME result *= inverse_delta_x2; break; #endif #ifdef HAVE_OP_DXY case 12: case 21: #include INTERP_DXY_COMPUTE_FILE_NAME result *= inverse_delta_x * inverse_delta_y; break; #endif #ifdef HAVE_OP_DXZ case 13: case 31: #include INTERP_DXZ_COMPUTE_FILE_NAME result *= inverse_delta_x * inverse_delta_z; break; #endif #ifdef HAVE_OP_DYY case 22: #include INTERP_DYY_COMPUTE_FILE_NAME result *= inverse_delta_y2; break; #endif #ifdef HAVE_OP_DYZ case 23: case 32: #include INTERP_DYZ_COMPUTE_FILE_NAME result *= inverse_delta_y * inverse_delta_z; break; #endif #ifdef HAVE_OP_DZZ case 33: #include INTERP_DZZ_COMPUTE_FILE_NAME result *= inverse_delta_z2; 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 ***/ /* end of switch (operation_codes[out]) */ } /* * ***store*** the result in the 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] = 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] = 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] = 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] = 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] = 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] = 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] = 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] = 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); /*NOTREACHED*/ 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 { /* * set up the Jacobian storage addressing */ const int J_offset = Jacobian_info->Jacobian_offset[out]; const int J_pt_stride = Jacobian_info->Jacobian_interp_point_stride; const int J_part_stride = Jacobian_info->Jacobian_part_stride; double *const Jacobian_ptr = & Jacobian_info->Jacobian_pointer[out][J_offset + pt*J_pt_stride + part*J_part_stride]; #if (N_DIMS >= 1) const int J_mx_stride = Jacobian_info->Jacobian_m_strides[X_AXIS]; #endif #if (N_DIMS >= 2) const int J_my_stride = Jacobian_info->Jacobian_m_strides[Y_AXIS]; #endif #if (N_DIMS >= 3) const int J_mz_stride = Jacobian_info->Jacobian_m_strides[Z_AXIS]; #endif #if (N_DIMS == 1) #define COEFF(mi) Jacobian_ptr[mi*J_mx_stride] #elif (N_DIMS == 2) #define COEFF(mi,mj) Jacobian_ptr[mi*J_mx_stride + mj*J_my_stride] #elif (N_DIMS == 3) #define COEFF(mi,mj,mk) \ Jacobian_ptr[mi*J_mx_stride + mj*J_my_stride + mk*J_mz_stride] #else #error "N_DIMS must be 1, 2, or 3!" #endif { fp factor; switch (operation_codes[out]) { #ifdef HAVE_OP_I case 0: #include COEFF_I_VAR_STORE_FILE_NAME break; #endif #ifdef HAVE_OP_DX case 1: factor = inverse_delta_x; #include COEFF_DX_VAR_STORE_FILE_NAME break; #endif #ifdef HAVE_OP_DY case 2: factor = inverse_delta_y; #include COEFF_DY_VAR_STORE_FILE_NAME break; #endif #ifdef HAVE_OP_DZ case 3: factor = inverse_delta_z; #include COEFF_DZ_VAR_STORE_FILE_NAME break; #endif #ifdef HAVE_OP_DXX case 11: factor = inverse_delta_x2; #include COEFF_DXX_VAR_STORE_FILE_NAME break; #endif #ifdef HAVE_OP_DXY case 12: case 21: factor = inverse_delta_x * inverse_delta_y; #include COEFF_DXY_VAR_STORE_FILE_NAME break; #endif #ifdef HAVE_OP_DXZ case 13: case 31: factor = inverse_delta_x * inverse_delta_z; #include COEFF_DXZ_VAR_STORE_FILE_NAME break; #endif #ifdef HAVE_OP_DYY case 22: factor = inverse_delta_y2; #include COEFF_DYY_VAR_STORE_FILE_NAME break; #endif #ifdef HAVE_OP_DYZ case 23: case 32: factor = inverse_delta_y * inverse_delta_z; #include COEFF_DYZ_VAR_STORE_FILE_NAME break; #endif #ifdef HAVE_OP_DZZ case 33: factor = inverse_delta_z2; #include COEFF_DZZ_VAR_STORE_FILE_NAME 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 ***/ /* 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 ***/ } } }