diff options
Diffstat (limited to 'src/Interpolate.c')
| -rw-r--r-- | src/Interpolate.c | 615 |
1 files changed, 0 insertions, 615 deletions
diff --git a/src/Interpolate.c b/src/Interpolate.c deleted file mode 100644 index b46cb40..0000000 --- a/src/Interpolate.c +++ /dev/null @@ -1,615 +0,0 @@ - /*@@ - @file Interpolate.c - @date Wed 17 Jan 2001 - @author Thomas Radke - @desc - Interpolation of arrays to arbitrary points - - This interpolator is based on the Cactus 3.x Fortran version - written by Paul Walker. It also contains some nice optimization - features from Erik Schnetter. Jonathan Thornburg added some - additional comments in October 2001. - @enddesc - - @history - @date Wed 17 Jan 2001 - @author Thomas Radke - @hdesc Translation from Fortran to C - @date Thu 18 Oct 2001 - @author Jonathan Thornburg - @hdesc Add lots of comments, PUGHINTERP_VERBOSE_DEBUG debugging code - @endhistory - @version $Id$ - @@*/ - -#include <assert.h> -#include <math.h> -#include <stdlib.h> -#include <string.h> -#include <stdio.h> - -#include "cctk.h" -#include "cctk_Parameters.h" -#include "pughInterpGH.h" - -/* the rcs ID and its dummy function to use it */ -static const char *rcsid = "$Header$"; -CCTK_FILEVERSION(CactusPUGH_PUGHInterp_Interpolate_c) - -#ifdef PUGHINTERP_VERBOSE_DEBUG - /* if this is >= 0, we print verbose debugging information at this point */ - int PUGHInterp_verbose_debug_n = -1; -#endif - -/* the highest order of interpolation we support so far */ -#define MAXORDER 3 - -/* the highest dimension for variables we can deal with (so far) */ -#define MAXDIM 3 - -/******************************************************************************/ - -/*@@ - @routine INTERPOLATE (macro) - @date 18 Oct 2001 - @author code by ???, these comments by Jonathan Thornburg - @desc - This macro does the interpolation of in_array[] to compute - a single value out_array[n] (actually out_array[n]subpart ; - see the comments below for details). The data to be interpolated - must be real numbers of some type, i.e. if the arrays are - complex this macro must be called separately for the real - and imaginary parts. - @enddesc - - @var cctk_type - @vdesc C type of input and output array elements (might be complex) - @endvar - - @var cctk_subtype - @vdesc C type of actual numbers being interpolated (must be real) - @endvar - - @var subpart - @vdesc string to be suffixed to input/output array element to get - to get real number, i.e. empty string if cctk_type is real, - .Re or .Im as appropriate if cctk_type is complex - @endvar - - @var in_array - @vdesc A pointer to array to be interpolated (strictly speaking, to - the array's [0][0]...[0] element); this is typically passed - as a void * pointer so we typecast it as necessary. - @endvar - - @var out_array - @vdesc A 1-dimensional array where the interpolation result should be - stored; this is typically passed as a void * pointer so we - typecast it as necessary. - @endvar - - @var order - @vdesc The order of the interpolation (1=linear, 2=quadratic, 3=cubic, ...) - @endvar - - @var point - @vdesc [MAXDIM] array of integers giving the integer grid coordinates - of the closest grid point to the interpolation point; the - interpolation molecule/stencil is centered at this point. - @endvar - - @var dims - @vdesc [MAXDIM] array of integers giving the dimensions of in_array . - @endvar - - @var n - @vdesc Position in out_array where we should store the interpolation - result. - @endvar - - @var coeff - @vdesc [MAXDIM][MAX_ORDER+1] array of (floating-point) interpolation - coefficients; detailed semantics are that coeff[axis][m] is the - coefficient of y[m] when the 1-dimensional Lagrange interpolation - polynomial passing through the order+1 points - {(0,y[0]), (1,y[1]), ..., (order,y[order])} - is evaluated at the position x=offset[axis]. - @endvar - @@*/ -/* - * The basic idea here is that conceptually we first interpolate the - * (say) 3D gridfn in the x direction at each y and z grid point, - * then interpolate that 2D plane of values in the y direction at - * each z grid point, and finally interpolate that 1D line of values - * in the z direction. The implementation actually interleaves the - * different directions' interpolations so that only 3 scalar temporaries - * are needed. - */ -#define INTERPOLATE(cctk_type, cctk_subtype, subpart, in_array, out_array, \ - order, point, dims, n, coeff) \ - { \ - int ii, jj, kk; \ - const cctk_type *fi; \ - /* for some reason (probably loop unrolling & pipelining) the compiler \ - produces faster code if arrays are used instead of scalars */ \ - cctk_subtype interp_result, fj[1], fk[1]; \ - \ - \ - interp_result = 0; \ - \ - /* NOTE-MAXDIM: support >3D arrays by adding more loops */ \ - for (kk = 0; kk <= order; kk++) \ - /*if (fabs(coeff[2][kk]) > 1e-6)*/ \ - { \ - fk[0] = 0; \ - for (jj = 0; jj <= order; jj++) \ - /*if (fabs(coeff[1][jj]) > 1e-6)*/ \ - { \ - /* NOTE-MAXDIM: for >3D arrays adapt the index calculation here */ \ - fi = (const cctk_type *) in_array + \ - point[0] + dims[0]*(point[1]+jj + dims[1]*(point[2]+kk)); \ - \ - fj[0] = 0; \ - for (ii = 0; ii <= order; ii++) \ - /*if (fabs(coeff[0][ii]) > 1e-6)*/ \ - { \ - fj[0] += fi[ii]subpart * coeff[0][ii]; \ - } \ - /* at this point we have just computed */ \ - /* fj[0] = in_array[*][jj][kk] interpolated to x=offset[0] */ \ - \ - fk[0] += fj[0] * coeff[1][jj]; \ - } \ - /* at this point we have just computed */ \ - /* fk[0] = fj[0][*][kk] interpolated to y=offset[1] */ \ - /* = in_array[*][*][kk] interpolated to */ \ - /* x=offset[0], y=offset[1] */ \ - \ - interp_result += fk[0] * coeff[2][kk]; \ - } \ - /* at this point we have just computed */ \ - /* interp_result = fk[0][*] interpolated to z=offset[2] */ \ - /* = in_array[*][*][*] interpolated to */ \ - /* x=offset[0], y=offset[1], z=offset[2] */ \ - \ - /* assign the result */ \ - ((cctk_type *) out_array)[n]subpart = interp_result; \ - } /* end of macro */ - -/* this is needed by some preprocessors to pass into INTERPOLATE - as a dummy macro */ -#define NOTHING - -/******************************************************************************/ - -/*@@ - @routine PUGHInterp_Interpolate - @date Wed 17 Jan 2001 - @author Thomas Radke - @desc - This routine interpolates a set of input arrays - to a set of output arrays (one-to-one) at arbitrary points - which are given by their coordinates and the underlying - regular, uniform grid. - - Current limitations of this implementation are: - - arrays up to three (MAXDIM) dimensions only can be handled - - interpolation orders up to three (MAXORDER) only are supported - - coordinates must be given as CCTK_REAL types - - input and output array types must be the same - (no type casting of interpolation results supported) - - Despite of these limitations, the code was programmed almost - generically in that it can easily be extended to support - higher-dimensional arrays or more interpolation orders. - Places where the code would need to be changed to do this, - are marked with NOTE-MAXDIM and/or NOTE-MAXORDER comments - as appropriate. - @enddesc - - @var num_points - @vdesc number of points to interpolate at - @vtype int - @vio in - @endvar - @var num_dims - @vdesc dimensionality of the input arrays - @vtype int - @vio in - @endvar - @var num_arrays - @vdesc number of input/output arrays - @vtype int - @vio in - @endvar - @var dims - @vdesc dimensions of the input arrays - @vtype int[ num_dims ] - @vio in - @endvar - @var coord - @vdesc coordinates to interpolate at - @vtype CCTK_REAL coord[ num_dims * num_points ] - @vio in - @endvar - @var origin - @vdesc origin of the underlying grid - @vtype CCTK_REAL origin[ num_dims ] - @vio in - @endvar - @var delta - @vdesc deltas of the underlying grid - @vtype CCTK_REAL delta[ num_dims ] - @vio in - @endvar - @var in_types - @vdesc CCTK variable types of input arrays - @vtype int in_types[ num_arrays ] - @vio in - @endvar - @var in_arrays - @vdesc list of input arrays - @vtype void *in_arrays[ num_arrays ] - @vio in - @endvar - @var out_types - @vdesc CCTK variable types of output arrays - @vtype int out_types[ num_arrays ] - @vio in - @endvar - @var out_arrays - @vdesc list of output arrays - @vtype void *out_arrays[ num_arrays ] - @vio out - @endvar - - @returntype int - @returndesc - 0 - successful interpolation - -1 - in case of any errors - @endreturndesc -@@*/ -int PUGHInterp_Interpolate (int order, - int num_points, - int num_dims, - int num_arrays, - const int dims[], - const CCTK_REAL coord[], - const CCTK_REAL origin[], - const CCTK_REAL delta[], - const int in_types[], - const void *const in_arrays[], - const int out_types[], - void *const out_arrays[]) -{ - int retval; - int i, a, n, shift; -#if 0 - int out_of_bounds; -#endif - int max_dims[MAXDIM], point[MAXDIM]; - CCTK_REAL delta_inv[MAXDIM]; - CCTK_REAL below[MAXDIM]; - CCTK_REAL offset[MAXDIM]; - CCTK_REAL coeff[MAXDIM][MAXORDER + 1]; - - - /* verify parameters and check against our restrictions */ - retval = -1; - if (num_dims < 1) - { - CCTK_WARN (1, "Number of dimensions must be positive"); - } - else if (num_dims > MAXDIM) - { - CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, - "Interpolation of %d-dimensional arrays not implemented", - num_dims); - } - else if (order < 1) - { - CCTK_WARN (1, "Inperpolation order must be positive"); - } - else if (order > MAXORDER) - { - CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, - "Interpolation order %d not implemented", order); - } - else if (num_points < 0) - { - CCTK_WARN (1, "Negative number of points given"); - } - else - { - retval = 0; - } - - /* immediately return in case of errors */ - if (retval) - { - return (retval); - } - - /* also immediately return if there's nothing to do */ - if (num_points == 0) - { - return (retval); - } - - /* avoid divisions by delta later on */ - for (i = 0; i < num_dims; i++) - { - delta_inv[i] = 1.0 / delta[i]; - } - - /* duplicate the dims[] vector into one with MAXDIM-1 elements - (with the remaining elements zeroed out) - so that we can use nested loops over MAXDIM dimensions later on */ - memset (max_dims, 0, sizeof (max_dims)); - memcpy (max_dims, dims, (num_dims - 1) * sizeof (int)); - - /* zero out the coefficients and set the elements with index 'order' to one - so that we can use nested loops over MAXDIM dimensions later on */ - memset (coeff, 0, sizeof (coeff)); - for (i = num_dims; i < MAXDIM; i++) - { - coeff[i][0] = 1; - } - - /* zero out the iterator */ - memset (point, 0, sizeof (point)); - - /* loop over all points to interpolate at */ - for (n = 0; n < num_points; n++) - { -#if 0 - /* reset the out-of-bounds flag */ - out_of_bounds = 0; -#endif - - /* loop over all dimensions */ - for (i = 0; i < num_dims; i++) - { - /* closest grid point for stencil/molecule */ - point[i] = floor ((coord[num_dims*n + i] - origin[i]) * delta_inv[i] - - 0.5 * (order - 1)); - -#if 0 - /* test bounds */ - out_of_bounds |= point[i] < 0 || point[i]+order >= dims[i]; -#endif - - /* check that the stencil/molecule isn't bigger than the grid */ - if (order+1 > dims[i]) /* stencil/molecule size = order+1 */ - { - return -1; - } - - /* if beyond lower bound shift the grid point to the right */ - shift = point[i]; - if (shift < 0) - { - point[i] -= shift; - } - - /* if beyond upper bound shift the grid point to the left */ - shift = point[i] + order - (dims[i] - 1); - if (shift > 0) - { - point[i] -= shift; - } - - assert (point[i] >= 0 && point[i]+order < dims[i]); - - /* physical coordinate of that grid point */ - below[i] = origin[i] + point[i] * delta[i]; - - /* offset from that grid point, in fractions of grid points */ - offset[i] = (coord[num_dims*n + i] - below[i]) * delta_inv[i]; - -#if 0 - /* This test ensures that it is really an interpolation, and not - an extrapolation */ - assert (offset[i]>=0.5*(order-1) && offset[i]<=0.5*(order-1)+1.0); -#endif - } - -#ifdef PUGHINTERP_VERBOSE_DEBUG -if (n == PUGHInterp_verbose_debug_n) - { - int ii; - printf("out_of_bounds = %d\n", out_of_bounds); - for (ii = 0 ; ii < num_dims ; ++ii) - { - printf("offset[%d] = %g\n", ii, (double) offset[ii]); - } - } -#endif /* PUGHINTERP_VERBOSE_DEBUG */ - -#if 0 - /* check bounds */ - if (out_of_bounds) - { - char *msg; - - - /* put all information into a single message string for output */ - msg = (char *) malloc (100 + num_dims*(10 + 4*20)); - sprintf (msg, "Interpolation stencil/molecule out of bounds at grid point [%d", - point[0]); - for (i = 1; i < num_dims; i++) - { - sprintf (msg, "%s, %d", msg, point[i]); - } - sprintf (msg, "%s]\nrange would be min/max [%f / %f", msg, - (double) below[0], (double) (below[0] + offset[0])); - for (i = 1; i < num_dims; i++) - { - sprintf (msg, "%s, %f / %f", msg, - (double) below[i], (double) (below[i] + offset[i])); - } - sprintf (msg, "%s]\ngrid is min/max [%f / %f", msg, - (double) origin[0], (double) (origin[0] + (dims[0]-1)*delta[0])); - for (i = 1; i < num_dims; i++) - { - sprintf (msg, "%s, %f / %f", msg, - (double)origin[i], (double)(origin[i] + (dims[i]-1)*delta[i])); - } - sprintf (msg, "%s]", msg); - CCTK_WARN (1, msg); - free (msg); - - continue; - } -#endif - - /* - * *** compute the interpolation coefficients according to the order *** - * - * (Thanks to Erik for formulating these so nicely.) - * - * These formulas are "just" the coefficients of the classical - * Lagrange interpolation polynomials along each dimension. - * For example, in 1 dimension the unique quadratic passing - * through the 3 points {(x0,y0), (x1,y1), (x2,y2)} is: - * ( x-x1)( x-x2) ( x-x0)( x-x2) ( x-x0)( x-x1) - * -------------- y0 + -------------- y1 + -------------- y2 - * (x0-x1)(x0-x2) (x1-x0)(x1-x2) (x2-x0)(x2-x1) - * (It's easy to see this: each of the terms is yi if x=xi, or - * zero if x=any other xj.) To get the formulas below, just negate - * each (x-x) factor, and substitute the values xi=i. - */ - /* - * NOTE-MAXORDER: support higher interpolation orders by adding the - * appropriate coefficients in another else branch - */ - switch (order) - { - case 1: - /* first order (linear) 1D interpolation */ - for (i = 0; i < num_dims; i++) - { - coeff[i][0] = 1 - offset[i]; - coeff[i][1] = offset[i]; - } - break; - case 2: - /* second order (quadratic) 1D interpolation */ - for (i = 0; i < num_dims; i++) - { - coeff[i][0] = (1-offset[i]) * (2-offset[i]) / ( 2 * 1 ); - coeff[i][1] = ( -offset[i]) * (2-offset[i]) / ( 1 * (-1)); - coeff[i][2] = ( -offset[i]) * (1-offset[i]) / ((-1) * (-2)); - } - break; - case 3: - /* third order (cubic) 1D interpolation */ - for (i = 0; i < num_dims; i++) - { - coeff[i][0] = (1-offset[i]) * (2-offset[i]) * (3-offset[i]) / - ( 3 * 2 * 1 ); - coeff[i][1] = ( -offset[i]) * (2-offset[i]) * (3-offset[i]) / - ( 2 * 1 * (-1)); - coeff[i][2] = ( -offset[i]) * (1-offset[i]) * (3-offset[i]) / - ( 1 * (-1) * (-2)); - coeff[i][3] = ( -offset[i]) * (1-offset[i]) * (2-offset[i]) / - ((-1) * (-2) * (-3)); - } - break; - default: - CCTK_WARN (0, "Implementation error"); - } - -#ifdef PUGHINTERP_VERBOSE_DEBUG -if (n == PUGHInterp_verbose_debug_n) - { - int ii,mm; - for (ii = 0 ; ii < num_dims ; ++ii) - { - for (mm = 0 ; mm <= order ; ++mm) - { - printf("coeff[%d][%d] = %g\n", - ii, mm, (double) coeff[ii][mm]); - } - } - } -#endif /* PUGHINTERP_VERBOSE_DEBUG */ - - /* now loop over all arrays to interpolate at the current point */ - for (a = 0; a < num_arrays; a++) - { - /* sorry, for now input and output arrays must be of same type */ - if (in_types[a] != out_types[a]) - { - CCTK_WARN (1, "Type casting of interpolation results not implemented"); - continue; - } - - /* now do the interpolation according to the array type - we support all kinds of CCTK_REAL* and CCTK_COMPLEX* types here */ - if (in_types[a] == CCTK_VARIABLE_REAL) - { - INTERPOLATE (CCTK_REAL, CCTK_REAL, NOTHING, in_arrays[a], - out_arrays[a], order, point, max_dims, n, coeff); - } - else if (in_types[a] == CCTK_VARIABLE_COMPLEX) - { - INTERPOLATE (CCTK_COMPLEX, CCTK_REAL, .Re, in_arrays[a], - out_arrays[a], order, point, max_dims, n, coeff); - INTERPOLATE (CCTK_COMPLEX, CCTK_REAL, .Im, in_arrays[a], - out_arrays[a], order, point, max_dims, n, coeff); - } -#ifdef CCTK_REAL4 - else if (in_types[a] == CCTK_VARIABLE_REAL4) - { - INTERPOLATE (CCTK_REAL4, CCTK_REAL4, NOTHING, in_arrays[a], - out_arrays[a], order, point, max_dims, n, coeff); - } - else if (in_types[a] == CCTK_VARIABLE_COMPLEX8) - { - INTERPOLATE (CCTK_COMPLEX8, CCTK_REAL4, .Re, in_arrays[a], - out_arrays[a], order, point, max_dims, n, coeff); - INTERPOLATE (CCTK_COMPLEX8, CCTK_REAL4, .Im, in_arrays[a], - out_arrays[a], order, point, max_dims, n, coeff); - } -#endif -#ifdef CCTK_REAL8 - else if (in_types[a] == CCTK_VARIABLE_REAL8) - { - INTERPOLATE (CCTK_REAL8, CCTK_REAL8, NOTHING, in_arrays[a], - out_arrays[a], order, point, max_dims, n, coeff); - } - else if (in_types[a] == CCTK_VARIABLE_COMPLEX16) - { - INTERPOLATE (CCTK_COMPLEX16, CCTK_REAL8, .Re, in_arrays[a], - out_arrays[a], order, point, max_dims, n, coeff); - INTERPOLATE (CCTK_COMPLEX16, CCTK_REAL8, .Im, in_arrays[a], - out_arrays[a], order, point, max_dims, n, coeff); - } -#endif -#ifdef CCTK_REAL16 - else if (in_types[a] == CCTK_VARIABLE_REAL16) - { - INTERPOLATE (CCTK_REAL16, CCTK_REAL16, NOTHING, in_arrays[a], - out_arrays[a], order, point, max_dims, n, coeff); - } - else if (in_types[a] == CCTK_VARIABLE_COMPLEX32) - { - INTERPOLATE (CCTK_COMPLEX32, CCTK_REAL16, .Re, in_arrays[a], - out_arrays[a], order, point, max_dims, n, coeff); - INTERPOLATE (CCTK_COMPLEX32, CCTK_REAL16, .Im, in_arrays[a], - out_arrays[a], order, point, max_dims, n, coeff); - } -#endif - else - { - CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, - "Unsupported variable type %d", (int) in_types[a]); - } - } /* end of loop over all arrays */ - - } /* end of loop over all points to interpolate at */ - - /* we're done */ - return (retval); -} |