Erik's patch which will

  (a) fix GetLocalCoords to check for too few ghostzones, but with a kludge,
      because I don't know how to get other processors' bbox information from PUGH.
  (b) change the interpolator to not shift the stencil any more.
This closes PR CactusPUGH/1202.


git-svn-id: http://svn.cactuscode.org/arrangements/CactusPUGH/PUGHInterp/trunk@31 1c20744c-e24a-42ec-9533-f5004cb800e5

2 files changed, 64 insertions, 80 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
diff --git a/src/Operator.c b/src/Operator.c
index 90a15be..fc6e2cf 100644
--- a/src/Operator.c
+++ b/src/Operator.c
@@ -771,44 +771,7 @@ static int GetLocalCoords (cGH *GH,
   CCTK_REAL *range_min, *range_max;
   CCTK_REAL *origin, *delta;
   CCTK_REAL *proc_coords;
-#define FUDGE 0.00001
-
-
-#if 0
-  /* TR 03/05/01
-     Following flipperoni code was commented out because it doesn't seem
-     to be used anywhere.
-     By this we can also get rid of USING the grid::domain parameter.
-   */
-
-  /* Do a fliperooni in octant mode */
-  /* This flips the sign around the origin only */
-  /* if the left side (cx0) is bigger than the position to interpolate at */
-  /* BoxInBox needs this for interpolating onto grids */ 
-  if (CCTK_Equals (domain, "octant"))
-  {
-    for (point = 0; point < num_points; point++)
-    {
-      tmp = 0;
-      for (dim = 0; dim < extras->dim; dim++)
-      {
-        if (coords[dim][point] <= GH->cctk_origin_space[dim])
-        {
-          tmp = 1;                 /* mark flipping */
-          break;
-        }
-      }
-/*** FIXME *** really change input arrays here ??? */
-      if (tmp)
-      {
-        for (dim = 0; dim < extras->dim; dim++)
-        {
-          coords[dim][point] = fabs (coords[dim][point]);
-        }
-      }
-    }
-  }
-#endif
+#define FUDGE 0.0
 
   /* get GH extension handles for PUGHInterp and PUGH */
   myGH = (pughInterpGH *) CCTK_GHExtension (GH, "PUGHInterp");
@@ -856,8 +819,12 @@ static int GetLocalCoords (cGH *GH,
     for (dim = 0; dim < extras->dim; dim++)
     {
       /* compute the coordinate ranges */
-      range_min[dim] = origin[dim] + (extras->lb[proc][dim] - FUDGE)*delta[dim];
-      range_max[dim] = origin[dim] + (extras->ub[proc][dim] + FUDGE)*delta[dim];
+      /* TODO: use bbox instead -- but the bboxes of other processors
+	 are now known */
+      int const has_lower = extras->lb[proc][dim] == 0;
+      int const has_upper = extras->ub[proc][dim] == GH->cctk_gsh[dim]-1;
+      range_min[dim] = origin[dim] + (extras->lb[proc][dim] + (!has_lower) * (extras->nghostzones[dim]-0.5) - FUDGE)*delta[dim];
+      range_max[dim] = origin[dim] + (extras->ub[proc][dim] - (!has_upper) * (extras->nghostzones[dim]-0.5) + FUDGE)*delta[dim];
     }
 
     /* and now which point will be processed by what processor */