Don't get hyperslabs from coordinate grid variables.

Get the bounding box information directly from coordinate tables.
This should fix wrong bboxes in the output.


git-svn-id: http://svn.cactuscode.org/arrangements/CactusIO/IOSDF/trunk@5 3af55ef0-e5e4-43b4-b020-ca5761ff09b2

4 files changed, 113 insertions, 200 deletions

diff --git a/src/Write1D.c b/src/Write1D.c
index debf166..31073b4 100644
--- a/src/Write1D.c
+++ b/src/Write1D.c
@@ -76,22 +76,19 @@ static char *OpenFile (const cGH *GH,
 int IOSDF_Write1D (const cGH *GH, int vindex, const char *alias)
 {
   ioSDFGH *myGH;
-  int do_dir[4], coord_index[3];
+  int do_dir[4];
   int i, dir, myproc, gindex, have_coords, mapping;
-  int num_requested_hslabs, num_returned_hslabs;
   int *extent_int;
   cGroup gdata;
   char *fullname, *groupname, *filename;
   CCTK_INT coord_system_handle, coord_handles[3];
-  double offset;
-  CCTK_REAL coord_lower[3];
+  CCTK_REAL minext[3], delta[3];
   CCTK_INT downsample[4];
   CCTK_INT *origin, *direction;
   CCTK_INT hsize, extent;
-  CCTK_INT vindices[2];
-  double *hdata[2];
+  double *hdata;
+  double bbox[2];
   const double dtime = GH->cctk_time;
-  const char *coordnames[] = {"x", "y", "z", "d"};
   DECLARE_CCTK_PARAMETERS
 
 
@@ -141,24 +138,6 @@ int IOSDF_Write1D (const cGH *GH, int vindex, const char *alias)
   coord_system_handle = Coord_GroupSystem (GH, groupname);
   free (groupname);
 
-  dir = gdata.dim < 3 ? gdata.dim : 3;
-
-  have_coords = coord_system_handle >= 0 &&
-                Util_TableGetIntArray (coord_system_handle, dir,
-                                       coord_handles, "COORDINATES") >= 0;
-  if (have_coords)
-  {
-    /* get the coordinate functions and coordinate physical minimum */
-    for (i = 0; i < dir; i++)
-    {
-      coord_index[i] = -1;
-      coord_lower[i] = 0;
-      Util_TableGetInt (coord_handles[i], &coord_index[i], "GAINDEX");
-      Util_TableGetReal (coord_handles[i], &coord_lower[i], "COMPMIN");
-      have_coords &= coord_index[i] >= 0;
-    }
-  }
-
   myproc = CCTK_MyProc (GH);
 
   origin     = calloc (2*gdata.dim, sizeof (CCTK_INT));
@@ -171,6 +150,22 @@ int IOSDF_Write1D (const cGH *GH, int vindex, const char *alias)
   downsample[2] = out_downsample_z;
   downsample[3] = 1;
 
+  dir = gdata.dim < 3 ? gdata.dim : 3;
+
+  have_coords = coord_system_handle >= 0 &&
+                Util_TableGetIntArray (coord_system_handle, dir,
+                                       coord_handles, "COORDINATES") >= 0;
+  if (have_coords)
+  {
+    /* get the coordinates ranges */
+    for (i = 0; i < dir; i++)
+    {
+      Util_TableGetReal (coord_handles[i], &minext[i], "COMPMIN");
+      Util_TableGetReal (coord_handles[i], &delta[i], "DELTA");
+      delta[i] *= downsample[i];
+    }
+  }
+
   /* get the variable's extents, compute the extent for 3D-diagonals as the
      minimum of grid points in each direction */
   CCTK_GroupgshVI (GH, gdata.dim, extent_int, vindex);
@@ -192,13 +187,8 @@ int IOSDF_Write1D (const cGH *GH, int vindex, const char *alias)
   /* now do the actual I/O looping over all directions */
   for (dir = 0; dir < 4; dir++)
   {
-    if (hsize <= 0)
-    {
-      continue;
-    }
-
     /* skip empty slices */
-    if (! do_dir[dir])
+    if (hsize <= 0 || ! do_dir[dir])
     {
       continue;
     }
@@ -210,10 +200,6 @@ int IOSDF_Write1D (const cGH *GH, int vindex, const char *alias)
       filename = OpenFile (GH, fullname, alias, &gdata, dir);
     }
 
-    /* get the number of hyperslabs to extract
-       (ie. whether to get a coordinate hyperslab too or not) */
-    num_requested_hslabs = have_coords && dir < 3 ? 2 : 1;
-
     /* set the direction vector */
     for (i = 0; i < gdata.dim; i++)
     {
@@ -233,9 +219,9 @@ int IOSDF_Write1D (const cGH *GH, int vindex, const char *alias)
       extent -= origin[dir];
 
       /* correct extent in the case of staggered grids */
-      if (CCTK_StaggerDirIndex(dir,gdata.stagtype)==1)
+      if (CCTK_StaggerDirIndex (dir, gdata.stagtype) == 1)
       {
-        --extent;
+        extent--;
       }
     }
     else    /* origin for CCTK_ARRAYS is always (0, 0, 0) */
@@ -265,26 +251,17 @@ int IOSDF_Write1D (const cGH *GH, int vindex, const char *alias)
       continue;
     }
 
-    /* allocate hyperslab buffers on I/O processor */
-    hdata[0] = hdata[1] = NULL;
-    if (myproc == 0)
-    {
-      hdata[0] = malloc (hsize * sizeof (double));
-      hdata[1] = have_coords ? malloc (hsize * sizeof (double)) : NULL;
-    }
+    /* allocate hyperslab buffer on I/O processor */
+    hdata = myproc == 0 ? malloc (hsize * sizeof (double)) : NULL;
 
-    /* get the hyperslabs */
-    vindices[0] = vindex;
-    vindices[1] = coord_index[dir];
-    num_returned_hslabs = Hyperslab_GetList (GH, mapping, num_requested_hslabs,
-                                             NULL, vindices, NULL, NULL,
-                                             (void **) hdata, NULL);
+    /* get the hyperslab */
+    i = Hyperslab_Get (GH, mapping, 0, vindex, 0, gdata.vartype, hdata);
 
     /* release the mapping structure */
     Hyperslab_FreeMapping (mapping);
 
     /* And dump the data to file */
-    if (num_returned_hslabs != num_requested_hslabs)
+    if (i != 0)
     {
       CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
                   "IOSDF_Write1D: Failed to extract hyperslab for "
@@ -292,35 +269,17 @@ int IOSDF_Write1D (const cGH *GH, int vindex, const char *alias)
     }
     else if (filename)
     {
+      bbox[0] = -1; bbox[1] = +1;
       if (have_coords)
       {
-        if (dir < 3)
-        {
-          /* get the staggering offset for the xyz coordinates */
-          offset = 0.5 * GH->cctk_delta_space[dir] *
-                   CCTK_StaggerDirIndex (dir, gdata.stagtype);
-          for (i = 0; i < hsize; i++)
-          {
-            hdata[1][i] += offset;
-          }
-        }
-        else
-        {
-          /* calculate the diagonal coordinates */
-          offset = GH->cctk_delta_space[0] * sqrt (3);
-          for (i = 0; i < hsize; i++)
-          {
-            hdata[1][i] = coord_lower[0]*sqrt (3) + i*offset;
-          }
-        }
+        bbox[0] = minext[dir];
+        bbox[1] = dir < 3 ? bbox[0] + delta[dir]*(hsize-1) :
+                  (minext[0] + hsize*GH->cctk_delta_space[0]) * sqrt (3);
       }
+      gft_out_set_bbox (bbox, 1);
 
       extent_int[0] = hsize;
-      i = have_coords ? gft_out_full (filename, dtime, extent_int,
-                                      coordnames[dir], 1, hdata[1], hdata[0]):
-                        gft_out_brief (filename, dtime, extent_int,
-                                       1, hdata[0]);
-      if (i <= 0)
+      if (gft_out_brief (filename, dtime, extent_int, 1, hdata) <= 0)
       {
         CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
                     "Error writing to 1D IOSDF output file '%s'", filename);
@@ -328,8 +287,8 @@ int IOSDF_Write1D (const cGH *GH, int vindex, const char *alias)
     }
 
     /* clean up */
-    free (hdata[0]);
-    free (hdata[1]);
+    free (hdata);
+
   } /* end of loop through all directions */
 
   /* free allocated resources */
diff --git a/src/Write2D.c b/src/Write2D.c
index 551e0c8..2906c79 100644
--- a/src/Write2D.c
+++ b/src/Write2D.c
@@ -43,7 +43,7 @@ static char **OpenFile (const cGH *GH,
    @enddesc
    @calls     Hyperslab_GlobalMappingByIndex
               Hyperslab_FreeMapping
-              Hyperslab_GetList
+              Hyperslab_Get
               OpenFile
               WriteData
 
@@ -75,21 +75,19 @@ static char **OpenFile (const cGH *GH,
 int IOSDF_Write2D (const cGH *GH, int vindex, const char *alias)
 {
   ioSDFGH *myGH;
-  int i, total_hsize, num_requested_hslabs, num_returned_hslabs;
+  int i, total_hsize;
   int dir, dir_i, dir_j, maxdir, myproc, gindex, have_coords;
   int mapping;
   cGroup gdata;
-  int coord_index[3];
   CCTK_INT coord_system_handle, coord_handles[3];
   char *fullname, *groupname;
-  int extent_int[3];
-  double offset[2];
-  CCTK_INT vindices[3], origin[3], extent[2], direction[6], downsample[2],
-           hsize[2];
-  double *hdata[3];
+  int full_extent[3], slice_extent[2];
+  double bbox[2*2];
+  CCTK_INT origin[3], extent[2], direction[6], downsample[2], hsize[2];
+  CCTK_REAL minext[2], delta[2];
+  double *hdata;
   char **filenames;
   const double dtime = GH->cctk_time;
-  const char *coordnames[] = {"x|y", "x|z", "y|z"};
   DECLARE_CCTK_PARAMETERS
 
 
@@ -124,17 +122,6 @@ int IOSDF_Write2D (const cGH *GH, int vindex, const char *alias)
   have_coords = coord_system_handle >= 0 &&
                 Util_TableGetIntArray (coord_system_handle, dir,
                                        coord_handles, "COORDINATES") >= 0;
-  if (have_coords)
-  {
-    /* get the coordinate functions and coordinate physical minimum */
-    for (i = 0; i < dir; i++)
-    {
-      coord_index[i] = -1;
-      Util_TableGetInt (coord_handles[i], &coord_index[i], "GAINDEX");
-      have_coords &= coord_index[i] >= 0;
-    }
-  }
-  num_requested_hslabs = have_coords ? 3 : 1;
 
   /* processor 0 opens the files on the first trip through */
   filenames = NULL;
@@ -145,12 +132,12 @@ int IOSDF_Write2D (const cGH *GH, int vindex, const char *alias)
   }
 
   /* get the extents of the variable */
-  CCTK_GroupgshVI (GH, gdata.dim, extent_int, vindex);
+  CCTK_GroupgshVI (GH, gdata.dim, full_extent, vindex);
 
   /* get the total number of grid points to check for zero-sized variables */
   for (dir = 0, hsize[0] = 1; dir < gdata.dim; dir++)
   {
-    hsize[0] *= extent_int[dir];
+    hsize[0] *= full_extent[dir];
   }
 
   /* now do the actual I/O looping over all directions */
@@ -179,8 +166,8 @@ int IOSDF_Write2D (const cGH *GH, int vindex, const char *alias)
     }
 
     /* set the extent vector */
-    extent[0] = extent_int[dir_i];
-    extent[1] = extent_int[dir_j];
+    extent[0] = full_extent[dir_i];
+    extent[1] = full_extent[dir_j];
 
     /* set the origin using the slice center from IOUtil */
     memset (origin, 0, sizeof (origin));
@@ -216,22 +203,30 @@ int IOSDF_Write2D (const cGH *GH, int vindex, const char *alias)
       continue;
     }
 
-    hdata[0] = hdata[1] = hdata[2];
-    if (myproc == 0)
+    /* get the bounding box information */
+    bbox[0] = bbox[2] = -1;
+    bbox[1] = bbox[3] = +1;
+    if (have_coords)
     {
-      /* allocate hyperslab buffers */
-      hdata[0] = malloc (total_hsize * sizeof (double));
-      hdata[1] = have_coords ? malloc (2 * total_hsize * sizeof(double)) : NULL;
-      hdata[2] = hdata[1] + total_hsize;
+      Util_TableGetReal (coord_handles[dir_i], &minext[0], "COMPMIN");
+      Util_TableGetReal (coord_handles[dir_i], &delta[0], "DELTA");
+      Util_TableGetReal (coord_handles[dir_j], &minext[1], "COMPMIN");
+      Util_TableGetReal (coord_handles[dir_j], &delta[1], "DELTA");
+
+      delta[0] *= downsample[0];
+      delta[1] *= downsample[1];
+      bbox[0] = minext[0];
+      bbox[1] = bbox[0] + delta[0]*(hsize[0]-1);
+      bbox[2] = minext[1];
+      bbox[3] = bbox[2] + delta[1]*(hsize[1]-1);
     }
+    gft_out_set_bbox (bbox, 2);
 
-    /* get the hyperslabs */
-    vindices[0] = vindex;
-    vindices[1] = coord_index[dir_i];
-    vindices[2] = coord_index[dir_j];
-    num_returned_hslabs = Hyperslab_GetList (GH, mapping, num_requested_hslabs,
-                                             NULL, vindices, NULL, NULL,
-                                             (void **) hdata, NULL);
+    /* allocate hyperslab buffer */
+    hdata = myproc == 0 ? malloc (total_hsize * sizeof (double)) : NULL;
+
+    /* get the hyperslab */
+    i = Hyperslab_Get (GH, mapping, 0, vindex, 0, gdata.vartype, hdata);
 
     /* release the mapping structure */
     Hyperslab_FreeMapping (mapping);
@@ -239,30 +234,12 @@ int IOSDF_Write2D (const cGH *GH, int vindex, const char *alias)
     /* and dump the data to file */
     if (filenames)
     {
-      if (num_returned_hslabs == num_requested_hslabs)
+      if (i == 0)
       {
-        if (have_coords)
-        {
-          /* get the staggering offset for the coordinates */
-          offset[0] = 0.5 * GH->cctk_delta_space[dir_i] *
-                      CCTK_StaggerDirIndex (dir_i, gdata.stagtype);
-          offset[1] = 0.5 * GH->cctk_delta_space[dir_j] *
-                      CCTK_StaggerDirIndex (dir_j, gdata.stagtype);
-          for (i = 0; i < total_hsize; i++)
-          {
-            hdata[1][i] += offset[0];
-            hdata[2][i] += offset[1];
-          }
-        }
-
-        extent[0] = hsize[0];
-        extent[1] = hsize[1];
+        slice_extent[0] = hsize[0];
+        slice_extent[1] = hsize[1];
 
-        i = have_coords ? gft_out_full (filenames[dir], dtime, extent_int,
-                                        coordnames[dir], 2, hdata[1], hdata[0]):
-                          gft_out_brief (filenames[dir], dtime, extent_int,
-                                         2, hdata[0]);
-        if (i <= 0)
+        if (gft_out_brief (filenames[dir], dtime, slice_extent, 2, hdata) <= 0)
         {
           CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
                       "Error writing to 2D IOSDF output file '%s'",
@@ -277,8 +254,7 @@ int IOSDF_Write2D (const cGH *GH, int vindex, const char *alias)
       }
 
       /* clean up */
-      free (hdata[0]);
-      free (hdata[1]);
+      free (hdata);
 
     } /* end of outputting the data by processor 0 */
 
diff --git a/src/Write3D.c b/src/Write3D.c
index f07ef55..314146b 100644
--- a/src/Write3D.c
+++ b/src/Write3D.c
@@ -38,7 +38,7 @@ static char *OpenFile (const cGH *GH, const char *fullname, const char *alias);
    @enddesc
    @calls     Hyperslab_GlobalMappingByIndex
               Hyperslab_FreeMapping
-              Hyperslab_GetList
+              Hyperslab_Get
               OpenFile
 
    @var       GH
@@ -69,16 +69,16 @@ static char *OpenFile (const cGH *GH, const char *fullname, const char *alias);
 int IOSDF_Write3D (const cGH *GH, int vindex, const char *alias)
 {
   int i, total_hsize;
-  int myproc, gindex, have_coords;
-  int num_requested_hslabs, num_returned_hslabs;
+  int myproc, gindex;
   cGroup gdata;
   CCTK_INT coord_system_handle, coord_handles[3];
   char *fullname, *groupname, *filename;
-  double *hdata[4];
+  double *hdata;
   int extent_int[3];
-  double offset[3];
   int mapping;
-  CCTK_INT vindices[4], extent[3], downsample[3], hsize[3];
+  CCTK_INT extent[3], downsample[3], hsize[3];
+  double bbox_dbl[2*3];
+  CCTK_REAL bbox[2*3], delta[3];
   const double dtime = GH->cctk_time;
   const CCTK_INT origin[] = {0, 0, 0},
                  direction[] = {1, 0, 0, 0, 1, 0, 0, 0, 1};
@@ -100,32 +100,10 @@ int IOSDF_Write3D (const cGH *GH, int vindex, const char *alias)
   }
 
   /* get the coordinate system associated with this grid variable */
-  vindices[0] = vindex;
   groupname = CCTK_GroupName (gindex);
   coord_system_handle = Coord_GroupSystem (GH, groupname);
   free (groupname);
 
-  have_coords = coord_system_handle >= 0 &&
-                Util_TableGetIntArray (coord_system_handle, 3,
-                                       coord_handles, "COORDINATES") >= 0;
-  if (have_coords)
-  {
-    /* get the coordinate functions and coordinate physical minimum */
-    for (i = 1; i <= 3; i++)
-    {
-      vindices[i] = -1;
-      Util_TableGetInt (coord_handles[i-1], &vindices[i], "GAINDEX");
-      have_coords &= vindices[i] >= 0;
-    }
-  }
-  num_requested_hslabs = have_coords ? 4 : 1;
-
-  /* What processor are we on? */
-  myproc = CCTK_MyProc (GH);
-
-  /* Open the file on processor 0 */
-  filename = myproc == 0 ? OpenFile (GH, fullname, alias) : NULL;
-
   /* get the total number of grid points to check for zero-sized variables */
   /* set the extent vector (copy from 'int' to 'CCTK_INT') */
   CCTK_GroupgshVI (GH, 3, extent_int, vindex);
@@ -144,6 +122,12 @@ int IOSDF_Write3D (const cGH *GH, int vindex, const char *alias)
   downsample[1] = out_downsample_y;
   downsample[2] = out_downsample_z;
 
+  /* What processor are we on? */
+  myproc = CCTK_MyProc (GH);
+
+  /* Open the file on processor 0 */
+  filename = myproc == 0 ? OpenFile (GH, fullname, alias) : NULL;
+
   /* get the hyperslab mapping */
   mapping = Hyperslab_GlobalMappingByIndex (GH, vindex, 3,
                                             direction, origin, extent,
@@ -174,46 +158,40 @@ int IOSDF_Write3D (const cGH *GH, int vindex, const char *alias)
     return (-1);
   }
 
-  hdata[0] = hdata[1] = hdata[2] = hdata[3];
-  if (myproc == 0)
+  /* get the bounding box information */
+  for (i = 0; i < 3; i++)
   {
-    /* allocate hyperslab buffers */
-    hdata[0] = malloc (total_hsize * sizeof (double));
-    hdata[1] = have_coords ? malloc (3 * total_hsize * sizeof (double)) : NULL;
-    hdata[2] = hdata[1] + 1*total_hsize;
-    hdata[3] = hdata[1] + 2*total_hsize;
+    bbox_dbl[2*i + 0] = -1;
+    bbox_dbl[2*i + 1] = +1;
   }
+  if (coord_system_handle >= 0 &&
+      Util_TableGetIntArray (coord_system_handle, 3, coord_handles,
+                             "COORDINATES") >= 0)
+  {
+    for (i = 0; i < 3; i++)
+    {
+      Util_TableGetReal (coord_handles[i], &bbox[2*i + 0], "COMPMIN");
+      Util_TableGetReal (coord_handles[i], &delta[i], "DELTA");
+
+      delta[i] *= downsample[i];
+      bbox_dbl[2*i + 0] = bbox[2*i + 0];
+      bbox_dbl[2*i + 1] = bbox_dbl[2*i + 0] + delta[i]*(hsize[i]-1);
+    }
+  }
+  gft_out_set_bbox (bbox_dbl, 3);
+
+  /* allocate hyperslab buffer */
+  hdata = myproc == 0 ? malloc (total_hsize * sizeof (double)) : NULL;
  
-  /* get the hyperslabs */
-  num_returned_hslabs = Hyperslab_GetList (GH, mapping, num_requested_hslabs,
-                                           NULL, vindices, NULL, NULL,
-                                           (void **) hdata, NULL);
+  /* get the hyperslab */
+  i = Hyperslab_Get (GH, mapping, 0, vindex, 0, gdata.vartype, hdata);
 
   /* And dump the data to file */
   if (myproc == 0)
   {
-    if (num_returned_hslabs == num_requested_hslabs)
+    if (i == 0)
     {
-      if (have_coords)
-      {
-        /* get the staggering offset for the coordinates */
-        for (i = 0; i < 3; i++)
-        {
-          offset[i] = 0.5 * GH->cctk_delta_space[i] *
-                      CCTK_StaggerDirIndex (i, gdata.stagtype);
-        }
-        for (i = 0; i < total_hsize; i++)
-        {
-          hdata[1][i] += offset[0];
-          hdata[2][i] += offset[1];
-          hdata[3][i] += offset[2];
-        }
-      }
-
-      i = have_coords ? gft_out_full (filename, dtime, extent_int, "x|y|z", 3,
-                                      hdata[1], hdata[0]) :
-                        gft_out_brief (filename, dtime, extent_int, 3,hdata[0]);
-      if (i <= 0)
+      if (gft_out_brief (filename, dtime, extent_int, 3, hdata) <= 0)
       {
         CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
                     "Error writing to 3D IOSDF output file '%s'", filename);
@@ -227,8 +205,7 @@ int IOSDF_Write3D (const cGH *GH, int vindex, const char *alias)
     }
 
     /* clean up */
-    free (hdata[0]);
-    free (hdata[1]);
+    free (hdata);
   } /* end of outputting the data by processor 0 */
 
   free (fullname);
diff --git a/src/ioSDFGH.h b/src/ioSDFGH.h
index e31867b..3706d07 100644
--- a/src/ioSDFGH.h
+++ b/src/ioSDFGH.h
@@ -68,6 +68,7 @@ int IOSDF_Output3DVarAs (const cGH *GH, const char *var, const char *alias);
 
 /* other function prototypes */
 void IOSDF_Startup (void);
+void IOSDF_Terminate (void);
 
 int IOSDF_Write1D (const cGH *GH, int vindex, const char *alias);
 int IOSDF_Write2D (const cGH *GH, int vindex, const char *alias);