Use ceil() when computing the slice centers for 1D lines/2D planes.

This should fix PR CactusBase/536.

Also added missing grdoc and removed source file ChooseOutput.c
which isn't used here.


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

1 files changed, 206 insertions, 189 deletions

diff --git a/src/Utils.c b/src/Utils.c
index 338bed9..b965d3b 100644
--- a/src/Utils.c
+++ b/src/Utils.c
@@ -1,267 +1,284 @@
- /*@@
+/*@@
    @file      Utils.c
    @date      Tue 4th July 2000
    @author    Gabrielle Allen
-   @desc 
-   Utility routines which may be called by other IO thorns
+   @desc
+              Utility routines which may be called by other IO thorns
    @enddesc 
-   @history
-   @endhistory
+   @version   $Id$
  @@*/
 
-/*#define DEBUG_IOUTIL 1*/
+
+#include <math.h>
+#include <stdlib.h>
 
 #include "cctk.h"
 #include "cctk_Parameters.h"
 #include "ioutil_Utils.h"
+#include "ioutil_CheckpointRecovery.h"
 
+/* the rcs ID and its dummy function to use it */
+static char *rcsid = "$Header$";
+CCTK_FILEVERSION(CactusBase_IOUtil_Utils_c)
 
- /*@@
+
+/* uncomment this if you want some debugging output */
+/* #define DEBUG_IOUTIL 1 */
+
+
+/*@@
    @routine    IOUtil_1DLines
    @date       July 4 2000
    @author     Gabrielle Allen, Gerd Lanfermann, Thomas Radke
-   @desc 
-   Fills out an array determining where to position 1D lines for output 
-   on a multidimensional regular Cartesian unigrid.
-   The first slot of the array specifies the 1D line direction, the
-   second slot fixes the indices of the starting point of that line on the
-   grid.
-   @enddesc 
-   @calls     
-   @calledby   
-   @history 
- 
-   @endhistory 
+   @desc
+               Fills out an array determining where to position 1D lines
+               for output on a multidimensional regular Cartesian unigrid.
+               The first slot of the array specifies the 1D line direction,
+               the second slot fixes the indices of the starting point
+               of that line on the grid.
+   @enddesc
 
-@@*/
+   @calls      CCTK_CoordRange
 
-int IOUtil_1DLines (cGH *GH, 
-                    int dimension, 
-                    int **index, 
-                    CCTK_REAL **coord, 
-                    int **spxyz)
+   @var        GH
+   @vdesc      Pointer to CCTK GH
+   @vtype      cGH *
+   @vio        in
+   @endvar
+   @var        num_dims
+   @vdesc      number of dimensions of underlying grid
+   @vtype      int
+   @vio        in
+   @endvar
+   @var        origin_index
+   @vdesc      origin for 1D lines given by grid point indices
+   @vtype      const int [num_dims][num_dims]
+   @vio        in
+   @endvar
+   @var        origin_phys
+   @vdesc      origin for 1D lines given by physical coordinates
+   @vtype      const CCTK_REAL [num_dims][num_dims]
+   @vio        in
+   @endvar
+   @var        slice_center
+   @vdesc      resulting 1D line slice center
+   @vtype      int [num_dims][num_dims]
+   @vio        out
+   @endvar
+
+   @returntype int
+   @returndesc
+               0  for success
+               -1 no coordinate system of given dimensions found
+   @endreturndesc
+ @@*/
+int IOUtil_1DLines (cGH *GH,
+                    int num_dims,
+                    int *const origin_index[],
+                    CCTK_REAL *const origin_phys[],
+                    int *slice_center[])
 {
   DECLARE_CCTK_PARAMETERS
-  int ierr=0;
   int dim, dir;
-  CCTK_REAL lower,upper;
-  const char *name=NULL;
-  CCTK_REAL lowercoords[3]={0,0,0};
-  CCTK_REAL uppercoords[3]={0,0,0};
+  char system_name[20];
+  CCTK_REAL *lower_range, *upper_range;
+
 
+  /* get the appropriate coordinate system name */
+  sprintf (system_name, "cart%dd", num_dims);
 
-  /* Get the appropriate coordinate system */
-  switch (dimension)
+  if (CCTK_CoordSystemHandle (system_name) < 0)
   {
-    case 1:
-      name = "cart1d";
-      ierr = CCTK_CoordRange(GH, &lower, &upper, 1, NULL, name);
-      if (ierr>-1)
-      {
-        lowercoords[0] = lower;
-        uppercoords[0] = upper;
-      }
-      break;
-    case 2:
-      name = "cart2d";
-      ierr = CCTK_CoordRange(GH, &lower, &upper, 1, NULL, name);
-      if (ierr>-1)
-      {
-        lowercoords[0] = lower;
-        uppercoords[0] = upper;
-      }
-      ierr = CCTK_CoordRange(GH, &lower, &upper, 2, NULL, name);
-      if (ierr>-1)
-      {
-        lowercoords[1] = lower;
-        uppercoords[1] = upper;
-      }
-      break;
-    case 3:
-      name = "cart3d";
-      ierr = CCTK_CoordRange(GH, &lower, &upper, 1, NULL, name);
-      if (ierr>-1)
-      {
-        lowercoords[0] = lower;
-        uppercoords[0] = upper;
-      }
-      ierr = CCTK_CoordRange(GH, &lower, &upper, 2, NULL, name);
-      if (ierr>-1)
-      {
-        lowercoords[1] = lower;
-        uppercoords[1] = upper;
-      }
-      ierr = CCTK_CoordRange(GH, &lower, &upper, 3, NULL, name);
-      if (ierr>-1)
-      {
-        lowercoords[2] = lower;
-        uppercoords[2] = upper;
-      }
-      break;
-    default:
-      CCTK_VWarn(4,__LINE__,__FILE__,"IOASCII",
-                 "IOUtil_1DLines: Support only for dim<=3");
-      ierr = -1;
-      break;
+    CCTK_VWarn (4, __LINE__, __FILE__, CCTK_THORNSTRING,
+                "IOUtil_1DLines: Cartesian coordinate system '%s' not found",
+                system_name);
+    memset (slice_center, 0, num_dims * num_dims * sizeof (int));
+    return (-1);
+  }
+
+  /* allocate arrays for ranges */
+  lower_range = (CCTK_REAL *) 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, system_name) < 0)
+    {
+      CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+                  "IOUtil_1DLines: Could not get ranges for %c-direction "
+                  "of coordinate system '%s'",
+                  'x' + dir, system_name);
+    }
   }
 
-  for (dir = 0; dir < dimension; dir++)
+  /* now set the slice center for each line
+     according to origin_index[] or origin_phys[] */
+  for (dir = 0; dir < num_dims; dir++)
   {
-    for (dim = 0; dim < dimension; dim++)
+    for (dim = 0; dim < num_dims; dim++)
     {
       if (dim == dir)
       {
-        spxyz[dir][dim] = 0;
+        slice_center[dir][dim] = 0;
       }
-      else if (index && index[dir][dim]>-1)
+      else if (origin_index && origin_index[dir][dim] >= 0)
       {
-        spxyz[dir][dim] = index[dir][dim];
+        slice_center[dir][dim] = origin_index[dir][dim];
       }
-      else if (ierr < 0)
+      else if (lower_range[dim] > origin_phys[dir][dim] || 
+               upper_range[dim] < origin_phys[dir][dim])
       {
-        CCTK_VWarn(4, __LINE__,__FILE__,"IOASCII",
-                   "IOUtil_1DLines: Cartesian coordinate system %s not found",
-                   name);
-        spxyz[dir][dim] = 0;
-      }
-      else if (lowercoords[dim]>coord[dir][dim] || 
-               uppercoords[dim]<coord[dir][dim])
-      {
-        CCTK_VWarn(2,__LINE__,__FILE__,"IOUtil",
-                   "IOUtil_1DLines: Coordinate in direction %d (%f,%f)"
-                   " doesn't contain %f",
-                   dim,lowercoords[dim],uppercoords[dim],coord[dir][dim]);
-        spxyz[dir][dim] = 0;
+        CCTK_VWarn (2, __LINE__, __FILE__, CCTK_THORNSTRING,
+                    "IOUtil_1DLines: Coordinate in %c-direction (%f, %f) "
+                    "doesn't contain %f",
+                    'x' + dim, (double) lower_range[dim],
+                    (double) upper_range[dim], (double) origin_phys[dir][dim]);
+        slice_center[dir][dim] = 0;
       }
       else 
       {
         /* Find index for first point above the chosen coordinate */
-        spxyz[dir][dim] = 0.5-(lowercoords[dim]-coord[dir][dim])/
-          (GH->cctk_delta_space[dim]);
-        
+        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 %d-coord of lines in %d-direction is %d\n",
-               dim,dir,spxyz[dir][dim]);
+               dim,dir,slice_center[dir][dim]);
 #endif
       }
     }
   }
 
-  return 0;
+  /* free allocated resources */
+  free (lower_range);
+
+  return (0);
 }
 
 
+/*@@
+   @routine    IOUtil_2DPlanes
+   @date       July 4 2000
+   @author     Gabrielle Allen, Gerd Lanfermann, Thomas Radke
+   @desc
+               Fills out an array determining where to position 2D planes
+               for output on a multidimensional regular Cartesian unigrid.
+   @enddesc
 
-int IOUtil_2DPlanes (cGH *GH, 
-                     int dimension, 
-                     int *index, 
-                     CCTK_REAL *coord, 
-                     int *sp2xyz)
+   @calls      CCTK_CoordRange
+
+   @var        GH
+   @vdesc      Pointer to CCTK GH
+   @vtype      cGH *
+   @vio        in
+   @endvar
+   @var        num_dims
+   @vdesc      number of dimensions of underlying grid
+   @vtype      int
+   @vio        in
+   @endvar
+   @var        origin_index
+   @vdesc      origin for 2D planes given by grid point indices
+   @vtype      const int [num_dims]
+   @vio        in
+   @endvar
+   @var        origin_phys
+   @vdesc      origin for 2D planes given by physical coordinates
+   @vtype      const CCTK_REAL [num_dims]
+   @vio        in
+   @endvar
+   @var        slice_center
+   @vdesc      resulting 2D plane slice center
+   @vtype      int [num_dims]
+   @vio        out
+   @endvar
+
+   @returntype int
+   @returndesc
+               0  for success
+               -1 no coordinate system of given dimensions found
+   @endreturndesc
+ @@*/
+int IOUtil_2DPlanes (cGH *GH,
+                     int num_dims,
+                     const int origin_index[],
+                     const CCTK_REAL origin_phys[],
+                     int slice_center[])
 {
   DECLARE_CCTK_PARAMETERS
-  int ierr=0;
   int dir;
-  CCTK_REAL lower,upper;
-  const char *name=NULL;
-  CCTK_REAL lowercoords[3]={0,0,0};
-  CCTK_REAL uppercoords[3]={0,0,0};
+  char system_name[20];
+  CCTK_REAL *lower_range, *upper_range;
 
 
-  /* Get the appropriate coordinate system */
-  /* Get the appropriate coordinate system */
-  switch (dimension)
+  /* get the appropriate coordinate system name */
+  sprintf (system_name, "cart%dd", num_dims);
+
+  if (CCTK_CoordSystemHandle (system_name) < 0)
   {
-    case 1:
-      name = "cart1d";
-      ierr = CCTK_CoordRange(GH, &lower, &upper, 1, NULL, name);
-      if (ierr>-1)
-      {
-        lowercoords[0] = lower;
-        uppercoords[0] = upper;
-      }
-      break;
-    case 2:
-      name = "cart2d";
-      ierr = CCTK_CoordRange(GH, &lower, &upper, 1, NULL, name);
-      if (ierr>-1)
-      {
-        lowercoords[0] = lower;
-        uppercoords[0] = upper;
-      }
-      ierr = CCTK_CoordRange(GH, &lower, &upper, 2, NULL, name);
-      if (ierr>-1)
-      {
-        lowercoords[1] = lower;
-        uppercoords[1] = upper;
-      }
-      break;
-    case 3:
-      name = "cart3d";
-      ierr = CCTK_CoordRange(GH, &lower, &upper, 1, NULL, name);
-      if (ierr>-1)
-      {
-        lowercoords[0] = lower;
-        uppercoords[0] = upper;
-      }
-      ierr = CCTK_CoordRange(GH, &lower, &upper, 2, NULL, name);
-      if (ierr>-1)
-      {
-        lowercoords[1] = lower;
-        uppercoords[1] = upper;
-      }
-      ierr = CCTK_CoordRange(GH, &lower, &upper, 3, NULL, name);
-      if (ierr>-1)
-      {
-        lowercoords[2] = lower;
-        uppercoords[2] = upper;
-      }
-      break;
-    default:
-      CCTK_VWarn(4,__LINE__,__FILE__,"IOASCII",
-                 "IOUtil_2DPlanes: Support only for dim<=3");
-      ierr = -1;
-      break;
+    CCTK_VWarn (4, __LINE__, __FILE__, CCTK_THORNSTRING,
+                "IOUtil_1DLines: Cartesian coordinate system '%s' not found",
+                system_name);
+    memset (slice_center, 0, num_dims * sizeof (int));
+    return (-1);
   }
 
+  /* allocate arrays for ranges */
+  lower_range = (CCTK_REAL *) calloc (2 * num_dims, sizeof (CCTK_REAL));
+  upper_range = lower_range + num_dims;
 
-  for (dir = 0; dir < dimension; dir++)
+  /* get the ranges in every direction */
+  for (dir = 0; dir < num_dims; dir++)
   {
-    if (index && index[dir]>-1)
+    if (CCTK_CoordRange (GH, &lower_range[dir], &upper_range[dir],
+                         dir + 1, NULL, system_name) < 0)
     {
-      sp2xyz[dir] = index[dir];
+      CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+                  "IOUtil_2DPlanes: Could not get ranges for %c-direction "
+                  "of coordinate system '%s'",
+                  'x' + dir, system_name);
     }
-    else if (ierr < 0) 
+  }
+
+  /* now set the slice center for each line
+     according to origin_index[] or origin_phys[] */
+  for (dir = 0; dir < num_dims; dir++)
+  {
+    if (origin_index && origin_index[dir] >= 0)
     {
-      CCTK_VWarn(4, __LINE__,__FILE__,"IOUtil",
-                 "IOUtil_2DPlanes: Cartesian coordinate system %s not found",
-                 name);
-      sp2xyz[dir] = 0;
+      slice_center[dir] = origin_index[dir];
     }
-    else if (lowercoords[dir]>coord[dir] || 
-             uppercoords[dir]<coord[dir])
+    else if (lower_range[dir] > origin_phys[dir] || 
+             upper_range[dir] < origin_phys[dir])
     {
-      CCTK_VWarn(2,__LINE__,__FILE__,"IOUtil",
-                 "IOUtil_2DPlanes: Coordinate in direction %d (%lf,%lf)"
-                 " doesn't contain %f",
-                 dir,lowercoords[dir],uppercoords[dir],coord[dir]);
-      sp2xyz[dir] = 0;
+      CCTK_VWarn (2, __LINE__, __FILE__, CCTK_THORNSTRING,
+                  "IOUtil_2DPlanes: Coordinate in %c-direction (%f, %f) "
+                  "doesn't contain %f",
+                  'x' + dir, (double) lower_range[dir],
+                  (double) upper_range[dir], (double) origin_phys[dir]);
+      slice_center[dir] = 0;
     }
     else
     {
       /* Find index for first point above the chosen coordinate */
-      sp2xyz[dir] = 0.5-(lowercoords[dir]-coord[dir])/
-        (GH->cctk_delta_space[dir]);
-      
+      slice_center[dir] = ceil ((origin_phys[dir] - lower_range[dir]) /
+                                GH->cctk_delta_space[dir] - 1e-6);
+
 #ifdef DEBUG_IOUTIL
       printf("sp2xyz for planes perpendicular to"
              " %d-direction is %d\n",
-             dir,sp2xyz[dir]);
+             dir,slice_center[dir]);
 #endif
     }
   }
 
-  return 0;
+  /* free allocated resources */
+  free (lower_range);
 
+  return (0);
 }