Set the N_boundary_points_to_omit[], off_centering_tolerance[], and

extrapolation_tolerance[] arrays in the options table to be passed down to
the local interpolator.


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

1 files changed, 147 insertions, 17 deletions

diff --git a/src/InterpGridArrays.c b/src/InterpGridArrays.c
index 369510f..d9e9abe 100644
--- a/src/InterpGridArrays.c
+++ b/src/InterpGridArrays.c
@@ -78,8 +78,9 @@ typedef struct
 /********************************************************************
  ********************    Internal Routines   ************************
  ********************************************************************/
-static int EvaluateParameterTable (int param_table_handle, int N_input_arrays,
-                                   CCTK_INT *input_array_time_levels);
+static int PrepareParameterTable (const cGH *GH, int param_table_handle,
+                                  int N_dims, int N_input_arrays,
+                                  CCTK_INT *input_array_time_levels);
 
 
 /*@@
@@ -271,8 +272,9 @@ int PUGHInterp_InterpGridArrays (const cGH *GH,
   input_array_time_levels = input_array_dims + N_dims + N_input_arrays;
 
   /* evaluate the options from the user-supplied parameter table */
-  retval = EvaluateParameterTable (param_table_handle, N_input_arrays,
-                                   input_array_time_levels);
+  retval = PrepareParameterTable (GH, param_table_handle,
+                                  N_dims, N_input_arrays,
+                                  input_array_time_levels);
   if (retval < 0)
   {
     free (origin_local);
@@ -352,8 +354,28 @@ int PUGHInterp_InterpGridArrays (const cGH *GH,
     input_array_types[i] = CCTK_VarTypeI (input_array_indices[i]);
   }
 
+#if 0
+  /* Check whether the number of ghostzones is sufficient for the local
+     interpolator to handle points on the processor's local bounding box.
+     This is done as a query, simply by setting N_output_arrays = 0. */
+  if (retval == 0)
+  {
+    for (i = 0; i < N_dims; i++)
+    {
+      bbox_interp_coords[i] = origin_local;
+    retval = CCTK_InterpLocalUniform (N_dims, local_interp_handle,
+                                      param_table_handle,
+                                      origin_local, delta_local, 2,
+                                      interp_coords_type, bbox_interp_coords,
+                                      N_input_arrays, input_array_dims,
+                                      input_array_types, input_arrays,
+                                      0, output_array_types,
+                                      output_arrays);
+  }
+#endif
+
   /* single-processor case directly calls the local interpolator */
-  if (retval >= 0 && pughGH->nprocs == 1)
+  if (retval == 0 && pughGH->nprocs == 1)
   {
     retval = CCTK_InterpLocalUniform (N_dims, local_interp_handle,
                                       param_table_handle,
@@ -899,18 +921,30 @@ int PUGHInterp_InterpGridArrays (const cGH *GH,
  ********************    Internal Routines   ************************
  ********************************************************************/
  /*@@
-   @routine    EvaluateParameterTable
+   @routine    PrepareParameterTable
    @date       Fri 20 Dec 2002
    @author     Thomas Radke
    @desc
-               Parses the options given in the user-supplied parameter table.
+               Parses the 'input_array_time_levels' option from the
+               user-supplied parameter table and sets up the bounding box
+               option arrays for the local interpolator
    @enddesc
 
+   @var        GH
+   @vdesc      pointer to the CCTK grid hierarchy
+   @vtype      const cGH *
+   @vio        in
+   @endvar
    @var        param_table_handle
    @vdesc      handle to the user-supplied parameter table
    @vtype      int
    @vio        in
    @endvar
+   @var        N_dims
+   @vdesc      number of dimensions in which to interpolate
+   @vtype      int
+   @vio        in
+   @endvar
    @var        N_input_arrays
    @vdesc      total number of input arrays
    @vtype      int
@@ -927,11 +961,16 @@ int PUGHInterp_InterpGridArrays (const cGH *GH,
                0 for success, or one of the UTIL_ERROR_TABLE_* error codes
    @endreturndesc
 @@*/
-static int EvaluateParameterTable (int param_table_handle, int N_input_arrays,
-                                   CCTK_INT *input_array_time_levels)
+static int PrepareParameterTable (const cGH *GH, int param_table_handle,
+                                  int N_dims, int N_input_arrays,
+                                  CCTK_INT *input_array_time_levels)
 {
-  int N_elements, retval;
+  int i, retval;
+  CCTK_INT *N_boundary_points_to_omit;
+  CCTK_REAL *off_centering_tolerance, *extrapolation_tolerance;
+
 
+  retval = 0;
 
   if (param_table_handle < 0)
   {
@@ -939,27 +978,118 @@ static int EvaluateParameterTable (int param_table_handle, int N_input_arrays,
     return (UTIL_ERROR_BAD_HANDLE);
   }
 
-  retval = 0;
-  N_elements = Util_TableGetIntArray (param_table_handle, N_input_arrays,
-                                      input_array_time_levels,
-                                      "input_array_time_levels");
-  if (N_elements == UTIL_ERROR_TABLE_NO_SUCH_KEY)
+  /* allocate table option arrays and initialize them to their defaults */
+  N_boundary_points_to_omit = malloc (2 * N_dims * sizeof (CCTK_INT));
+  off_centering_tolerance = malloc (2 * 2 * N_dims * sizeof (CCTK_REAL));
+  extrapolation_tolerance = off_centering_tolerance + 2 * N_dims;
+  for (i = 0; i < 2 * N_dims; i++)
+  {
+    N_boundary_points_to_omit[i] = 0;
+    off_centering_tolerance[i] = 0.0;
+    extrapolation_tolerance[i] = 1e-10;
+  }
+
+  /* read the 'N_boundary_points_to_omit' options array from the table
+     and verify its datatype and element count */
+  i = Util_TableGetIntArray (param_table_handle, 2 * N_dims,
+                             N_boundary_points_to_omit,
+                             "N_boundary_points_to_omit");
+  if (i < 0 && i != UTIL_ERROR_TABLE_NO_SUCH_KEY)
+  {
+    CCTK_WARN (1, "Array with key 'N_boundary_points_to_omit' must be of "
+                  "CCTK_INT datatype");
+    retval = UTIL_ERROR_BAD_INPUT;
+  }
+  else if (i >= 0 && i < 2 * N_dims)
+  {
+    CCTK_WARN (1, "Array with key 'N_boundary_points_to_omit' has fewer "
+                  "elements than interpolation bounding box");
+    retval = UTIL_ERROR_BAD_INPUT;
+  }
+
+  /* read the 'off_centering_tolerance' options array from the table 
+     and verify its datatype and element count */
+  i = Util_TableGetRealArray (param_table_handle, 2 * N_dims,
+                              off_centering_tolerance,
+                              "off_centering_tolerance");
+  if (i < 0 && i != UTIL_ERROR_TABLE_NO_SUCH_KEY)
+  {
+    CCTK_WARN (1, "Array with key 'off_centering_tolerance' must be of "
+                  "CCTK_REAL datatype");
+    retval = UTIL_ERROR_BAD_INPUT;
+  }
+  else if (i >= 0 && i < 2 * N_dims)
+  {
+    CCTK_WARN (1, "Array with key 'off_centering_tolerance' has fewer elements "
+                  "than interpolation bounding box");
+    retval = UTIL_ERROR_BAD_INPUT;
+  }
+
+  /* read the 'extrapolation_tolerance' options array from the table
+     and verify its datatype and element count */
+  i = Util_TableGetRealArray (param_table_handle, 2 * N_dims,
+                              extrapolation_tolerance, 
+                              "extrapolation_tolerance");
+  if (i < 0 && i != UTIL_ERROR_TABLE_NO_SUCH_KEY)
+  {
+    CCTK_WARN (1, "Array with key 'extrapolation_tolerance' must be of "
+                  "CCTK_REAL datatype");
+    retval = UTIL_ERROR_BAD_INPUT;
+  }
+  else if (i >= 0 && i < 2 * N_dims)
+  {
+    CCTK_WARN (1, "Array with key 'extrapolation_tolerance' has fewer elements "
+                  "than interpolation bounding box");
+    retval = UTIL_ERROR_BAD_INPUT;
+  }
+
+  /* read the 'input_array_time_levels' options array from the table 
+     and verify its datatype and element count */
+  i = Util_TableGetIntArray (param_table_handle, N_input_arrays,
+                             input_array_time_levels,
+                             "input_array_time_levels");
+  if (i == UTIL_ERROR_TABLE_NO_SUCH_KEY)
   {
     /* if no such option was specified, take the current timelevel as default */
     memset (input_array_time_levels, 0, N_input_arrays * sizeof (CCTK_INT));
   }
-  else if (N_elements < 0)
+  else if (i < 0)
   {
     CCTK_WARN (1, "Array with key 'input_array_time_levels' must be of "
                   "CCTK_INT datatype");
     retval = UTIL_ERROR_BAD_INPUT;
   }
-  else if (N_elements < N_input_arrays)
+  else if (i < N_input_arrays)
   {
     CCTK_WARN (1, "Array with key 'input_array_time_levels' has fewer elements "
                   "than total number of input arrays");
     retval = UTIL_ERROR_BAD_INPUT;
   }
 
+  /* complete the bounding box arrays and set them in the parameter table */
+  for (i = 0; i < 2 * N_dims; i++)
+  {
+    if (! GH->cctk_bbox[i])
+    {
+      N_boundary_points_to_omit[i] = 0;
+      off_centering_tolerance[i] = 0.0;
+      extrapolation_tolerance[i] = 0.0;
+    }
+  }
+  if (! retval)
+  {
+    Util_TableSetIntArray (param_table_handle, 2 * N_dims,
+                           N_boundary_points_to_omit,
+                           "N_boundary_points_to_omit");
+    Util_TableSetRealArray (param_table_handle, 2 * N_dims,
+                            off_centering_tolerance, "off_centering_tolerance");
+    Util_TableSetRealArray (param_table_handle, 2 * N_dims,
+                            extrapolation_tolerance, "extrapolation_tolerance");
+  }
+
+  /* free allocated resources */
+  free (N_boundary_points_to_omit);
+  free (off_centering_tolerance);
+
   return (retval);
 }