diff options
Diffstat (limited to 'src/Interpolate.c')
-rw-r--r-- | src/Interpolate.c | 97 |
1 files changed, 57 insertions, 40 deletions
diff --git a/src/Interpolate.c b/src/Interpolate.c index 0ddd1b3..025784c 100644 --- a/src/Interpolate.c +++ b/src/Interpolate.c @@ -22,6 +22,7 @@ @version $Id$ @@*/ +#include <assert.h> #include <math.h> #include <stdlib.h> #include <string.h> @@ -138,9 +139,11 @@ CCTK_FILEVERSION(CactusPUGH_PUGHInterp_Interpolate_c) \ /* 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 + \ @@ -148,6 +151,7 @@ CCTK_FILEVERSION(CactusPUGH_PUGHInterp_Interpolate_c) \ fj[0] = 0; \ for (ii = 0; ii <= order; ii++) \ + /*if (fabs(coeff[0][ii]) > 1e-6)*/ \ { \ fj[0] += fi[ii]subpart * coeff[0][ii]; \ } \ @@ -281,6 +285,7 @@ int PUGHInterp_Interpolate (int order, int retval; int i, a, n, out_of_bounds, shift; 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]; @@ -328,6 +333,12 @@ int PUGHInterp_Interpolate (int order, 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 */ @@ -355,9 +366,12 @@ int PUGHInterp_Interpolate (int order, for (i = 0; i < num_dims; i++) { /* closest grid point for stencil */ - point[i] = floor ((coord[num_dims*n + i] - origin[i]) / delta[i] + point[i] = floor ((coord[num_dims*n + i] - origin[i]) * delta_inv[i] - 0.5 * (order - 1)); + /* test bounds */ + out_of_bounds |= point[i] < 0 || point[i]+order >= dims[i]; + /* if beyond lower bound shift the grid point to the right */ shift = point[i]; if (shift < 0) @@ -371,15 +385,20 @@ int PUGHInterp_Interpolate (int order, { 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[i]; + offset[i] = (coord[num_dims*n + i] - below[i]) * delta_inv[i]; - /* test bounds */ - out_of_bounds |= point[i] < 0 || point[i]+order >= dims[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 @@ -449,43 +468,41 @@ if (n == PUGHInterp_verbose_debug_n) * NOTE-MAXORDER: support higher interpolation orders by adding the * appropriate coefficients in another else branch */ - if (order == 1) + switch (order) { - /* 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"); + 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 |