Fixed wrong defaults for 2D slice center setup.

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

1 files changed, 21 insertions, 20 deletions

diff --git a/src/Utils.c b/src/Utils.c
index 7c06e83..4d1cd1b 100644
--- a/src/Utils.c
+++ b/src/Utils.c
@@ -303,7 +303,7 @@ int IOUtil_1DLines (const cGH *GH,
   int dim, dir, coord_system_handle, have_coords, have_origin_index, len;
   char *coord_system_name;
   CCTK_INT *coord_handles;
-  CCTK_REAL *lower, *upper;
+  CCTK_REAL *lower, *upper, *delta;
 
 
   /* get the default coordinate system associated with this grid dimension */
@@ -356,14 +356,16 @@ int IOUtil_1DLines (const cGH *GH,
   }
 
   /* get the ranges in every direction */
-  lower = calloc (2 * num_dims, sizeof (CCTK_REAL));
-  upper = lower + num_dims;
+  lower = calloc (3 * num_dims, sizeof (CCTK_REAL));
+  upper = lower + 1*num_dims;
+  delta = lower + 2*num_dims;
   if (! have_origin_index)
   {
     for (dir = 0; dir < num_dims; dir++)
     {
       if (Util_TableGetReal (coord_handles[dir], &lower[dir], "COMPMIN") < 0 ||
-          Util_TableGetReal (coord_handles[dir], &upper[dir], "COMPMAX") < 0)
+          Util_TableGetReal (coord_handles[dir], &upper[dir], "COMPMAX") < 0 ||
+          Util_TableGetReal (coord_handles[dir], &delta[dir], "DELTA") < 0)
       {
         CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING,
                     "IOUtil_1DLines: Could not get ranges for %c-direction "
@@ -408,8 +410,7 @@ int IOUtil_1DLines (const cGH *GH,
       {
         /* Find index for first point above the chosen coordinate */
         slice_center[dir][dim] =
-          ceil ((origin_phys[dir][dim] - lower[dim]) /
-                GH->cctk_delta_space[dim] - 1e-6);
+          ceil ((origin_phys[dir][dim] - lower[dim]) / delta[dim] - 1e-6);
 
 #ifdef DEBUG_IOUTIL
         printf("spxyz for %c-coord of lines in %c-direction is %d\n",
@@ -479,7 +480,7 @@ int IOUtil_2DPlanes (const cGH *GH,
   int have_coords, have_origin_index, len, coord_system_handle, dir;
   char *coord_system_name;
   CCTK_INT *coord_handles;
-  CCTK_REAL *lower, *upper;
+  CCTK_REAL *lower, *upper, *delta;
 
 
   /* get the default coordinate system associated with this grid dimension */
@@ -525,18 +526,20 @@ int IOUtil_2DPlanes (const cGH *GH,
   }
 
   /* get the ranges in every direction */
-  lower = calloc (2 * num_dims, sizeof (CCTK_REAL));
-  upper = lower + num_dims;
+  lower = calloc (3 * num_dims, sizeof (CCTK_REAL));
+  upper = lower + 1*num_dims;
+  delta = lower + 2*num_dims;
   if (have_coords)
   {
     for (dir = 0; dir < num_dims; dir++)
     {
       if (Util_TableGetReal (coord_handles[dir], &lower[dir], "COMPMIN") < 0 ||
-          Util_TableGetReal (coord_handles[dir], &upper[dir], "COMPMAX") < 0)
+          Util_TableGetReal (coord_handles[dir], &upper[dir], "COMPMAX") < 0 ||
+          Util_TableGetReal (coord_handles[dir], &delta[dir], "DELTA") < 0)
       {
         CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING,
-                    "IOUtil_2DPlanes: Could not get ranges for %c-direction "
-                    "of associated coordinate system '%s'",
+                    "IOUtil_2DPlanes: Could not get ranges/deltas for "
+                    "%c-direction of associated coordinate system '%s'",
                     'x' + (num_dims - dir - 1), coord_system_name);
       }
     }
@@ -550,14 +553,14 @@ int IOUtil_2DPlanes (const cGH *GH,
     {
       slice_center[dir] = origin_index[dir];
     }
-    else if (lower[dir] > origin_phys[dir] ||
-             upper[dir] < origin_phys[dir])
+    else if (lower[num_dims-1-dir] > origin_phys[dir] ||
+             upper[num_dims-1-dir] < origin_phys[dir])
     {
       CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
                   "IOUtil_2DPlanes: %c-coordinate for slice center of 2D "
                   "planes (%f) is out of grid coordinates range (%f, %f)",
                   'x' + (num_dims - dir - 1), (double) origin_phys[dir],
-                  (double) lower[dir], (double) upper[dir]);
+                  (double) lower[num_dims-1-dir], (double) upper[num_dims-1-dir]);
       CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
                   "IOUtil_2DPlanes: slice center will default to %c-index 0",
                   'x' + (num_dims - dir - 1));
@@ -566,13 +569,11 @@ int IOUtil_2DPlanes (const cGH *GH,
     else
     {
       /* Find index for first point above the chosen coordinate */
-      slice_center[dir] = ceil ((origin_phys[dir] - lower[dir]) /
-                                GH->cctk_delta_space[(num_dims - dir - 1)] -
-                                1e-6);
-
+      slice_center[dir] = ceil ((origin_phys[dir] - lower[num_dims-1-dir]) /
+                                delta[num_dims-1-dir] - 1e-6);
 #ifdef DEBUG_IOUTIL
       printf("sp2xyz for planes perpendicular to %d-direction is %d\n",
-             (num_dims - dir - 1), (CCTK_INT) slice_center[dir]);
+             dir, slice_center[dir]);
 #endif
     }
   }