From 0e2bd788d3498bc4a674b4d68f0212b905c2c53d Mon Sep 17 00:00:00 2001 From: jthorn Date: Fri, 19 Oct 2001 16:25:13 +0000 Subject: There are two changes in this checkin: * Add a bunch of comments (some grdoc, some "just" ordinary C comments) documenting how the interpolator works, including grdoc comments describing all the arguments of the INTERPOLATE() macro. * #ifdef PUGHINTERP_VERBOSE_DEBUG Add some debugging code to print the interpolation coefficients etc at a single grid point; which grid point is specified by a global variable that the caller can set as appropriate. This is all inside the #ifdef, so in a normal compilation there's no overhead. #endif git-svn-id: http://svn.cactuscode.org/arrangements/CactusPUGH/PUGHInterp/trunk@19 1c20744c-e24a-42ec-9533-f5004cb800e5 --- src/Interpolate.c | 186 ++++++++++++++++++++++++++++++++++++++++++++++++------ 1 file changed, 166 insertions(+), 20 deletions(-) diff --git a/src/Interpolate.c b/src/Interpolate.c index 65e797a..3856ce2 100644 --- a/src/Interpolate.c +++ b/src/Interpolate.c @@ -6,13 +6,18 @@ 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. + 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$ @@*/ @@ -20,6 +25,7 @@ #include #include #include +#include #include "cctk.h" #include "cctk_Parameters.h" @@ -29,6 +35,10 @@ 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 @@ -36,7 +46,84 @@ CCTK_FILEVERSION(CactusPUGH_PUGHInterp_Interpolate_c) /* the highest dimension for variables we can deal with (so far) */ #define MAXDIM 3 -/* macro to do the interpolation of in_array[n] to out_array[n] at point[] */ +/******************************************************************************/ + +/*@@ + @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) \ { \ @@ -49,33 +136,47 @@ CCTK_FILEVERSION(CactusPUGH_PUGHInterp_Interpolate_c) \ interp_result = 0; \ \ - /* NOTE: support >3D arrays by adding more loops */ \ + /* NOTE-MAXDIM: support >3D arrays by adding more loops */ \ for (kk = 0; kk <= order; kk++) \ { \ fk[0] = 0; \ for (jj = 0; jj <= order; jj++) \ { \ - /* NOTE: for >3D arrays adapt the index calculation here */ \ + /* 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++) \ { \ 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[k] */ \ \ /* 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 @@ -83,21 +184,23 @@ CCTK_FILEVERSION(CactusPUGH_PUGHInterp_Interpolate_c) @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 + 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) + - 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 - generic in that it can easily be extended to support higher- - dimensional arrays or more interpolation orders. - Please see the NOTES in this source file for details. + 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 @@ -279,6 +382,18 @@ int PUGHInterp_Interpolate (int order, out_of_bounds |= point[i] < 0 || point[i]+order >= dims[i]; } +#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 */ + /* check bounds */ if (out_of_bounds) { @@ -314,10 +429,26 @@ int PUGHInterp_Interpolate (int order, continue; } - /* compute the interpolation coefficients according to the order - Thanks to Erik for formulating these so nicely. */ - /* NOTE: support higher interpolation orders by adding the - appropriate coefficients in another else branch */ + /* + * *** 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 + */ if (order == 1) { /* first order (linear) 1D interpolation */ @@ -357,6 +488,21 @@ int PUGHInterp_Interpolate (int order, 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++) { -- cgit v1.2.3