Moved check for interpolation points with a stencil outside of the grid

up from the local interpolation code into the global interpolator.
This means for the interpolation of local arrays, points are allowed to be
shifted back into the grid.


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

2 files changed, 136 insertions, 43 deletions

diff --git a/src/Interpolate.c b/src/Interpolate.c
index a097abf..146b61c 100644
--- a/src/Interpolate.c
+++ b/src/Interpolate.c
@@ -283,7 +283,10 @@ int PUGHInterp_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];
@@ -359,8 +362,10 @@ int PUGHInterp_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++)
@@ -369,8 +374,10 @@ int PUGHInterp_Interpolate (int order,
       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];
@@ -413,6 +420,7 @@ if (n == PUGHInterp_verbose_debug_n)
 	}
 #endif /* PUGHINTERP_VERBOSE_DEBUG */
 
+#if 0
     /* check bounds */
     if (out_of_bounds)
     {
@@ -447,6 +455,7 @@ if (n == PUGHInterp_verbose_debug_n)
 
       continue;
     }
+#endif
 
     /*
      * *** compute the interpolation coefficients according to the order ***
diff --git a/src/Operator.c b/src/Operator.c
index 8064fd3..0ae230d 100644
--- a/src/Operator.c
+++ b/src/Operator.c
@@ -19,6 +19,7 @@
 
 #include <stdlib.h>
 #include <string.h>
+#include <math.h>    /* floor(3) */
 
 #include "cctk.h"
 #include "CactusPUGH/PUGH/src/include/pugh.h"
@@ -70,6 +71,9 @@ static int CheckArguments (const cGH *GH,
                            int N_input_arrays,
                            int N_output_arrays,
                            const int interp_coord_array_types[]);
+static int CheckOutOfBounds (const cGH *GH, const char *coord_system_name,
+                             int order, int N_dims, int N_points,
+                             const int *dims, const CCTK_REAL *const *coords);
 #ifdef CCTK_MPI
 static int GetLocalCoords (const cGH *GH,
                            int N_points,
@@ -177,26 +181,18 @@ int PUGHInterp_InterpGV (cGH *GH,
                          void *const output_arrays[],
                          const int output_array_types[])
 {
-  const CCTK_REAL **data;
+  int i, nprocs, N_dims, point, array, retval;
   CCTK_REAL *interp_local_coords;
   CCTK_REAL *origin, *delta;
-  int nprocs;
-  int dim, N_dims, *dims;
-  int array;
-  int point;
-  int retval;
+  const CCTK_REAL **data;
   const void **input_arrays;
-  pGH *pughGH;
-  pGExtras *extras;
+  const pGH *pughGH;
+  const pGExtras *extras;
   cGroupDynamicData group_data;
 #ifdef CCTK_MPI
-  int N_local_points;
-  int offset;
-  int proc;
-  int type;
+  int offset, proc, type, maxtype, N_local_points;
   void **local_output_arrays;
   pughInterpGH *myGH;
-  int maxtype;
   type_desc_t *this, *type_desc;
 #endif
 
@@ -219,29 +215,29 @@ int PUGHInterp_InterpGV (cGH *GH,
     return (retval);
   }
 
+  /* get extension handle for PUGH */
+  pughGH = CCTK_GHExtension (GH, "PUGH");
+
+  /* get the extras pointer of the first coordinate
+     This is used later on to verify the layout of the input arrays as well
+     as for mapping points to processors. */
+  i = CCTK_CoordIndex (1, NULL, coord_system_name);
+  extras = ((const pGA *) pughGH->variables[i][0])->extras;
+
   /* get dimensions, origin, and delta of the processor-local grid */
   /* NOTE: getting the dimensions should be a flesh routine as well
            for now we get the dimensions of every coordinate and take the
            i'th element - this is inconsistent !! */
-  dims = malloc (2 * N_dims * sizeof (int));
   origin = malloc (2 * N_dims * sizeof (CCTK_REAL));
   delta  = origin + N_dims;
   input_arrays = malloc (N_input_arrays * sizeof (void *));
-  for (dim = 0; dim < N_dims; dim++)
+  for (i = 0; i < N_dims; i++)
   {
-    CCTK_GrouplshVI (GH, N_dims, dims + N_dims,
-                     CCTK_CoordIndex (dim + 1, NULL, coord_system_name));
-    dims[dim] = dims[dim + N_dims];
-
-    CCTK_CoordLocalRange (GH, &origin[dim], &delta[dim], dim + 1, NULL,
+    CCTK_CoordLocalRange (GH, &origin[i], &delta[i], i + 1, NULL,
                           coord_system_name);
-    delta[dim] = (delta[dim] - origin[dim]) / dims[dim];
+    delta[i] = (delta[i] - origin[i]) / extras->lnsize[i];
   }
 
-  /* get extension handle for PUGH */
-  pughGH = CCTK_GHExtension (GH, "PUGH");
-  extras = NULL;
-
   /* check that the input arrays dimensions match the coordinate system
      (but their dimensionality can be less) */
   for (array = 0; array < N_input_arrays; array++)
@@ -267,7 +263,7 @@ int PUGHInterp_InterpGV (cGH *GH,
       continue;
     }
 
-    if (memcmp (group_data.lsh, dims, group_data.dim * sizeof (int)))
+    if (memcmp (group_data.lsh, extras->lnsize, group_data.dim * sizeof (int)))
     {
       CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
                   "Dimensions of input array variable with index %d "
@@ -276,21 +272,24 @@ int PUGHInterp_InterpGV (cGH *GH,
       retval = -1;
     }
 
-    /* the extras pointer will refer to the first input array with the
-       highest dimension (need this information later on) */
-    if (! extras || extras->dim < group_data.dim)
-    {
-      extras = ((pGA *) pughGH->variables[input_array_indices[array]][0])->extras;
-    }
-
     /* get the data pointer to the input array (use current timelevel) */
     input_arrays[array] = CCTK_VarDataPtrI (GH, 0, input_array_indices[array]);
   }
+  if (retval >= 0)
+  {
+    /* check for out-of-bounds points
+       This check was originally in the local interpolator code but was disabled
+       there and moved up here instead. Now local arrays can have out-of-bounds
+       points, grid arrays cannot. */
+    retval = CheckOutOfBounds (GH, coord_system_name, order, N_dims, N_points,
+                               extras->nsize,
+                               (const CCTK_REAL **) interp_coord_arrays);
+  }
+
   if (retval < 0)
   {
     free (input_arrays);
     free (origin);
-    free (dims);
     return (retval);
   }
 
@@ -304,9 +303,9 @@ int PUGHInterp_InterpGV (cGH *GH,
     data = (const CCTK_REAL **) interp_coord_arrays;
     for (point = 0; point < N_points; point++)
     {
-      for (dim = 0; dim < N_dims; dim++)
+      for (i = 0; i < N_dims; i++)
       {
-        *interp_local_coords++ = data[dim][point];
+        *interp_local_coords++ = data[i][point];
       }
     }
     interp_local_coords -= N_dims * N_points;
@@ -314,7 +313,7 @@ int PUGHInterp_InterpGV (cGH *GH,
     /* call the interpolator function */
     retval = PUGHInterp_Interpolate (order,
                                      N_points, N_dims, N_output_arrays,
-                                     dims, interp_local_coords,
+                                     extras->lnsize, interp_local_coords,
                                      origin, delta,
                                      output_array_types, input_arrays,
                                      output_array_types, output_arrays);
@@ -323,7 +322,6 @@ int PUGHInterp_InterpGV (cGH *GH,
     free (interp_local_coords);
     free (input_arrays);
     free (origin);
-    free (dims);
 
     return (retval);
   }
@@ -405,7 +403,6 @@ int PUGHInterp_InterpGV (cGH *GH,
   if (retval <= 0)
   {
     free (input_arrays);
-    free (dims);
     free (origin);
     free (type_desc);
     return (-1);
@@ -423,7 +420,6 @@ int PUGHInterp_InterpGV (cGH *GH,
   {
     free (input_arrays);
     free (origin);
-    free (dims);
     free (type_desc);
     return (retval);
   }
@@ -478,7 +474,7 @@ int PUGHInterp_InterpGV (cGH *GH,
        output arrays for this processor */
     PUGHInterp_Interpolate (order,
                             myGH->N_points_to[proc], N_dims, N_output_arrays,
-                            dims, interp_local_coords, origin, delta,
+                            extras->lnsize, interp_local_coords, origin, delta,
                             output_array_types, input_arrays,
                             output_array_types, local_output_arrays);
 
@@ -496,7 +492,6 @@ int PUGHInterp_InterpGV (cGH *GH,
   free (local_output_arrays);
   free (input_arrays);
   free (origin);
-  free (dims);
 
   /* now send the interpolation results back to the processors they were
      requested from, also receive my own results that were computed
@@ -1038,3 +1033,92 @@ static int CheckArguments (const cGH *GH,
 
   return (1);
 }
+
+
+static int CheckOutOfBounds (const cGH *GH, const char *coord_system_name,
+                             int order, int N_dims, int N_points,
+                             const int *dims, const CCTK_REAL *const *coords)
+{
+  int i, p, point, out_of_bounds, retval;
+  CCTK_REAL *origin, *delta, *delta_inv, *below, *offset;
+  char *msg;
+
+
+  msg    = malloc (100 + N_dims*(10 + 4*30));
+  origin = malloc (5 * N_dims * sizeof (CCTK_REAL));
+  delta     = origin + 1*N_dims;
+  delta_inv = origin + 2*N_dims;
+  below     = origin + 3*N_dims;
+  offset    = origin + 4*N_dims;
+
+  /* get the global origin and delta of the coordinate system */
+  for (i = 0; i < N_dims; i++)
+  {
+    CCTK_CoordRange (GH, &origin[i], &delta[i], i+1, NULL, coord_system_name);
+    delta[i] = (delta[i] - origin[i]) / (dims[i]-1);
+
+    /* avoid expensive divisions by delta later on */
+    delta_inv[i] = 1.0 / delta[i];
+  }
+
+  retval = 0;
+
+  for (p = 0; p < N_points; p++)
+  {
+    /* reset the out-of-bounds flag */
+    out_of_bounds = 0;
+
+    /* loop over all dimensions */
+    for (i = 0; i < N_dims; i++)
+    {
+      /* closest grid point for stencil */
+      point = floor ((coords[i][p] - origin[i]) * delta_inv[i]
+                     - 0.5 * (order - 1));
+
+      /* test bounds */
+      out_of_bounds |= point < 0 || point+order >= dims[i];
+
+      /* physical coordinate of that grid point */
+      below[i] = origin[i] + point * delta[i];
+
+      /* offset from that grid point, in fractions of grid points */
+      offset[i] = (coords[i][p] - below[i]) * delta_inv[i];
+    }
+
+    /* check bounds */
+    if (out_of_bounds)
+    {
+      /* put all information into a single message string for output */
+      sprintf (msg, "Interpolation stencil out of bounds at grid point [%f",
+               (double) coords[0][p]);
+      for (i = 1; i < N_dims; i++)
+      {
+        sprintf (msg, "%s, %f", msg, (double) coords[i][p]);
+      }
+      sprintf (msg, "%s]\nrange would be min/max [%f / %f", msg,
+               (double) below[0], (double) (below[0] + offset[0]));
+      for (i = 1; i < N_dims; i++)
+      {
+        sprintf (msg, "%s, %f / %f", msg,
+                 (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]-1)*delta[0]));
+      for (i = 1; i < N_dims; i++)
+      {
+        sprintf (msg, "%s, %f / %f", msg,
+                 (double)origin[i], (double)(origin[i] + (dims[i]-1)*delta[i]));
+      }
+      sprintf (msg, "%s]", msg);
+      CCTK_WARN (1, msg);
+
+      retval--;
+    }
+  }
+
+  /* free allocated resources */
+  free (origin);
+  free (msg);
+
+  return (retval);
+}