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

1 files 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 <math.h>
 #include <stdlib.h>
 #include <string.h>
+#include <stdio.h>
 
 #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++)
     {