aboutsummaryrefslogtreecommitdiff
path: root/src/Interpolate.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/Interpolate.c')
-rw-r--r--src/Interpolate.c609
1 files changed, 609 insertions, 0 deletions
diff --git a/src/Interpolate.c b/src/Interpolate.c
new file mode 100644
index 0000000..6e9af89
--- /dev/null
+++ b/src/Interpolate.c
@@ -0,0 +1,609 @@
+/*@@
+ @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, LOCALINTERP_VERBOSE_DEBUG debugging code
+ @date 22 Jan 2002
+ @author Jonathan Thornburg
+ @hdesc Move all local-interpolation code from LocalInterp to here
+ @endhistory
+
+ @version $Header$
+ @@*/
+
+#include <math.h>
+#include <stdlib.h>
+#include <string.h>
+#include <stdio.h>
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+#include "Interpolate.h"
+
+/* the rcs ID and its dummy function to use it */
+static const char *rcsid = "$Header$";
+CCTK_FILEVERSION(CactusBase_LocalInterp_Interpolate_c)
+
+#ifdef LOCALINTERP_VERBOSE_DEBUG
+ /* if this is >= 0, we print verbose debugging information at this point */
+ int LocalInterp_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 stencil/molecule 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++) \
+ { \
+ fk[0] = 0; \
+ for (jj = 0; jj <= order; jj++) \
+ { \
+ /* 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[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 LocalInterp_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 LocalInterp_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;
+ }
+
+ /* 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];
+ }
+
+#ifdef LOCALINTERP_VERBOSE_DEBUG
+if (n == LocalInterp_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 /* LOCALINTERP_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
+ */
+ if (order == 1)
+ {
+ /* first order (linear) 1D interpolation */
+ for (i = 0; i < num_dims; i++)
+ {
+ coeff[i][0] = 1 - offset[i];
+ coeff[i][1] = offset[i];
+ }
+ }
+ else if (order == 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));
+ }
+ }
+ else if (order == 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));
+ }
+ }
+ else
+ {
+ CCTK_WARN (0, "Implementation error");
+ }
+
+#ifdef LOCALINTERP_VERBOSE_DEBUG
+if (n == LocalInterp_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 /* LOCALINTERP_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);
+}
generated by cgit v1.2.3 (git 2.39.1) at 2024-05-17 23:46:01 +0000