diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/UniformCartesian/Interpolate.c | 32 |
1 files changed, 24 insertions, 8 deletions
diff --git a/src/UniformCartesian/Interpolate.c b/src/UniformCartesian/Interpolate.c index 95b44ce..52b8fce 100644 --- a/src/UniformCartesian/Interpolate.c +++ b/src/UniformCartesian/Interpolate.c @@ -283,8 +283,12 @@ int LocalInterp_Interpolate (int order, void *const out_arrays[]) { int retval; - int i, a, n, out_of_bounds, shift; + 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]; @@ -331,6 +335,12 @@ int LocalInterp_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 */ @@ -351,16 +361,23 @@ int LocalInterp_Interpolate (int order, /* 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 */ - 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)); +#if 0 + /* test bounds */ + out_of_bounds |= point[i] < 0 || point[i]+order >= dims[i]; +#endif + /* if beyond lower bound shift the grid point to the right */ shift = point[i]; if (shift < 0) @@ -379,10 +396,7 @@ int LocalInterp_Interpolate (int order, 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]; + offset[i] = (coord[num_dims*n + i] - below[i]) * delta_inv[i]; } #ifdef LOCALINTERP_VERBOSE_DEBUG @@ -397,6 +411,7 @@ if (n == LocalInterp_verbose_debug_n) } #endif /* LOCALINTERP_VERBOSE_DEBUG */ +#if 0 /* check bounds */ if (out_of_bounds) { @@ -419,11 +434,11 @@ if (n == LocalInterp_verbose_debug_n) (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])); + (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]*delta[i])); + (double)origin[i], (double)(origin[i] + (dims[i]-1)*delta[i])); } sprintf (msg, "%s]", msg); CCTK_WARN (1, msg); @@ -431,6 +446,7 @@ if (n == LocalInterp_verbose_debug_n) continue; } +#endif /* * *** compute the interpolation coefficients according to the order *** |