/*@@ @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. @enddesc @history @date Wed 17 Jan 2001 @author Thomas Radke @hdesc Translation from Fortran to C @endhistory @version $Id$ @@*/ #include #include #include #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) /* 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 /* macro to do the interpolation of in_array[n] to out_array[n] at point[] */ #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: 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 */ \ 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]; \ } \ fk[0] += fj[0] * coeff[1][jj]; \ } \ interp_result += fk[0] * coeff[2][kk]; \ } \ \ /* assign the result */ \ ((cctk_type *) out_array)[n]subpart = interp_result; \ } /* 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 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. @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, out_of_bounds, shift; int max_dims[MAXDIM], point[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); } /* 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++) { /* reset the out-of-bounds flag */ out_of_bounds = 0; /* loop over all dimensions */ for (i = 0; i < num_dims; i++) { /* closest grid point for stencil */ point[i] = floor ((coord[num_dims*n + i] - origin[i]) / delta[i] - 0.5 * (order - 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[i]; /* test bounds */ out_of_bounds |= point[i] < 0 || point[i]+order >= dims[i]; } /* 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 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]*delta[0])); for (i = 1; i < num_dims; i++) { sprintf (msg, "%s, %f / %f", msg, (double) origin[i], (double) (origin[i] + dims[i]*delta[i])); } sprintf (msg, "%s]", msg); CCTK_WARN (1, msg); free (msg); 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 */ 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"); } /* 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); }