diff options
Diffstat (limited to 'archive/2d.cube.order2.smooth0.c')
| -rw-r--r-- | archive/2d.cube.order2.smooth0.c | 534 |
1 files changed, 534 insertions, 0 deletions
diff --git a/archive/2d.cube.order2.smooth0.c b/archive/2d.cube.order2.smooth0.c new file mode 100644 index 0000000..34e17ca --- /dev/null +++ b/archive/2d.cube.order2.smooth0.c @@ -0,0 +1,534 @@ + /*@@ + @file 2d.cube.order2.smooth0.c + @date 23 Oct 2001 + @author Jonathan Thornburg <jthorn@aei.mpg.de> + @desc + Generalized interpolation for 2d, hypercube-shaped molecules, + order=2, smoothing=0. For details, see the header comments + for "InterpLocalArrays.c" in this directory. + @enddesc + + @version $Id$ + @@*/ + +#include <math.h> +#include <limits.h> +#include <stdlib.h> +#include <string.h> +#include <stdio.h> + +#include "util_ErrorCodes.h" +#include "cctk.h" +#include "InterpLocalArrays.h" + +/* the rcs ID and its dummy function to use it */ +static const char *rcsid = "$Header$"; +CCTK_FILEVERSION( + CactusPUGH_LocalInterp_GeneralizedPolynomial_2d_cube_order2_smooth0_c + ) + +/******************************************************************************/ + +/*@@ + @routine LocalInterp_ILA_2d_cube_ord2_sm0 + @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 InterpLocalArrays() (in "InterpLocalArrays.c" + in this same directory). + + This function's arguments are all a subset of those of + InterpLocalArrays() ; the only difference is that this function + takes all its arguments explicitly, whereas InputLocalArrays() + takes some of them indirectly via a key/value parameter table. + + @returntype int + @returndesc This function's return result is the same as that of + InterpLocalArrays(): + 0 ==> successful interpolation + -1 ==> in case of any errors + @endreturndesc + + @@*/ +int LocalInterp_ILA_2d_cube_ord2_sm0 + (int param_table_handle, + const CCTK_REAL coord_system_origin[], + const CCTK_REAL grid_spacing[], + int N_interp_points, + const CCTK_INT interp_coord_type_codes[], + const void *const interp_coords[], + 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[], + int N_output_arrays, + const CCTK_INT output_array_type_codes[], + void *const output_arrays[], + const CCTK_INT operand_indices[], const CCTK_INT opcodes[]) +{ +/* + * Implementation notes: + * + * The basic outline of this function is as follows: + * + * compute "which derivatives are wanted" flags + * for each interpolation point + * { + * declare all the coefficients + * declare all the data-values variables + * ***fetch*** interpolation point coordinates + * compute coefficients for all derivatives which are wanted + * for each output array + * { + * int part; + * for (part = 0 ; part <= 1 ; ++part) + * { + * if (this output array is computed + * using a different input array + * than the previous one || part != 0) + * then ***fetch*** the input array values + * into local variables + * { + * fp result; + * switch (opcode) + * { + * case 0: + * result = evaluate the interpolant + * break; + * case 1: + * result = evaluate the interpolant + * break; + * case ... + * } + * ***store*** result in output array + * bool complex_flag = is datatype complex? + * if (! complex_flag) + * then break; + * } + * } + * } + * } + * + * Here "***fetch***" and "***store***" are all actually switches on + * the various array datatypes. For complex datatypes they offset the + * 1D array position by part to handle 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: + * input, output = 0-origin indices each selecting an input/output array + * point = 0-origin index selecting an interpolation point + */ + +/* + * these are compile-time constants here; InterpLocalArrays() decoded + * them and called us (as opposed to another function) based in part + * on these values + */ +#define N_DIMS 2 +#define MOLECULE_SIZE 3 + +/* layout of axes in N_dims-element arrays */ +#define X_AXIS 0 +#define Y_AXIS 1 + +/* input array size, strides, and subscripting computation */ +const int stride_i = input_array_strides[X_AXIS]; +const int stride_j = input_array_strides[Y_AXIS]; +#define SUB2(i,j) (i*stride_i + j*stride_j) + +/* macros used by machine-generated interpolation coefficient expressions */ +#define RATIONAL(num,den) (num/den) + +/* + * compute flags specifying which derivatives are wanted + */ +bool want_I = false; +bool want_dx = false, want_dy = false; + { +int output; + for (output = 0 ; output < N_output_arrays ; ++output) + { + switch (opcodes[output]) + { + case 0: want_I = true; break; + case 1: want_dx = true; break; + case 2: want_dy = true; break; + default: + CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, + "Generalized interpolation opcode %d not supported", + opcodes[output]); /*NOTREACHED*/ + return UTIL_ERROR_BAD_INPUT; /*** ERROR RETURN ***/ + } + } + } + +/* + * interpolate at each point + */ + { +int point; + for (point = 0 ; point < N_interp_points ; ++point) + { + /* declare all the interpolation coefficients */ + #include "coeffs/2d.cube.order2.smooth0.I.dcl.c" + #include "coeffs/2d.cube.order2.smooth0.dx.dcl.c" + #include "coeffs/2d.cube.order2.smooth0.dy.dcl.c" + + /* declare all the data-values variables */ + #include "coeffs/2d.cube.size3.data-var.dcl.c" + + /* + * ***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_coord_type_codes[axis]) + { +#ifdef CCTK_REAL + 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[point]; + break; + } +#endif +#ifdef 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[point]; + break; + } +#endif +#ifdef 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[point]; + break; + } +#endif + default: + CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, + "interp-coords datatype %d not supported", + interp_coord_type_codes[axis]); + /*NOTREACHED*/ + return UTIL_ERROR_BAD_INPUT; /*** ERROR RETURN ***/ + } + } + + /* + * locate interpolation molecules with respect to the grid, + * i.e. compute interp_rel_(x,y) + */ + { + fp interp_rel_x, interp_rel_y; /* (x,y) coordinates of interpolation */ + /* point relative to molecule center */ + /* (in units of the grid spacing) */ + const int center_i + = LocalInterp_molecule_posn(coord_system_origin[X_AXIS], + grid_spacing[X_AXIS], + input_array_min_subscripts[X_AXIS], + input_array_max_subscripts[X_AXIS], + MOLECULE_SIZE, + interp_coords_fp[X_AXIS], + &interp_rel_x, + (int *) NULL, (int *) NULL); + const int center_j + = LocalInterp_molecule_posn(coord_system_origin[Y_AXIS], + grid_spacing[Y_AXIS], + input_array_min_subscripts[Y_AXIS], + input_array_max_subscripts[Y_AXIS], + MOLECULE_SIZE, + interp_coords_fp[Y_AXIS], + &interp_rel_y, + (int *) NULL, (int *) NULL); + const int center_sub = SUB2(center_i, center_j); + + /* + * compute the coefficients for whichever derivatives are wanted + * using machine-generated coefficient expressions + * ... these expressions are polynomials in (x,y) + * ==> we need these names for the relative coordinates + * (copying to fresh local variables will likely also + * give better register allocation, since [xy]_rel + * had their addresses taken and so probably won't be + * register-allocated) + */ + const fp x = interp_rel_x; + const fp y = interp_rel_y; + if (want_I) + then { + #include "coeffs/2d.cube.order2.smooth0.I.coeff.c" + } + if (want_dx) + then { + #include "coeffs/2d.cube.order2.smooth0.dx.coeff.c" + } + if (want_dy) + then { + #include "coeffs/2d.cube.order2.smooth0.dy.coeff.c" + } + + /* + * compute each output array at this point + */ + { + int output; + const void *input_array_ptr = NULL; + for (output = 0 ; output < N_output_arrays ; ++output) + { + const int input = operand_indices[output]; + const int input_offset = input_array_offsets[input]; + + /* + * for each real/imag part of complex data values + * ... for real we'll break out of this loop at the bottom + * after only a single iteration; + * for complex we'll do both iterations + */ + int part; + for (part = 0 ; part <= 1 ; ++part) + { + if ( (input_arrays[input] != input_array_ptr) + || (part != 0) ) + then { + /* + * ***fetch*** the input array values + * into local variables + */ + input_array_ptr = input_arrays[input]; + switch (input_array_type_codes[input]) + { +#ifdef CCTK_REAL +case CCTK_VARIABLE_REAL: + { + const CCTK_REAL *const input_array_ptr_real + = (const CCTK_REAL *) input_array_ptr; + #undef DATA + #define DATA(i,j) \ + input_array_ptr_real[input_offset + center_sub + SUB2(i,j)] + #include "coeffs/2d.cube.size3.data-var.assign.c" + break; + } +#endif +#ifdef CCTK_REAL4 +case CCTK_VARIABLE_REAL4: + { + const CCTK_REAL4 *const input_array_ptr_real4 + = (const CCTK_REAL4 *) input_array_ptr; + #undef DATA + #define DATA(i,j) \ + input_array_ptr_real4[input_offset + center_sub + SUB2(i,j)] + #include "coeffs/2d.cube.size3.data-var.assign.c" + break; + } +#endif +#ifdef CCTK_REAL8 +case CCTK_VARIABLE_REAL8: + { + const CCTK_REAL8 *const input_array_ptr_real8 + = (const CCTK_REAL8 *) input_array_ptr; + #undef DATA + #define DATA(i,j) \ + input_array_ptr_real8[input_offset + center_sub + SUB2(i,j)] + #include "coeffs/2d.cube.size3.data-var.assign.c" + break; + } +#endif +#ifdef CCTK_COMPLEX +case CCTK_VARIABLE_COMPLEX: + { + const CCTK_COMPLEX *const input_array_ptr_complex + = (const CCTK_COMPLEX *) input_array_ptr; + #undef DATA + #define DATA(i,j) \ + ( (const CCTK_REAL *) \ + & input_array_ptr_complex[input_offset + center_sub + SUB2(i,j)] ) [part] + #include "coeffs/2d.cube.size3.data-var.assign.c" + break; + } +#endif +#ifdef CCTK_COMPLEX8 +case CCTK_VARIABLE_COMPLEX8: + { + const CCTK_COMPLEX8 *const input_array_ptr_complex8 + = (const CCTK_COMPLEX8 *) input_array_ptr; + #undef DATA + #define DATA(i,j) \ + ( (const CCTK_REAL4 *) \ + & input_array_ptr_complex8[input_offset + center_sub + SUB2(i,j)] ) [part] + #include "coeffs/2d.cube.size3.data-var.assign.c" + break; + } +#endif +#ifdef CCTK_COMPLEX16 +case CCTK_VARIABLE_COMPLEX16: + { + const CCTK_COMPLEX16 *const input_array_ptr_complex16 + = (const CCTK_COMPLEX16 *) input_array_ptr; + #undef DATA + #define DATA(i,j) \ + ( (const CCTK_REAL8 *) \ + & input_array_ptr_complex16[input_offset + center_sub + SUB2(i,j)] ) [part] + #include "coeffs/2d.cube.size3.data-var.assign.c" + break; + } +#endif +default: + CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, + "input datatype %d not supported", + input_array_type_codes[input]); /*NOTREACHED*/ + return UTIL_ERROR_BAD_INPUT; /*** ERROR RETURN ***/ + } + } + + /* + * evaluate the interpolant + * as a linear combination of the data variables + */ + { + fp result; + switch (opcodes[output]) + { + case 0: + result = + #include "coeffs/2d.cube.order2.smooth0.I.eval.c" + break; + case 1: + result = + #include "coeffs/2d.cube.order2.smooth0.dx.eval.c" + break; + case 2: + result = + #include "coeffs/2d.cube.order2.smooth0.dy.eval.c" + break; + default: + CCTK_VWarn(1, __LINE__, __FILE__, + CCTK_THORNSTRING, + "opcode %d not supported", + opcodes[output]); /*NOTREACHED*/ + return UTIL_ERROR_BAD_INPUT; + /*** ERROR RETURN ***/ + } + + /* + * ***store*** the result in the output array + */ + switch (output_array_type_codes[output]) + { +case CCTK_VARIABLE_REAL: + { + CCTK_REAL *const output_array_ptr_real + = (CCTK_REAL *) output_arrays[output]; + output_array_ptr_real[point] = result; + break; + } +case CCTK_VARIABLE_REAL4: + { + CCTK_REAL4 *const output_array_ptr_real4 + = (CCTK_REAL4 *) output_arrays[output]; + output_array_ptr_real4[point] = result; + break; + } +case CCTK_VARIABLE_REAL8: + { + CCTK_REAL8 *const output_array_ptr_real8 + = (CCTK_REAL8 *) output_arrays[output]; + output_array_ptr_real8[point] = result; + break; + } +case CCTK_VARIABLE_COMPLEX: + { + CCTK_COMPLEX *const output_array_ptr_complex + = (CCTK_COMPLEX *) output_arrays[output]; + ((CCTK_REAL *) & output_array_ptr_complex[point]) [part] + = result; + break; + } +case CCTK_VARIABLE_COMPLEX8: + { + CCTK_COMPLEX8 *const output_array_ptr_complex8 + = (CCTK_COMPLEX8 *) output_arrays[output]; + ((CCTK_REAL4 *) & output_array_ptr_complex8[point]) [part] + = result; + break; + } +case CCTK_VARIABLE_COMPLEX16: + { + CCTK_COMPLEX16 *const output_array_ptr_complex16 + = (CCTK_COMPLEX16 *) output_arrays[output]; + ((CCTK_REAL8 *) & output_array_ptr_complex16[point]) [part] + = result; + break; + } +default: + CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, + "output datatype %d not supported", + output_array_type_codes[output]); /*NOTREACHED*/ + return UTIL_ERROR_BAD_INPUT; /*** ERROR RETURN ***/ + } + + /* decode datatype: is it real or complex? */ + { + bool complex_flag; + switch (output_array_type_codes[output]) + { + case CCTK_VARIABLE_REAL: + case CCTK_VARIABLE_REAL4: + case CCTK_VARIABLE_REAL8: + complex_flag = false; + break; + case CCTK_VARIABLE_COMPLEX: + case CCTK_VARIABLE_COMPLEX8: + case CCTK_VARIABLE_COMPLEX16: + complex_flag = true; + break; + default: + CCTK_VWarn(1, __LINE__, __FILE__, + CCTK_THORNSTRING, + "output datatype %d not supported", + output_array_type_codes[output]); + /*NOTREACHED*/ + return UTIL_ERROR_BAD_INPUT; + /*** ERROR RETURN ***/ + } + + /* skip part=1 (imaginary part) for real datatypes */ + if (! complex_flag) + then break; + } + } + /* end of for (part = ...) loop */ + } + /* end of for (output = ...) loop */ + } + } + } + + /* end of for (point = ...) loop */ + } + +return 0; /*** NORMAL RETURN ***/ + } +} |