Commented out the check for points where the interpolation stencil is outside

of the grid. This check moved up into the global interpolator in PUGHInterp.
Now, for the interpolation of local arrays, points are allowed to be shifted
back into the grid.

You will also need to update PUGHInterp now to get this check.


git-svn-id: http://svn.cactuscode.org/arrangements/CactusBase/LocalInterp/trunk@121 df1f8a13-aa1d-4dd4-9681-27ded5b42416

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 ***