diff options
Diffstat (limited to 'Carpet/CarpetIOHDF5/src/util/hdf5tobinary_slicer.cc')
| -rw-r--r-- | Carpet/CarpetIOHDF5/src/util/hdf5tobinary_slicer.cc | 814 |
1 files changed, 814 insertions, 0 deletions
diff --git a/Carpet/CarpetIOHDF5/src/util/hdf5tobinary_slicer.cc b/Carpet/CarpetIOHDF5/src/util/hdf5tobinary_slicer.cc new file mode 100644 index 000000000..d1d75c169 --- /dev/null +++ b/Carpet/CarpetIOHDF5/src/util/hdf5tobinary_slicer.cc @@ -0,0 +1,814 @@ + /*@@ + @file hdf5tobinary_slicer.cc + @date April 23, 2009 + @author Erik Schnetter, based on hdf5toascii_slicer.cc by Thomas Radke + @desc + Utility program to extract a 1D line or 2D slice from 3D datasets + in datafiles generated by CarpetIOHDF5. + The extracted slab is written to a binary file. + @enddesc + @@*/ + +#include <algorithm> +#include <iostream> +#include <iomanip> +#include <limits> +#include <sstream> +#include <string> +#include <vector> +#include <cassert> +#include <cmath> +#include <cstdio> +#include <cstring> +#include <regex.h> + +#define H5_USE_16_API 1 +#include <hdf5.h> + +using namespace std; + +/***************************************************************************** + ************************* Macro Definitions *************************** + *****************************************************************************/ + +// uncomment the following line to get some debug output +// #define DEBUG 1 + +// value for an unset parameter +#define PARAMETER_UNSET -424242.0 + +// fuzzy factor for comparing dataset timesteps with user-specified value +#define FUZZY_FACTOR 1e-6 + +// the datatype for the 'start' argument in H5Sselect_hyperslab() +#if (H5_VERS_MAJOR == 1 && \ + (H5_VERS_MINOR < 6 || (H5_VERS_MINOR == 6 && H5_VERS_RELEASE < 4))) +#define h5size_t hssize_t +#else +#define h5size_t hsize_t +#endif + +// macro to check the return code of calls to the HDF5 library +#define CHECK_HDF5(fn_call) { \ + const int _retval = fn_call; \ + if (_retval < 0) { \ + cerr << "HDF5 call " << #fn_call \ + << " in file " << __FILE__ << " line " << __LINE__ \ + << " returned error code " << _retval << endl; \ + } \ + } + + +/***************************************************************************** + ************************** Type Definitions *************************** + *****************************************************************************/ + +typedef struct { + hid_t file; + string name; + bool is_real; + hsize_t dims[3]; // C order + + int map, mglevel, component, iteration, timelevel, rflevel; + double time; + int iorigin[3]; // Fortran order + double origin[3], delta[3]; // Fortran order +} patch_t; + + +/***************************************************************************** + ************************* Global Data ************************* + *****************************************************************************/ + +// the slab coordinate as selected by the user (Fortan order) +static double slab_coord[3] = {PARAMETER_UNSET, PARAMETER_UNSET, PARAMETER_UNSET}; + +// flag for outputting data in full 3D +static bool out3D = false; + +// the specific timestep selected by the user +static double timestep = PARAMETER_UNSET; + +// the regular expression to match against dataset names +static const char* regex = NULL; +static regex_t preg; + +// the dataset name, if given explicitly +static const char* datasetname = NULL; + +// the list of all patches +static vector<patch_t> patchlist; + +// maximum refinement level across all patches +static int max_rflevel = -1; + +// whether to print omitted patches +static bool verbose = false; + +// pointer to array that is filled with the data (Fortran order) +static int datadims[3] = {-1,-1,-1}; +static float* dataptr = NULL; +static unsigned char* bdataptr = NULL; + +// output file name +static char* basename1 = NULL; + +// output format, either raw binary or HDF5 +static bool out_hdf5 = false; +static int out_hdf5_4d_index = -1; + +// output type (float or byte) +static bool out_float = false; +static float data_min = -1.0; +static float data_max = +1.0; + +/***************************************************************************** + ************************* Function Prototypes ************************* + *****************************************************************************/ +static herr_t FindPatches (hid_t group, const char *name, void *_file); +static void ReadPatch (const patch_t& patch, int last_iteration); +static bool ComparePatches (const patch_t& a, const patch_t& b); + + + /*@@ + @routine main + @date Sun 30 July 2006 + @author Thomas Radke + @desc + Evaluates command line options and opens the input files. + @enddesc + + @var argc + @vdesc number of command line parameters + @vtype int + @vio in + @endvar + @var argv + @vdesc array of command line parameters + @vtype char* const [] + @vio in + @endvar + @@*/ +int main (int argc, char *const argv[]) +{ + int i; + bool help = false; +#if 0 + int iorigin[3] = {-1, -1, -1}; +#endif + int num_slab_options = 0; + + + // evaluate command line parameters + for (i = 1; i < argc; i++) { + if (strcmp (argv[i], "--help") == 0) { + help = true; break; + } else if (strcmp (argv[i], "--verbose") == 0) { + verbose = true; + } else if (strcmp (argv[i], "--out3d-cube") == 0) { + out3D = true; + } else if (strcmp (argv[i], "--timestep") == 0 and i+1 < argc) { + timestep = atof (argv[++i]); + } else if (strcmp (argv[i], "--match") == 0 and i+1 < argc) { + regex = argv[++i]; + if (regcomp (&preg, regex, REG_EXTENDED)) { + cerr << "Error: invalid regular expression '" << regex << "' given" + << endl << endl; + return (-1); + } + } else if (strcmp (argv[i], "--dataset") == 0 and i+1 < argc) { + datasetname = argv[++i]; + } else if (strcmp (argv[i], "--out-precision") == 0 and i+1 < argc) { + cout << setprecision (atoi (argv[++i])); + } else if (strcmp (argv[i], "--out0d-point") == 0 and i+3 < argc) { + slab_coord[0] = atof (argv[++i]); + slab_coord[1] = atof (argv[++i]); + slab_coord[2] = atof (argv[++i]); num_slab_options++; + } 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) { + slab_coord[0] = atof (argv[++i]); num_slab_options++; + } else if (strcmp (argv[i], "--out2d-xzplane-y") == 0 and i+1 < argc) { + slab_coord[1] = atof (argv[++i]); num_slab_options++; + } else if (strcmp (argv[i], "--out2d-xyplane-z") == 0 and i+1 < argc) { + 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_slab_options++; + } else if (strcmp (argv[i], "--out2d-xzplane-yi") == 0 and i+1 < argc) { + 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_slab_options++; +#endif + } else if (strcmp (argv[i], "--size") == 0 and i+3 < argc) { + datadims[0] = atof (argv[++i]); + datadims[1] = atof (argv[++i]); + datadims[2] = atof (argv[++i]); + } else if (strcmp (argv[i], "--basename") == 0 and i+1 < argc) { + basename1 = argv[++i]; + } else if (strcmp (argv[i], "--out-hdf5") == 0) { + out_hdf5 = true; + } else if (strcmp (argv[i], "--out-float") == 0) { + out_float = true; + } else if (strcmp (argv[i], "--data-min") == 0 and i+1 < argc) { + data_min = atof (argv[++i]); + } else if (strcmp (argv[i], "--data-max") == 0 and i+1 < argc) { + data_max = atof (argv[++i]); + } else if (strcmp (argv[i], "--out-hdf5-4d-index") == 0 and i+1 < argc) { + out_hdf5_4d_index = atof (argv[++i]); + } else { + break; + } + } + + /* give some help if called with incorrect number of parameters */ + if (help or i >= argc or num_slab_options != (out3D ? 0 : 1)) { + const string indent (strlen (argv[0]) + 1, ' '); + cerr << endl << " ---------------------------" + << endl << " Carpet HDF5 to ASCII Slicer" + << endl << " ---------------------------" << endl + << endl + << "Usage: " << endl + << argv[0] << " [--help]" << endl + << indent << "[--match <regex string>]" << endl + << indent << "[--dataset <string (can contain %p)>]" << endl + << indent << "[--out-precision <digits>]" << endl + << indent << "[--timestep <cctk_time value>]" << endl + << indent << "[--verbose]" << endl + << indent << "<--out0d-point value value value> | <--out1d-line value value> | <--out2d-plane value> | <out3d-cube>" << endl + << indent << "--size <ni> <nj> <nk>" << endl + << indent << "--basename <filename>" << endl + << indent << "[--out-hdf5]" << endl + << indent << "[--out-float]" << endl + << indent << "[--out-hdf5-4d-index <t index>]" << endl + << indent << "<hdf5_infiles>" << endl << endl + << " where" << endl + << " [--help] prints this help" << endl + << " [--match <regex string>] selects HDF5 datasets by their names" << endl + << " matching a regex string using POSIX" << endl + << " Extended Regular Expression syntax" << 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 + << endl + << " and either --out0d-point value value value> or <--out1d-line value value> or <--out2d-plane value> or <--out3d-cube>" << endl + << " must be specified as in the following:" << endl + << " --out0d-point <origin_x> <origin_y> <origin_z>" << endl + << 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 + << endl + << " --out3d-cube" << endl +#if 0 + << " --out2d-yzplane-xi <origin_xi>" << endl + << " --out2d-xzplane-yi <origin_yi>" << endl + << " --out2d-xyplane-zi <origin_zi>" << endl +#endif + << endl + << " eg, to extract a 2D xy-plane at z = 0:" << endl + << " " << argv[0] << " --out2d-xyplane-z 0 alp.file_*.h5" << endl + << endl + << " or the same plane but only for datasets of refinement level 0:" << endl + << " " << argv[0] << " --match 'ADMBASE::alp it=[0-9]+ tl=0 rl=0' --out2d-xyplane-z 0 alp.file_*.h5" << endl; + return (0); + } + + if (not basename1) { + cerr << "Parameter --basename is not set" << endl; + return 1; + } + + if (regex and datasetname) { + cerr << "Parameters --match and --dataset cannot be used together" << endl; + return 1; + } + + { + if (datadims[0]<0 || datadims[1]<0 || datadims[2]<0) { + cerr << "Parameter --size is not set" << endl; + return 1; + } + size_t const npoints = (size_t)datadims[0]*datadims[1]*datadims[2]; + dataptr = new float[npoints]; + assert (dataptr); + bdataptr = new unsigned char[npoints]; + assert (bdataptr); + } + + vector<hid_t> filelist; + for (; i < argc; i++) { + hid_t file; + + H5E_BEGIN_TRY { + file = H5Fopen (argv[i], H5F_ACC_RDONLY, H5P_DEFAULT); + } H5E_END_TRY; + + if (file < 0) { + cerr << "Could not open Carpet HDF5 input file '" << argv[i] << "'" + << endl << endl; + continue; + } + + if (not datasetname) { + cout << "File " << argv[i] << ":" << endl; + size_t const before = patchlist.size(); + CHECK_HDF5 (H5Giterate (file, "/", NULL, FindPatches, &file)); + size_t const after = patchlist.size(); + size_t const count = after - before; + cout << " will read " << count << " " + << (count==1 ? "patch" : "patches") << endl; + } else { + const string ppmarker ("%p"); + const string filemarker (".file_"); + string name (datasetname); + const size_t pppos = name.find(ppmarker); + if (pppos != string::npos) { + const string filename (argv[i]); + const size_t procpos = filename.find(filemarker); + assert (procpos != string::npos); + const int proc = atoi(filename.c_str() + procpos + filemarker.length()); + ostringstream namebuf; + namebuf << name.substr(0,pppos) << proc + << name.substr(pppos + ppmarker.length()); + name = namebuf.str(); + } + cout << " Examining dataset " << name << endl; + hid_t group; + CHECK_HDF5 (group = H5Dopen (file, name.c_str())); + FindPatches (group, name.c_str(), &file); + CHECK_HDF5 (H5Dclose (group)); + } + + filelist.push_back (file); + } + + if (! patchlist.empty()) { + /* now sort the patches by their time attribute */ + sort (patchlist.begin(), patchlist.end(), ComparePatches); + + cout << "# 1D ASCII output created by"; + for (int j = 0; j < argc; j++) cout << " " << argv[j]; + cout << endl << "#" << endl; + + int last_iteration = patchlist[0].iteration - 1; + for (size_t j = 0; j < patchlist.size(); j++) { + ReadPatch (patchlist[j], last_iteration); + last_iteration = patchlist[j].iteration; + } + } else { + cerr << "No valid datasets found" << endl << endl; + } + + // find minimum and maximum + { + float const float_max = numeric_limits<float>::max(); + float data_min1 = +float_max; + float data_max1 = -float_max; + double data_sum = 0.0; + double data_sum2 = 0.0; + double data_count = 0.0; + size_t const npoints = (size_t)datadims[0]*datadims[1]*datadims[2]; + for (size_t n=0; n<npoints; ++n) { + float const data = dataptr[n]; + if (data > -sqrt(float_max) and data < sqrt(float_max)) { + data_min1 = min(data_min1, data); + data_max1 = max(data_max1, data); + data_sum += (double)data; + data_sum2 += pow((double)data, 2.0); + ++ data_count; + } + } + double const data_avg = data_sum / data_count; + double const data_stddev = + sqrt(max(0.0, data_sum2/data_count - pow(data_avg, 2))); + + cerr << "Minimum: " << data_min1 << endl + << "Maximum: " << data_max1 << endl + << "Average: " << data_avg << endl + << "Standard deviation: " << data_stddev << endl + << "Valid points: " << data_count << endl; + } + + // convert from float to byte + if (not out_float) { + unsigned char const byte_min = 0; + unsigned char const byte_max = 255; + float const scale_factor = + ((float)byte_max - (float)byte_min + 1) / (data_max - data_min); + size_t const npoints = (size_t)datadims[0]*datadims[1]*datadims[2]; + for (size_t n=0; n<npoints; ++n) { + float const data = dataptr[n]; + float const data_scaled = + (float)byte_min + floor((data - data_min) * scale_factor); + float const data_clamped = + max((float)byte_min, min((float)byte_max, data_scaled)); + unsigned char const bdata = + (unsigned char)data_clamped; + bdataptr[n] = bdata; + } + } + + hid_t const out_type = out_float ? H5T_NATIVE_FLOAT : H5T_NATIVE_UCHAR; + size_t const out_size = out_float ? sizeof *dataptr : sizeof *bdataptr; + void const * const out_data = + out_float ? + static_cast<void const*>(dataptr) : + static_cast<void const*>(bdataptr); + + // write output files + if (out_hdf5) { + // HDF5 output + if (out_hdf5_4d_index == -1) { + // 3D output + char filename[1000]; + snprintf (filename, sizeof filename, "%s.h5", basename1); + hid_t const outfile = H5Fcreate (filename, H5F_ACC_TRUNC, + H5P_DEFAULT, H5P_DEFAULT); + assert (outfile >= 0); + hsize_t const dims[3] = { datadims[2], datadims[1], datadims[0] }; + hid_t const dataspace = H5Screate_simple (3, dims, NULL); + assert (dataspace >= 0); + hid_t const dataset_properties = H5Pcreate (H5P_DATASET_CREATE); + assert (dataset_properties >= 0); + hid_t const dataset = H5Dcreate (outfile, "data", + out_type, dataspace, + dataset_properties); + assert (dataset >= 0); + hid_t const transfer_properties = H5Pcreate (H5P_DATASET_XFER); + assert (transfer_properties >= 0); + { + herr_t const herr = H5Dwrite (dataset, out_type, + H5S_ALL, H5S_ALL, + transfer_properties, out_data); + assert (herr >= 0); + } + CHECK_HDF5 (H5Dclose (dataset)); + CHECK_HDF5 (H5Sclose (dataspace)); + CHECK_HDF5 (H5Fclose (outfile)); + } else { + // 4d output + char filename[1000]; + snprintf (filename, sizeof filename, "%s.h5", basename1); + htri_t is_hdf5; + H5E_BEGIN_TRY { + is_hdf5 = H5Fis_hdf5 (filename); + } H5E_END_TRY; + bool const file_is_new = not (is_hdf5 > 0); + hid_t const outfile = + file_is_new ? + H5Fcreate (filename, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT) : + H5Fopen (filename, H5F_ACC_RDWR, H5P_DEFAULT); + assert (outfile >= 0); + + hsize_t const dims[4] = + { 0, datadims[2], datadims[1], datadims[0] }; + hsize_t const maxdims[4] = + { H5S_UNLIMITED, datadims[2], datadims[1], datadims[0] }; + hid_t const dataspace = H5Screate_simple (4, dims, maxdims); + assert (dataspace >= 0); + hid_t const dataset_properties = H5Pcreate (H5P_DATASET_CREATE); + assert (dataset_properties >= 0); + // hsize_t const chunk_dims[4] = {1, 32, 32, 32 }; + // Chunks must be < 4 GByte + // hsize_t const chunk_dims[4] = + // {1, datadims[2], datadims[1], datadims[0] }; + hsize_t const chunk_dims[4] = + {1, min(128,datadims[2]), min(128,datadims[1]), min(128,datadims[0]) }; + CHECK_HDF5 (H5Pset_chunk (dataset_properties, 4, chunk_dims)); + int const compression_level = 1; + CHECK_HDF5 (H5Pset_deflate (dataset_properties, compression_level)); + hid_t const dataset = + file_is_new ? + H5Dcreate (outfile, "data", + out_type, dataspace, dataset_properties) : + H5Dopen (outfile, "data"); + assert (dataset >= 0); + // TODO: assert that the dataset's dataspace has the correct + // dimensions + + hsize_t const extended_dims[4] = + { out_hdf5_4d_index+1, datadims[2], datadims[1], datadims[0] }; + CHECK_HDF5 (H5Dextend (dataset, extended_dims)); + + // for debugging only + // CHECK_HDF5 (H5Fflush (outfile, H5F_SCOPE_GLOBAL)); + + hsize_t const memory_dims[4] = + { 1, datadims[2], datadims[1], datadims[0] }; + hid_t const memory_dataspace = H5Screate_simple (4, memory_dims, NULL); + assert (memory_dataspace >= 0); + + hid_t const file_dataspace = H5Screate_simple (4, extended_dims, NULL); + assert (file_dataspace >= 0); + // hssize_t const file_offset[4] = + // { out_hdf5_4d_index, 0, 0, 0 }; + // CHECK_HDF5 (H5Soffset_simple (file_dataspace, file_offset)); + hsize_t const file_offset[4] = { out_hdf5_4d_index, 0, 0, 0 }; + CHECK_HDF5 (H5Sselect_hyperslab (file_dataspace, H5S_SELECT_SET, + file_offset, NULL, memory_dims, NULL)); + + hid_t const transfer_properties = H5Pcreate (H5P_DATASET_XFER); + assert (transfer_properties >= 0); + CHECK_HDF5 (H5Dwrite (dataset, + out_type, memory_dataspace, file_dataspace, + transfer_properties, out_data)); + + CHECK_HDF5 (H5Sclose (file_dataspace)); + CHECK_HDF5 (H5Sclose (memory_dataspace)); + CHECK_HDF5 (H5Dclose (dataset)); + CHECK_HDF5 (H5Sclose (dataspace)); + CHECK_HDF5 (H5Fclose (outfile)); + } + } else { + // raw binary output + { + char filename[1000]; + snprintf (filename, sizeof filename, "%s.bin", basename1); + FILE* const outfile = fopen(filename, "w"); + assert (outfile); + size_t const npoints = (size_t)datadims[0]*datadims[1]*datadims[2]; + size_t const icnt = fwrite (out_data, out_size, npoints, outfile); + assert (icnt == npoints); + int const ierr = fclose (outfile); + assert (not ierr); + } + } + + { + delete[] dataptr; + dataptr = NULL; + delete[] bdataptr; + bdataptr = NULL; + } + + // close all input files + for (size_t j = 0; j < filelist.size(); j++) { + if (filelist[j] >= 0) { + CHECK_HDF5 (H5Fclose (filelist[j])); + } + } + + return (0); +} + + + /*@@ + @routine FindPatches + @date Tue 8 August 2006 + @author Thomas Radke + @desc + The worker routine which is called by H5Giterate(). + It checks whether the current HDF5 object is a patch matching + the user's slab criteria, and adds it to the patches list. + @enddesc + + @var group + @vdesc HDF5 object to start the iteration + @vtype hid_t + @vio in + @endvar + @var name + @vdesc name of the object at the current iteration + @vtype const char * + @vio in + @endvar + @var _file + @vdesc pointer to the descriptor of the currently opened file + @vtype void * + @vio in + @endvar +@@*/ +static herr_t FindPatches (hid_t group, const char *name, void *_file) +{ + int ndims; + hid_t dataset, datatype, attr; + H5G_stat_t object_info; + H5T_class_t typeclass; + patch_t patch; + + + /* we are interested in datasets only - skip anything else */ + CHECK_HDF5 (H5Gget_objinfo (group, name, 0, &object_info)); + if (object_info.type != H5G_DATASET) { + return (0); + } + + /* check the dataset's datatype - make sure it is either integer or real */ + CHECK_HDF5 (dataset = H5Dopen (group, name)); + patch.name = name; + + CHECK_HDF5 (datatype = H5Dget_type (dataset)); + CHECK_HDF5 (typeclass = H5Tget_class (datatype)); + CHECK_HDF5 (H5Tclose (datatype)); + + // read the timestep + CHECK_HDF5 (attr = H5Aopen_name (dataset, "time")); + CHECK_HDF5 (H5Aread (attr, H5T_NATIVE_DOUBLE, &patch.time)); + CHECK_HDF5 (H5Aclose (attr)); + + // read the dimensions + hid_t dataspace = -1; + CHECK_HDF5 (dataspace = H5Dget_space (dataset)); + ndims = H5Sget_simple_extent_ndims (dataspace); + + bool is_okay = false; + if (typeclass != H5T_FLOAT and typeclass != H5T_INTEGER) { + if (verbose) { + cerr << "skipping dataset '" << name << "':" << endl + << " is not of integer or floating-point datatype" << endl; + } + } else if (ndims != 3) { + if (verbose) { + cerr << "skipping dataset '" << name << "':" << endl + << " dataset has " << ndims << " dimensions" << endl; + } + } else if (regex && regexec (&preg, name, 0, NULL, 0)) { + if (verbose) { + cerr << "skipping dataset '" << name << "':" << endl + << " name doesn't match regex" << endl; + } + } else if (timestep != PARAMETER_UNSET and + fabs (timestep - patch.time) > FUZZY_FACTOR) { + if (verbose) { + cerr << "skipping dataset '" << name << "':" << endl + << " timestep (" << patch.time << ") doesn't match" << endl; + } + } else { + if (not out3D) { + CHECK_HDF5 (attr = H5Aopen_name (dataset, "origin")); + CHECK_HDF5 (H5Aread (attr, H5T_NATIVE_DOUBLE, patch.origin)); + CHECK_HDF5 (H5Aclose (attr)); + CHECK_HDF5 (attr = H5Aopen_name (dataset, "delta")); + CHECK_HDF5 (H5Aread (attr, H5T_NATIVE_DOUBLE, patch.delta)); + CHECK_HDF5 (H5Aclose (attr)); + } + CHECK_HDF5 (H5Sget_simple_extent_dims (dataspace, patch.dims, NULL)); + + int i; + is_okay = out3D; + if (not out3D) { + for (i = 0; i < 3; 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]; + if (not is_okay) break; + } + } + } + if (not is_okay) { + if (verbose) { + cerr << "skipping dataset '" << name << "':" << endl + << " slab " << slab_coord[i] << " is out of dataset range [" + << patch.origin[i] << ", " + << patch.origin[i] + (patch.dims[2-i]-1)*patch.delta[i] << "]" + << endl; + } + } + } + CHECK_HDF5 (H5Sclose (dataspace)); + + if (not is_okay) { + CHECK_HDF5 (H5Dclose (dataset)); + return (0); + } + + patch.file = *(hid_t *) _file; + patch.is_real = typeclass == H5T_FLOAT; + + /* read attributes */ + CHECK_HDF5 (attr = H5Aopen_name (dataset, "timestep")); + CHECK_HDF5 (H5Aread (attr, H5T_NATIVE_INT, &patch.iteration)); + CHECK_HDF5 (H5Aclose (attr)); + CHECK_HDF5 (attr = H5Aopen_name (dataset, "level")); + CHECK_HDF5 (H5Aread (attr, H5T_NATIVE_INT, &patch.rflevel)); + CHECK_HDF5 (H5Aclose (attr)); + CHECK_HDF5 (attr = H5Aopen_name (dataset, "group_timelevel")); + CHECK_HDF5 (H5Aread (attr, H5T_NATIVE_INT, &patch.timelevel)); + CHECK_HDF5 (H5Aclose (attr)); + CHECK_HDF5 (attr = H5Aopen_name (dataset, "carpet_mglevel")); + CHECK_HDF5 (H5Aread (attr, H5T_NATIVE_INT, &patch.mglevel)); + CHECK_HDF5 (H5Aclose (attr)); + CHECK_HDF5 (attr = H5Aopen_name (dataset, "iorigin")); + CHECK_HDF5 (H5Aread (attr, H5T_NATIVE_INT, patch.iorigin)); + CHECK_HDF5 (H5Aclose (attr)); + + CHECK_HDF5 (H5Dclose (dataset)); + + // the component and map info are not contained in the HDF5 dataset + // and thus have to be parsed from the dataset's name + patch.map = patch.component = 0; + string::size_type loc = patch.name.find (" m=", 0); + if (loc != string::npos) patch.map = atoi (patch.name.c_str() + loc + 3); + loc = patch.name.find (" c=", 0); + if (loc != string::npos) patch.component = atoi(patch.name.c_str() + loc + 3); + + if (max_rflevel < patch.rflevel) max_rflevel = patch.rflevel; + + patchlist.push_back (patch); + + return (0); +} + + + /*@@ + @routine ReadPatch + @date Tue 8 August 2006 + @author Thomas Radke + @desc + Reads a slab from the given patch + and outputs it in CarpetIOASCII format. + @enddesc + + @var patch + @vdesc patch to output + @vtype const patch_t& + @vio in + @endvar + @var last_iteration + @vdesc iteration number of the last patch that was output + @vtype int + @vio in + @endvar +@@*/ +static void ReadPatch (const patch_t& patch, int last_iteration) +{ + if (patch.iteration != last_iteration) cout << endl; + + cout << "# iteration " << patch.iteration << endl + << "# refinement level " << patch.rflevel + << " multigrid level " << patch.mglevel + << " map " << patch.map + << " component " << patch.component + << " time level " << patch.timelevel + << endl; + + // C order + h5size_t slabstart[3] = + {patch.iorigin[2], patch.iorigin[1], patch.iorigin[0]}; + hsize_t slabcount[3] = {patch.dims[0], patch.dims[1], patch.dims[2]}; + hsize_t slabsize[3] = {datadims[2], datadims[1], datadims[0]}; + for (int d=0; d<3; ++d) { + if (slab_coord[d] != PARAMETER_UNSET) { + slabstart[2-d] = (h5size_t) ((slab_coord[d] - patch.origin[d]) / + patch.delta[d] + 0.5); + slabcount[2-d] = 1; + } + } + for (int d=0; d<3; ++d) { + cout << "slab[" << d << "]: lbnd=" << slabstart[2-d] << " lsh=" << slabcount[2-d] << " gsh=" << slabsize[2-d] << endl; + } + for (int d=0; d<3; ++d) { + assert (slabstart[2-d] >= 0); + assert (slabcount[2-d] >= 0); + assert (slabstart[2-d] + slabcount[2-d] <= slabsize[2-d]); + } + + hid_t dataset, filespace, slabspace; + CHECK_HDF5 (dataset = H5Dopen (patch.file, patch.name.c_str())); + CHECK_HDF5 (filespace = H5Dget_space (dataset)); +#if 0 + CHECK_HDF5 (slabspace = H5Screate_simple (3, slabcount, slabsize)); + CHECK_HDF5 (H5Soffset_simple (slabspace, (hssize_t*)slabstart)); +// CHECK_HDF5 (H5Sselect_hyperslab (filespace, H5S_SELECT_SET, +// slabstart, NULL, slabcount, NULL)); +#endif + CHECK_HDF5 (slabspace = H5Screate_simple (3, slabsize, NULL)); + CHECK_HDF5 (H5Sselect_hyperslab (slabspace, H5S_SELECT_SET, + slabstart, NULL, slabcount, NULL)); + assert (patch.is_real); + CHECK_HDF5 (H5Dread (dataset, + H5T_NATIVE_FLOAT, + slabspace, filespace, H5P_DEFAULT, + dataptr)); + CHECK_HDF5 (H5Sclose (slabspace)); + CHECK_HDF5 (H5Sclose (filespace)); + CHECK_HDF5 (H5Dclose (dataset)); +} + + +// comparison function to sort patches by iteration number, refinement level, +// multigrid level, map, and component +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.component != b.component) return (a.component < b.component); + return (a.timelevel < b.timelevel); +} |