CarpetIOHDF5: hdf5toascii_slicer now allows 1D line extraction from 3D HDF5 files

darcs-hash:20060828173429-776a0-b231124b73645983b8ed56efa35d5ec2c9354add.gz

1 files changed, 53 insertions, 35 deletions

diff --git a/Carpet/CarpetIOHDF5/src/util/hdf5toascii_slicer.cc b/Carpet/CarpetIOHDF5/src/util/hdf5toascii_slicer.cc
index 7efb1207e..4d2210f74 100644
--- a/Carpet/CarpetIOHDF5/src/util/hdf5toascii_slicer.cc
+++ b/Carpet/CarpetIOHDF5/src/util/hdf5toascii_slicer.cc
@@ -3,9 +3,9 @@
    @date      Sun 30 July 2007
    @author    Thomas Radke
    @desc
-              Utility program to extract a 2D slice from 3D datasets
+              Utility program to extract a 1D line or 2D slice from 3D datasets
               in datafiles generated by CarpetIOHDF5.
-              The extracted slice is written to stdout in CarpetIOASCII format.
+              The extracted slab is written to stdout in CarpetIOASCII format.
    @enddesc
  @@*/
 
@@ -71,8 +71,8 @@ typedef struct {
  *************************       Global Data         *************************
  *****************************************************************************/
 
-// the slice coordinate as selected by the user
-static double slice_coord[3] = {PARAMETER_UNSET, PARAMETER_UNSET, PARAMETER_UNSET};
+// the slab coordinate as selected by the user
+static double slab_coord[3] = {PARAMETER_UNSET, PARAMETER_UNSET, PARAMETER_UNSET};
 
 // the specific timestep selected by the user
 static double timestep = PARAMETER_UNSET;
@@ -120,7 +120,7 @@ int main (int argc, char *const argv[])
 #if 0
   int iorigin[3] = {-1, -1, -1};
 #endif
-  int num_slice_options = 0, num_slicei_options = 0;
+  int num_slab_options = 0;
 
 
   // evaluate command line parameters
@@ -131,21 +131,30 @@ int main (int argc, char *const argv[])
       verbose = true;
     } else if (strcmp (argv[i], "--timestep") == 0 and i+1 < argc) {
       timestep = atof (argv[++i]);
-    } else if (strcmp (argv[i], "--out_precision") == 0 and i+1 < argc) {
+    } else if (strcmp (argv[i], "--out-precision") == 0 and i+1 < argc) {
       cout << setprecision (atoi (argv[++i]));
+    } else if (strcmp (argv[i], "--out1d-xline-yz") == 0 and i+2 < argc) {
+      slab_coord[1] = atof (argv[++i]);
+      slab_coord[2] = atof (argv[++i]); num_slab_options++;
+    } else if (strcmp (argv[i], "--out1d-yline-xz") == 0 and i+2 < argc) {
+      slab_coord[0] = atof (argv[++i]);
+      slab_coord[2] = atof (argv[++i]); num_slab_options++;
+    } else if (strcmp (argv[i], "--out1d-zline-xy") == 0 and i+2 < argc) {
+      slab_coord[0] = atof (argv[++i]);
+      slab_coord[1] = atof (argv[++i]); num_slab_options++;
     } else if (strcmp (argv[i], "--out2d-yzplane-x") == 0 and i+1 < argc) {
-      slice_coord[0] = atof (argv[++i]); num_slice_options++;
+      slab_coord[0] = atof (argv[++i]); num_slab_options++;
     } else if (strcmp (argv[i], "--out2d-xzplane-y") == 0 and i+1 < argc) {
-      slice_coord[1] = atof (argv[++i]); num_slice_options++;
+      slab_coord[1] = atof (argv[++i]); num_slab_options++;
     } else if (strcmp (argv[i], "--out2d-xyplane-z") == 0 and i+1 < argc) {
-      slice_coord[2] = atof (argv[++i]); num_slice_options++;
+      slab_coord[2] = atof (argv[++i]); num_slab_options++;
 #if 0
     } else if (strcmp (argv[i], "--out2d-yzplane-xi") == 0 and i+1 < argc) {
-      iorigin[0] = atoi (argv[++i]); num_slicei_options++;
+      iorigin[0] = atoi (argv[++i]); num_slab_options++;
     } else if (strcmp (argv[i], "--out2d-xzplane-yi") == 0 and i+1 < argc) {
-      iorigin[1] = atoi (argv[++i]); num_slicei_options++;
+      iorigin[1] = atoi (argv[++i]); num_slab_options++;
     } else if (strcmp (argv[i], "--out2d-xyplane-zi") == 0 and i+1 < argc) {
-      iorigin[2] = atoi (argv[++i]); num_slicei_options++;
+      iorigin[2] = atoi (argv[++i]); num_slab_options++;
 #endif
     } else {
       break;
@@ -153,7 +162,7 @@ int main (int argc, char *const argv[])
   }
 
   /* give some help if called with incorrect number of parameters */
-  if (help or i >= argc or num_slice_options + num_slicei_options != 1) {
+  if (help or i >= argc or num_slab_options != 1) {
     const string indent (strlen (argv[0]) + 1, ' ');
     cerr << endl << "   ---------------------------"
          << endl << "   Carpet HDF5 to ASCII Slicer"
@@ -161,20 +170,24 @@ int main (int argc, char *const argv[])
          << endl
          << "Usage: " << endl
          << argv[0] << " [--help]" << endl
-         << indent << "[--out_precision <digits>]" << endl
+         << indent << "[--out-precision <digits>]" << endl
          << indent << "[--timestep <cctk_time value>]" << endl
          << indent << "[--verbose]" << endl
-         << indent << "<--out2d-plane value>" << endl
+         << indent << "<--out1d-line value value> | <--out2d-plane value>" << endl
          << indent << "<hdf5_infiles>" << endl << endl
          << "  where" << endl
          << "    [--help]                         prints this help" << endl
-         << "    [--out_precision <digits>]       sets the output precision" << endl
+         << "    [--out-precision <digits>]       sets the output precision" << endl
          << "                                     for floating-point numbers" << endl
          << "    [--timestep <cctk_time value>]   selects all HDF5 datasets which" << endl
          << "                                     (fuzzily) match the specified time" << endl
          << "    [--verbose]                      lists skipped HDF5 datasets on stderr" << endl
-         << "  and <--out2d-plane value> must be specified" << endl
-         << "  as exactly one of the following options:" << endl
+         << "  and either <--out1d-line value value> or <--out2d-plane value> triples" << endl
+         << "  must be specified as in the following:" << endl
+         << "    --out1d-xline-yz  <origin_y> <origin_z>" << endl
+         << "    --out1d-yline-xz  <origin_x> <origin_z>" << endl
+         << "    --out1d-zline-xy  <origin_x> <origin_y>" << endl
+         << endl
          << "    --out2d-yzplane-x  <origin_x>" << endl
          << "    --out2d-xzplane-y  <origin_y>" << endl
          << "    --out2d-xyplane-z  <origin_z>" << endl
@@ -242,7 +255,7 @@ int main (int argc, char *const argv[])
    @desc
                The worker routine which is called by H5Giterate().
                It checks whether the current HDF5 object is a patch matching
-               the user's slice criteria, and adds it to the patches list.
+               the user's slab criteria, and adds it to the patches list.
    @enddesc
 
    @var        group
@@ -306,7 +319,7 @@ static herr_t FindPatches (hid_t group, const char *name, void *_file)
            << "  dataset has " << ndims << " dimensions" << endl;
     }
   } else if (timestep != PARAMETER_UNSET and
-             fabs (timestep - patch.time) / FUZZY_FACTOR > 1) {
+             fabs (timestep - patch.time) > FUZZY_FACTOR) {
     if (verbose) {
       cerr << "skipping dataset '" << name << "':" << endl
            << "  timestep (" << patch.time << ") doesn't match" << endl;
@@ -322,17 +335,17 @@ static herr_t FindPatches (hid_t group, const char *name, void *_file)
 
     int i;
     for (i = 0; i < 3; i++) {
-      if (slice_coord[i] != PARAMETER_UNSET) {
-        is_okay = slice_coord[i] >= patch.origin[i] and
-                  slice_coord[i] <= patch.origin[i] +
+      if (slab_coord[i] != PARAMETER_UNSET) {
+        is_okay = slab_coord[i] >= patch.origin[i] and
+                  slab_coord[i] <= patch.origin[i] +
                                     (patch.dims[2-i]-1)*patch.delta[i];
-        break;
+        if (not is_okay) break;
       }
     }
     if (not is_okay) {
       if (verbose) {
         cerr << "skipping dataset '" << name << "':" << endl
-             << "  slice " << slice_coord[i] << " is out of dataset range ["
+             << "  slab " << slab_coord[i] << " is out of dataset range ["
              << patch.origin[i] << ", "
              << patch.origin[i] + (patch.dims[2-i]-1)*patch.delta[i] << "]"
              << endl;
@@ -388,7 +401,7 @@ static herr_t FindPatches (hid_t group, const char *name, void *_file)
    @date       Tue 8 August 2006
    @author     Thomas Radke
    @desc
-               Reads a 2D slice from the given patch
+               Reads a slab from the given patch
                and outputs it in CarpetIOASCII format.
    @enddesc
 
@@ -419,8 +432,8 @@ static void ReadPatch (const patch_t& patch, int last_iteration)
   h5size_t slabstart[3] = {0, 0, 0};
   hsize_t slabcount[3] = {patch.dims[0], patch.dims[1], patch.dims[2]};
   for (int i = 0; i < 3; i++) {
-    if (slice_coord[i] != PARAMETER_UNSET) {
-      slabstart[2-i] = (h5size_t) ((slice_coord[i] - patch.origin[i]) /
+    if (slab_coord[i] != PARAMETER_UNSET) {
+      slabstart[2-i] = (h5size_t) ((slab_coord[i] - patch.origin[i]) /
                                    patch.delta[i] + 0.5);
       slabcount[2-i] = 1;
     }
@@ -454,17 +467,22 @@ static void ReadPatch (const patch_t& patch, int last_iteration)
              << patch.iorigin[0] + i*(1 << (max_rflevel - patch.rflevel)) << " "
              << patch.iorigin[1] + j*(1 << (max_rflevel - patch.rflevel)) << " "
              << patch.iorigin[2] + k*(1 << (max_rflevel - patch.rflevel)) << "\t"
+             << setiosflags (ios_base::fixed)
              << patch.time << "\t"
              << patch.origin[0] + i*patch.delta[0] << " "
              << patch.origin[1] + j*patch.delta[1] << " "
              << patch.origin[2] + k*patch.delta[2] << "\t"
-             << *data++ << endl;
+             << resetiosflags (ios_base::fixed)
+             << setiosflags (ios_base::scientific)
+             << *data++
+             << resetiosflags (ios_base::scientific)
+             << endl;
       }
-      if (slice_coord[0] == PARAMETER_UNSET) cout << endl;
+      if (slab_coord[0] == PARAMETER_UNSET) cout << endl;
     }
-    if (slice_coord[1] == PARAMETER_UNSET) cout << endl;
+    if (slab_coord[1] == PARAMETER_UNSET) cout << endl;
   }
-  if (slice_coord[2] == PARAMETER_UNSET) cout << endl;
+  if (slab_coord[2] == PARAMETER_UNSET) cout << endl;
 
   data -= npoints;
   delete[] data;
@@ -476,9 +494,9 @@ static void ReadPatch (const patch_t& patch, int last_iteration)
 static bool ComparePatches (const patch_t& a, const patch_t& b)
 {
   if (a.iteration != b.iteration) return (a.iteration < b.iteration);
-  if (a.rflevel != b.rflevel) return (a.rflevel < b.rflevel);
-  if (a.mglevel != b.mglevel) return (a.mglevel < b.mglevel);
-  if (a.map != b.map) return (a.map < b.map);
+  if (a.rflevel   != b.rflevel)   return (a.rflevel   < b.rflevel);
+  if (a.mglevel   != b.mglevel)   return (a.mglevel   < b.mglevel);
+  if (a.map       != b.map)       return (a.map       < b.map);
   if (a.component != b.component) return (a.component < b.component);
   return (a.timelevel < b.timelevel);
 }