Fixed IOUtil_1DLines() and IOUtil_2DPlanes() to setup slice centers even if

no coordinate system is available but slice index coordinates were specified
by the user instead.
This closes PR CactusBase/1536.


git-svn-id: http://svn.cactuscode.org/arrangements/CactusBase/IOUtil/trunk@190 b32723a9-ab3a-4a60-88e2-2e5d99d7c17a

1 files changed, 124 insertions, 86 deletions

diff --git a/src/Utils.c b/src/Utils.c
index 4b559f6..4870d61 100644
--- a/src/Utils.c
+++ b/src/Utils.c
@@ -301,105 +301,122 @@ int IOUtil_1DLines (const cGH *GH,
                     CCTK_REAL *const *const origin_phys,
                     int *const *slice_center)
 {
-  int dim, dir;
+  int dim, dir, handle, no_origin_index;
   char coord_system_name[20];
   CCTK_REAL *lower_range, *upper_range;
-  int retval;
 
-  retval = 0;
-
-  /* allocate arrays for ranges */
-  lower_range = (CCTK_REAL *) calloc (2 * num_dims, sizeof (CCTK_REAL));
-  upper_range = lower_range + num_dims;
 
   /* get the appropriate coordinate system name */
   sprintf (coord_system_name, "cart%dd", num_dims);
-  if (CCTK_CoordSystemHandle (coord_system_name) >= 0)
+  handle = CCTK_CoordSystemHandle (coord_system_name);
+
+  /* if no coordinate system is available then origin_index[] must be given */
+  no_origin_index = origin_index == NULL;
+  if (handle < 0 && origin_index)
   {
-    /* get the ranges in every direction */
+    /* check that origin_index[] is valid */
     for (dir = 0; dir < num_dims; dir++)
     {
-      if (CCTK_CoordRange (GH, &lower_range[dir], &upper_range[dir],
-                           dir + 1, NULL, coord_system_name) < 0)
+      for (dim = 0; dim < num_dims; dim++)
       {
-        CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
-                    "IOUtil_1DLines: Could not get ranges for %c-direction "
-                    "of coordinate system '%s'",
-                    'x' + dir, coord_system_name);
+        if (origin_index[dir][dim] < 0)
+        {
+          no_origin_index = 1;
+          break;
+        }
       }
     }
   }
-  else
+
+  if (handle < 0 && no_origin_index)
   {
-    CCTK_VWarn (4, __LINE__, __FILE__, CCTK_THORNSTRING,
-                "IOUtil_1DLines: Cartesian coordinate system '%s' not found"
-                " - no slice center set up for output of 1D lines from %dD "
-                "variables", coord_system_name, num_dims);
+    CCTK_VWarn (3, __LINE__, __FILE__, CCTK_THORNSTRING,
+                "IOUtil_1DLines: Cartesian coordinate system '%s' not found "
+                "and no slice center index coordinates set - slice center will "
+                "not be set up for output of 1D lines from %dD variables",
+                coord_system_name, num_dims);
     for (dir = 0; dir < num_dims; dir++)
     {
-      for (dim = 0; dim < num_dims; dim++)
+      memset (slice_center[dir], 0, num_dims * sizeof (int));
+    }
+
+    return (-1);
+  }
+
+  lower_range = upper_range = NULL;
+  if (handle >= 0)
+  {
+    /* allocate arrays for ranges */
+    lower_range = calloc (2 * num_dims, sizeof (CCTK_REAL));
+    upper_range = lower_range + num_dims;
+
+    /* get the ranges in every direction */
+    for (dir = 0; dir < num_dims; dir++)
+    {
+      if (CCTK_CoordRange (GH, &lower_range[dir], &upper_range[dir],
+                           dir + 1, NULL, coord_system_name) < 0)
       {
-        slice_center[dir][dim] = 0;
+        CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+                    "IOUtil_1DLines: Could not get ranges for %c-direction "
+                    "of coordinate system '%s'",
+                    'x' + dir, coord_system_name);
       }
     }
-
-    retval = -1;
   }
 
   /* now set the slice center for each line
      according to origin_index[] or origin_phys[] */
-
-  if (! retval)
+  for (dir = 0; dir < num_dims; dir++)
   {
-    for (dir = 0; dir < num_dims; dir++)
+    for (dim = 0; dim < num_dims; dim++)
     {
-      for (dim = 0; dim < num_dims; dim++)
+      if (dim == dir)
       {
-        if (dim == dir)
-        {
-          /* line always starts at the first point */
-          slice_center[dir][dim] = 0;
-        }
-        else if (origin_index && origin_index[dir][dim] >= 0)
-        {
-          /* FIXME: check upper index bounds also ?? */
-          slice_center[dir][dim] = origin_index[dir][dim];
-        }
-        else if (lower_range[dim] > origin_phys[dir][dim] ||
-                 upper_range[dim] < origin_phys[dir][dim])
-        {
-          CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
-                      "IOUtil_1DLines: %c-coordinate (%f) for slice center of "
-                      "1D lines in %c-direction for %dD variables is out of "
-                      "grid coordinates range (%f, %f)",
-                      'x' + dim, (double) origin_phys[dir][dim],
-                      'x' + dir, num_dims,
-                      (double) lower_range[dim], (double) upper_range[dim]);
-          CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
-                      "IOUtil_1DLines: slice center will default to %c-index 0",
-                      'x' + dir);
-          slice_center[dir][dim] = 0;
-        }
-        else
-        {
-          /* Find index for first point above the chosen coordinate */
-          slice_center[dir][dim] =
-            ceil ((origin_phys[dir][dim] - lower_range[dim]) /
-                  GH->cctk_delta_space[dim] - 1e-6);
+        /* line always starts at the first point */
+        slice_center[dir][dim] = 0;
+      }
+      else if (origin_index && origin_index[dir][dim] >= 0)
+      {
+        /* FIXME: check upper index bounds also ?? */
+        slice_center[dir][dim] = origin_index[dir][dim];
+      }
+      else if (lower_range[dim] > origin_phys[dir][dim] ||
+               upper_range[dim] < origin_phys[dir][dim])
+      {
+        CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+                    "IOUtil_1DLines: %c-coordinate (%f) for slice center of "
+                    "1D lines in %c-direction for %dD variables is out of "
+                    "grid coordinates range (%f, %f)",
+                    'x' + dim, (double) origin_phys[dir][dim],
+                    'x' + dir, num_dims,
+                    (double) lower_range[dim], (double) upper_range[dim]);
+        CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+                    "IOUtil_1DLines: slice center will default to %c-index 0",
+                    'x' + dir);
+        slice_center[dir][dim] = 0;
+      }
+      else
+      {
+        /* Find index for first point above the chosen coordinate */
+        slice_center[dir][dim] =
+          ceil ((origin_phys[dir][dim] - lower_range[dim]) /
+                GH->cctk_delta_space[dim] - 1e-6);
 
 #ifdef DEBUG_IOUTIL
-          printf("spxyz for %c-coord of lines in %c-direction is %d\n",
-                 'x' + dim,'x' + dir, slice_center[dir][dim]);
+        printf("spxyz for %c-coord of lines in %c-direction is %d\n",
+               'x' + dim,'x' + dir, slice_center[dir][dim]);
 #endif
-        }
       }
     }
   }
 
   /* free allocated resources */
-  free (lower_range);
+  if (lower_range)
+  {
+    free (lower_range);
+  }
 
-  return (retval);
+  return (0);
 }
 
 
@@ -452,37 +469,57 @@ int IOUtil_2DPlanes (const cGH *GH,
                      const CCTK_REAL *origin_phys,
                      int *slice_center)
 {
-  int dir;
+  int handle, dir;
   char coord_system_name[20];
   CCTK_REAL *lower_range, *upper_range;
 
 
   /* get the appropriate coordinate system name */
   sprintf (coord_system_name, "cart%dd", num_dims);
+  handle = CCTK_CoordSystemHandle (coord_system_name);
 
-  if (CCTK_CoordSystemHandle (coord_system_name) < 0)
+  /* if no coordinate system is available then origin_index[] must be given */
+  if (handle < 0 && origin_index)
   {
-    CCTK_VWarn (4, __LINE__, __FILE__, CCTK_THORNSTRING,
-                "IOUtil_2DPlanes: Cartesian coordinate system '%s' not found"
-                " - no slice center set up for output of 2D planes from %dD "
-                "variables", coord_system_name, num_dims);
-    return (-1);
+    /* check that origin_index[] is valid */
+    for (dir = 0; dir < num_dims; dir++)
+    {
+      if (origin_index[dir] < 0)
+      {
+        origin_index = NULL;
+        break;
+      }
+    }
   }
 
-  /* allocate arrays for ranges */
-  lower_range = (CCTK_REAL *) calloc (2 * num_dims, sizeof (CCTK_REAL));
-  upper_range = lower_range + num_dims;
+  if (handle < 0 && ! origin_index)
+  {
+    CCTK_VWarn (3, __LINE__, __FILE__, CCTK_THORNSTRING,
+                "IOUtil_2DPlanes: Cartesian coordinate system '%s' not found "
+                "and no slice center index coordinates set - slice center will "
+                "not be set up for output of 2D planes from %dD variables",
+                coord_system_name, num_dims);
+    return (-1);
+  }
 
-  /* get the ranges in every direction */
-  for (dir = 0; dir < num_dims; dir++)
+  lower_range = upper_range = NULL;
+  if (handle >= 0)
   {
-    if (CCTK_CoordRange (GH, &lower_range[dir], &upper_range[dir],
-                         num_dims - dir, NULL, coord_system_name) < 0)
+    /* allocate arrays for ranges */
+    lower_range = calloc (2 * num_dims, sizeof (CCTK_REAL));
+    upper_range = lower_range + num_dims;
+
+    /* get the ranges in every direction */
+    for (dir = 0; dir < num_dims; dir++)
     {
-      CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
-                  "IOUtil_2DPlanes: Could not get ranges for %c-direction "
-                  "of coordinate system '%s'",
-                  'x' + (num_dims - dir - 1), coord_system_name);
+      if (CCTK_CoordRange (GH, &lower_range[dir], &upper_range[dir],
+                           num_dims - dir, NULL, coord_system_name) < 0)
+      {
+        CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+                    "IOUtil_2DPlanes: Could not get ranges for %c-direction "
+                    "of coordinate system '%s'",
+                    'x' + (num_dims - dir - 1), coord_system_name);
+      }
     }
   }
 
@@ -522,7 +559,10 @@ int IOUtil_2DPlanes (const cGH *GH,
   }
 
   /* free allocated resources */
-  free (lower_range);
+  if (lower_range)
+  {
+    free (lower_range);
+  }
 
   return (0);
 }
@@ -704,13 +744,11 @@ static void SetOutputVar (int vindex, const char *optstring, void *arg)
   char key[128];
   CCTK_INT type, nelems;
   char *fullname;
-  const info_t *info;
+  const info_t *info = arg;
   ioRequest *request;
   DECLARE_CCTK_PARAMETERS
 
 
-  info = (const info_t *) arg;
-
   /* allocate a new I/O request structure and initialize it with defaults */
   request = IOUtil_DefaultIORequest (info->GH, vindex, info->out_every_default);