From cb9aaede7365d60bc7e0543ce36873de7893607a Mon Sep 17 00:00:00 2001
From: tradke <tradke@df1f8a13-aa1d-4dd4-9681-27ded5b42416>
Date: Thu, 19 Dec 2002 11:11:52 +0000
Subject: 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
---
 src/UniformCartesian/Interpolate.c | 32 ++++++++++++++++++++++++--------
 1 file changed, 24 insertions(+), 8 deletions(-)

(limited to 'src')

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 ***
-- 
cgit v1.2.3