]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add HDF5 classes
authorDaniel Garcia-Sanchez <daniel.garcia-sanchez@insp.upmc.fr>
Thu, 2 Aug 2018 11:09:54 +0000 (13:09 +0200)
committerDaniel Garcia-Sanchez <daniel.garcia-sanchez@insp.upmc.fr>
Wed, 9 Jan 2019 17:08:36 +0000 (18:08 +0100)
include/deal.II/base/hdf5.h [new file with mode: 0644]
source/base/CMakeLists.txt
source/base/hdf5.cc [new file with mode: 0644]

diff --git a/include/deal.II/base/hdf5.h b/include/deal.II/base/hdf5.h
new file mode 100644 (file)
index 0000000..08e9c98
--- /dev/null
@@ -0,0 +1,363 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 1998 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+#ifndef dealii_hdf5_h
+#define dealii_hdf5_h
+
+#include <deal.II/base/config.h>
+
+#include <deal.II/lac/full_matrix.h>
+
+#include <hdf5.h>
+
+#include <vector>
+
+DEAL_II_NAMESPACE_OPEN
+
+/**
+ * Namespace containing the HDF5 wrappers.
+ *
+ * This set of classes can be used to store the data in the hdf5 file format.
+ *
+ * The classes can be used to improve the interaction with python scripts.
+ *
+ * This python script writes the parameters for a deal.ii simulation:
+ * @code
+ * h5_file = h5py.File('simulation.hdf5','w')
+ * data = h5_file.create_group('data')
+ * data.attrs['nb_frequency_points'] = 50 # int
+ * data.attrs['nb_slice_points'] = 10 # int
+ * data.attrs['rho'] = 2300.5 # double
+ * data.attrs['simulation_type'] = 'elastic_equation' # utf8 string
+ * @endcode
+ *
+ * C++ deal.ii simulation with MPI HDF5:
+ * @code
+ * hdf5::File data_file("simulation.hdf5", MPI_COMM_WORLD,
+ *                      hdf5::File::Mode::open);
+ * hdf5::Group data = data_file.group("data");
+ *
+ * // Read some parameters
+ * const int nb_frequency_points = data.attr<int>("nb_frequency_points");
+ * const int nb_slice_points = data.attr<int>("nb_cut_points");
+ * const double rho = data.attr<double>("rho");
+ * const std::string simulation_type =
+ * data.attr<std::string>("simulation_type");
+ *
+ * // Create some distributed datasets
+ * hdf5::DataSet<double> frequency_dataset(data.template
+ *                          create_dataset<double>("frequency",
+ *                          std::vector<hsize_t>{nb_frequency_points}));
+ *
+ * hdf5::DataSet<std::complex<double>> displacement(parameters.data.template
+ * create_dataset<std::complex<double>>("displacement",
+ *                                      std::vector<hsize_t>{nb_slice_points,
+ *                                      nb_frequency_points}));
+ *
+ * // ... do some calculations ...
+ * // The displacement data is distributed across the MPI processes.
+ * // The displacement_data vector is local. The coordinates vector corresponds
+ * // to the coordinates of the data in the hdf5 dataset.
+ * if (coordinates.size() > 0) {
+ *   displacement.write_data_selection(displacement_data, coordinates);
+ * } else {
+ *   // Even if this MPI process has no displacement data to write,
+ *   // it has to contribute to the collective call
+ *   displacement.write_data_none();
+ * }
+ *
+ * // ... do some calculations ...
+ * // The whole frequency dataset can be written by one single MPI process,
+ * // but the rest of the processes have to contribute to the collective call
+ * if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0) {
+ *   frequency_dataset.write_data(frequency);
+ * } else {
+ *   frequency_dataset.write_data_none();
+ * }
+ *
+ * // Write the simulation metadata
+ * data.write_attr("active_cells", triangulation.n_active_cells());
+ * data.write_attr("degrees_of_freedom", dof_handler.n_dofs());
+ * @endcode
+ *
+ * Read the simulation data with python:
+ * @code
+ * h5_file = h5py.File('simulation.hdf5','r+')
+ * data = h5_file['data']
+ * frequency = data['frequency'] # float64 dtype
+ * displacement = data['displacement'] # complex128 dtype
+ * print('Degrees of freedom: ' + repr(data.attrs['degrees_of_freedom']))
+ * @endcode
+ *
+ * @author Daniel Garcia-Sanchez, 2018
+ */
+namespace hdf5
+
+{
+  /**
+   * General class for the HDF5 objects.
+   *
+   * @author Daniel Garcia-Sanchez, 2018
+   */
+  class HDF5Object
+  {
+  protected:
+    HDF5Object(const std::string name, bool mpi);
+
+  public:
+    enum class Mode
+    {
+      create,
+      open
+    };
+
+    /**
+     * Reads an attribute. T can be double, int, unsigned int, bool,
+     * std::complex<double> or std::string. Note that the encoding of
+     * std::string is UTF8 in order to be compatible with python3.
+     *
+     * Datatype conversion takes place at the time of a read or write and is
+     * automatic. See the <a
+     * href="https://support.hdfgroup.org/HDF5/doc/UG/HDF5_Users_Guide-Responsive%20HTML5/index.html#t=HDF5_Users_Guide%2FDatatypes%2FHDF5_Datatypes.htm%23TOC_6_10_Data_Transferbc-26&rhtocid=6.5_2">"Data
+     * Transfer: Datatype Conversion and Selection"</a>  section in the HDF5
+     * User's Guide.
+     */
+    template <typename T>
+    T
+    attr(const std::string attr_name) const;
+
+    /**
+     * Writes an attribute. T can be double, int, unsigned int, bool,
+     * std::complex<double> or std::string. Note that the encoding of
+     * std::string is UTF8 in order to be compatible with python3.
+     *
+     * Datatype conversion takes place at the time of a read or write and is
+     * automatic. See the <a
+     * href="https://support.hdfgroup.org/HDF5/doc/UG/HDF5_Users_Guide-Responsive%20HTML5/index.html#t=HDF5_Users_Guide%2FDatatypes%2FHDF5_Datatypes.htm%23TOC_6_10_Data_Transferbc-26&rhtocid=6.5_2">"Data
+     * Transfer: Datatype Conversion and Selection"</a>  section in the HDF5
+     * User's Guide.
+     */
+    template <typename T>
+    void
+                      write_attr(const std::string attr_name, T value) const;
+    const std::string name;
+
+  protected:
+    std::shared_ptr<hid_t> hdf5_reference;
+    // If true use mpi hdf5. If false use serial hdf5.
+    const bool mpi;
+  };
+
+  /**
+   * This class implements a HDF5 DataSet.
+   *
+   * @author Daniel Garcia-Sanchez, 2018
+   */
+  template <typename T>
+  class DataSet : public HDF5Object
+  {
+    friend class Group;
+
+  protected:
+    DataSet(const std::string    name,
+            const hid_t &        parent_group_id,
+            std::vector<hsize_t> dimensions,
+            bool                 mpi,
+            const Mode           mode);
+
+  public:
+    ~DataSet();
+
+    /**
+     * Writes data in the dataset. T can be double, int, unsigned int, bool
+     * or std::complex<double>.
+     *
+     * Datatype conversion takes place at the time of a read or write and is
+     * automatic. See the <a
+     * href="https://support.hdfgroup.org/HDF5/doc/UG/HDF5_Users_Guide-Responsive%20HTML5/index.html#t=HDF5_Users_Guide%2FDatatypes%2FHDF5_Datatypes.htm%23TOC_6_10_Data_Transferbc-26&rhtocid=6.5_2">"Data
+     * Transfer: Datatype Conversion and Selection"</a>  section in the HDF5
+     * User's Guide.
+     */
+    void
+    write_data(const std::vector<T> &data) const;
+
+    /**
+     * Writes data in the dataset. T can be double, int, unsigned int, bool
+     * or std::complex<double>.
+     *
+     * Datatype conversion takes place at the time of a read or write and is
+     * automatic. See the <a
+     * href="https://support.hdfgroup.org/HDF5/doc/UG/HDF5_Users_Guide-Responsive%20HTML5/index.html#t=HDF5_Users_Guide%2FDatatypes%2FHDF5_Datatypes.htm%23TOC_6_10_Data_Transferbc-26&rhtocid=6.5_2">"Data
+     * Transfer: Datatype Conversion and Selection"</a>  section in the HDF5
+     * User's Guide.
+     */
+    void
+    write_data(const dealii::FullMatrix<T> &data) const;
+
+    /**
+     * Writes data to a subset of the dataset. T can be double, int, unsigned
+     * int, bool or std::complex<double>.
+     *
+     * The selected elements can be scattered and take any shape in the dataset.
+     * For examples, in the 4D case, we will be selecting three points and a 4D
+     * dataspace has rank 4, so the selection will be described in a 3-by-4
+     * array. To select the points (1,1,1,1), (14,6,12,18), and (8,22,30,22),
+     * the point selection array would be as follows:
+     *
+     *    0  0  0  0
+     *
+     *   13  5 11 17
+     *
+     *    7 21 29 21
+     *
+     * Datatype conversion takes place at the time of a read or write and is
+     * automatic. See the <a
+     * href="https://support.hdfgroup.org/HDF5/doc/UG/HDF5_Users_Guide-Responsive%20HTML5/index.html#t=HDF5_Users_Guide%2FDatatypes%2FHDF5_Datatypes.htm%23TOC_6_10_Data_Transferbc-26&rhtocid=6.5_2">"Data
+     * Transfer: Datatype Conversion and Selection"</a>  section in the HDF5
+     * User's Guide.
+     */
+    void
+    write_data_selection(const std::vector<T> &     data,
+                         const std::vector<hsize_t> coordinates) const;
+
+    /**
+     * Writes data to a subset of the dataset. T can be double, int, unsigned
+     * int, bool or std::complex<double>.
+     *
+     * The selected elements form a hyperslab in the dataset.
+     *
+     * Datatype conversion takes place at the time of a read or write and is
+     * automatic. See the <a
+     * href="https://support.hdfgroup.org/HDF5/doc/UG/HDF5_Users_Guide-Responsive%20HTML5/index.html#t=HDF5_Users_Guide%2FDatatypes%2FHDF5_Datatypes.htm%23TOC_6_10_Data_Transferbc-26&rhtocid=6.5_2">"Data
+     * Transfer: Datatype Conversion and Selection"</a>  section in the HDF5
+     * User's Guide.
+     */
+    void
+    write_data_hyperslab(const std::vector<T> &     data,
+                         const std::vector<hsize_t> offset,
+                         const std::vector<hsize_t> count) const;
+
+    /**
+     * Writes data to a subset of the dataset. T can be double, int, unsigned
+     * int, bool or std::complex<double>.
+     *
+     * The selected elements form a hyperslab in the dataset.
+     *
+     * Datatype conversion takes place at the time of a read or write and is
+     * automatic. See the <a
+     * href="https://support.hdfgroup.org/HDF5/doc/UG/HDF5_Users_Guide-Responsive%20HTML5/index.html#t=HDF5_Users_Guide%2FDatatypes%2FHDF5_Datatypes.htm%23TOC_6_10_Data_Transferbc-26&rhtocid=6.5_2">"Data
+     * Transfer: Datatype Conversion and Selection"</a>  section in the HDF5
+     * User's Guide.
+     */
+    void
+    write_data_hyperslab(const dealii::FullMatrix<T> &data,
+                         const std::vector<hsize_t>   offset,
+                         const std::vector<hsize_t>   count) const;
+    void
+                               write_data_none() const;
+    const unsigned int         rank;
+    const std::vector<hsize_t> dimensions;
+
+  private:
+    std::shared_ptr<hid_t> dataspace;
+    std::shared_ptr<hid_t> t_type;
+    unsigned int           total_size;
+  };
+
+  /**
+   * This class implements a HDF5 Group
+   *
+   * @author Daniel Garcia-Sanchez, 2018
+   */
+  class Group : public HDF5Object
+  {
+  protected:
+    Group(const std::string name,
+          const Group &     parent_group,
+          const bool        mpi,
+          const Mode        mode);
+    Group(const std::string name, const bool mpi);
+
+  public:
+    Group
+    group(const std::string name);
+
+    /**
+     * Creates a group.
+     */
+    Group
+    create_group(const std::string name);
+
+    /**
+     * Creates and writes data to a dataset. T can be double, int, unsigned int,
+     * bool or std::complex<double>.
+     *
+     * Datatype conversion takes place at the time of a read or write and is
+     * automatic. See the <a
+     * href="https://support.hdfgroup.org/HDF5/doc/UG/HDF5_Users_Guide-Responsive%20HTML5/index.html#t=HDF5_Users_Guide%2FDatatypes%2FHDF5_Datatypes.htm%23TOC_6_10_Data_Transferbc-26&rhtocid=6.5_2">"Data
+     * Transfer: Datatype Conversion and Selection"</a>  section in the HDF5
+     * User's Guide.
+     */
+    template <typename T>
+    void
+    write_dataset(const std::string name, const T &data) const;
+    template <typename T>
+
+    /**
+     * Creates a dataset. T can be double, int, unsigned int, bool or
+     * std::complex<double>.
+     *
+     * Datatype conversion takes place at the time of a read or write and is
+     * automatic. See the <a
+     * href="https://support.hdfgroup.org/HDF5/doc/UG/HDF5_Users_Guide-Responsive%20HTML5/index.html#t=HDF5_Users_Guide%2FDatatypes%2FHDF5_Datatypes.htm%23TOC_6_10_Data_Transferbc-26&rhtocid=6.5_2">"Data
+     * Transfer: Datatype Conversion and Selection"</a>  section in the HDF5
+     * User's Guide.
+     */
+    DataSet<T>
+    create_dataset(const std::string          name,
+                   const std::vector<hsize_t> dimensions) const;
+  };
+
+  /**
+   * This class implements a HDF5 Group
+   *
+   * @author Daniel Garcia-Sanchez, 2018
+   */
+  class File : public Group
+  {
+  private:
+    File(const std::string name,
+         const bool        mpi,
+         const MPI_Comm    mpi_communicator,
+         const Mode        mode);
+
+  public:
+    /**
+     * Creates or opens a hdf5 file.
+     */
+    File(const std::string name,
+         const MPI_Comm    mpi_communicator,
+         const Mode        mode = Mode::create);
+
+    /**
+     * Creates or opens a hdf5 file in parallel.
+     */
+    File(const std::string name, const Mode mode = Mode::create);
+  };
+} // namespace hdf5
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
index 1379d6b6e226b749d38f7e85ac2d0592382900f6..34b1e95033d8078c126e895d00d3a4922b555ccf 100644 (file)
@@ -43,6 +43,7 @@ SET(_unity_include_src
   index_set.cc
   job_identifier.cc
   logstream.cc
+  hdf5.cc
   mpi.cc
   multithread_info.cc
   named_selection.cc
diff --git a/source/base/hdf5.cc b/source/base/hdf5.cc
new file mode 100644 (file)
index 0000000..73384b4
--- /dev/null
@@ -0,0 +1,732 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2017 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+#include <deal.II/base/config.h>
+
+#include <deal.II/base/hdf5.h>
+
+#include <deal.II/lac/full_matrix.h>
+
+#include <hdf5.h>
+
+#include <memory>
+#include <numeric>
+#include <vector>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace hdf5
+
+{
+  HDF5Object::HDF5Object(const std::string name, bool mpi)
+    : name(name)
+    , mpi(mpi)
+  {}
+
+  template <typename T>
+  T
+  HDF5Object::attr(const std::string attr_name) const
+  {
+    std::shared_ptr<hid_t> t_type;
+    T                      value;
+    hid_t                  attr;
+    if (std::is_same<T, double>::value)
+      {
+        t_type  = std::shared_ptr<hid_t>(new hid_t);
+        *t_type = H5T_NATIVE_DOUBLE;
+      }
+    else if (std::is_same<T, int>::value)
+      {
+        t_type  = std::shared_ptr<hid_t>(new hid_t);
+        *t_type = H5T_NATIVE_INT;
+      }
+    else if (std::is_same<T, unsigned int>::value)
+      {
+        t_type  = std::shared_ptr<hid_t>(new hid_t);
+        *t_type = H5T_NATIVE_UINT;
+      }
+    else if (std::is_same<T, std::complex<double>>::value)
+      {
+        t_type  = std::shared_ptr<hid_t>(new hid_t, [](auto pointer) {
+          // Relase the HDF5 resource
+          H5Tclose(*pointer);
+          delete pointer;
+        });
+        *t_type = H5Tcreate(H5T_COMPOUND, sizeof(std::complex<double>));
+        //  The C++ standards committee agreed to mandate that the storage
+        //  format used for the std::complex type be binary-compatible with the
+        //  C99 type, i.e. an array T[2] with consecutive real [0] and imaginary
+        //  [1] parts.
+        H5Tinsert(*t_type, "r", 0, H5T_NATIVE_DOUBLE);
+        H5Tinsert(*t_type, "i", sizeof(double), H5T_NATIVE_DOUBLE);
+      }
+    else
+      {
+        Assert(false, ExcInternalError());
+      }
+
+    attr = H5Aopen(*hdf5_reference, attr_name.data(), H5P_DEFAULT);
+    H5Aread(attr, *t_type, &value);
+    H5Aclose(attr);
+    return value;
+  }
+
+  template <>
+  bool
+  HDF5Object::attr(const std::string attr_name) const
+  {
+    // The enum field generated by h5py can be casted to int
+    int   int_value;
+    hid_t attr;
+    attr = H5Aopen(*hdf5_reference, attr_name.data(), H5P_DEFAULT);
+    H5Aread(attr, H5T_NATIVE_INT, &int_value);
+    H5Aclose(attr);
+    // The int can be casted to a bool
+    bool bool_value = (bool)int_value;
+    return bool_value;
+  }
+
+  template <>
+  std::string
+  HDF5Object::attr(const std::string attr_name) const
+  {
+    // Reads a UTF8 python string
+    //
+    // code from
+    // https://support.hdfgroup.org/ftp/HDF5/examples/misc-examples/vlstratt.c
+    //
+    // In the case of a variable length string the user does not have to reserve
+    // memory for string_out. The call HAread will reserve the memory and the
+    // user has to free the memory.
+    //
+    // Todo:
+    // - In debug mode, check the type of the attribute. If it is not a variable
+    // length string raise an exception
+    // - Use valgrind to check that there are no memory leaks
+    // - Use H5Dvlen_reclaim instead of free
+    // - Use collective
+
+    char * string_out;
+    hid_t  attr;
+    hid_t  type;
+    herr_t ret;
+
+    /* Create a datatype to refer to. */
+    type = H5Tcopy(H5T_C_S1);
+    Assert(type >= 0, ExcInternalError());
+
+    // Python strings are encoded in UTF8
+    ret = H5Tset_cset(type, H5T_CSET_UTF8);
+    Assert(type >= 0, ExcInternalError());
+
+    ret = H5Tset_size(type, H5T_VARIABLE);
+    Assert(ret >= 0, ExcInternalError());
+
+    attr = H5Aopen(*hdf5_reference, attr_name.data(), H5P_DEFAULT);
+    Assert(attr >= 0, ExcInternalError());
+
+    ret = H5Aread(attr, type, &string_out);
+    Assert(ret >= 0, ExcInternalError());
+
+    std::string string_value(string_out);
+    // The memory of the variable length string has to be freed.
+    // H5Dvlen_reclaim could be also used
+    free(string_out);
+    H5Tclose(type);
+    H5Aclose(attr);
+    return string_value;
+  }
+
+  template <>
+  std::vector<std::vector<double>>
+  HDF5Object::attr(const std::string attr_name) const
+  {
+    hid_t attr;
+    hid_t attr_space;
+    attr       = H5Aopen(*hdf5_reference, attr_name.data(), H5P_DEFAULT);
+    attr_space = H5Aget_space(attr);
+    AssertDimension(H5Sget_simple_extent_ndims(attr_space), 2);
+    hsize_t dims[2];
+    H5Sget_simple_extent_dims(attr_space, dims, NULL);
+
+    std::vector<std::vector<double>> matrix_value(dims[0],
+                                                  std::vector<double>(dims[1]));
+
+    double *hdf5_data = (double *)malloc(dims[0] * dims[1] * sizeof(double));
+    H5Aread(attr, H5T_NATIVE_DOUBLE, hdf5_data);
+    for (unsigned int dim0_idx = 0; dim0_idx < dims[0]; ++dim0_idx)
+      {
+        for (unsigned int dim1_idx = 0; dim1_idx < dims[1]; ++dim1_idx)
+          {
+            matrix_value[dim0_idx][dim1_idx] =
+              hdf5_data[dim0_idx * dims[1] + dim1_idx];
+          }
+      }
+
+    H5Sclose(attr_space);
+    H5Aclose(attr);
+
+    free(hdf5_data);
+
+    return matrix_value;
+  }
+
+  template <typename T>
+  void
+  HDF5Object::write_attr(const std::string attr_name, T value) const
+  {
+    hid_t                  attr;
+    hid_t                  aid;
+    std::shared_ptr<hid_t> t_type;
+
+    if (std::is_same<T, double>::value)
+      {
+        t_type  = std::shared_ptr<hid_t>(new hid_t);
+        *t_type = H5T_NATIVE_DOUBLE;
+      }
+    else if (std::is_same<T, int>::value)
+      {
+        t_type  = std::shared_ptr<hid_t>(new hid_t);
+        *t_type = H5T_NATIVE_INT;
+      }
+    else if (std::is_same<T, unsigned int>::value)
+      {
+        t_type  = std::shared_ptr<hid_t>(new hid_t);
+        *t_type = H5T_NATIVE_UINT;
+      }
+    else if (std::is_same<T, std::complex<double>>::value)
+      {
+        t_type  = std::shared_ptr<hid_t>(new hid_t, [](auto pointer) {
+          // Relase the HDF5 resource
+          H5Tclose(*pointer);
+          delete pointer;
+        });
+        *t_type = H5Tcreate(H5T_COMPOUND, sizeof(std::complex<double>));
+        //  The C++ standards committee agreed to mandate that the storage
+        //  format used for the std::complex type be binary-compatible with the
+        //  C99 type, i.e. an array T[2] with consecutive real [0] and imaginary
+        //  [1] parts.
+        H5Tinsert(*t_type, "r", 0, H5T_NATIVE_DOUBLE);
+        H5Tinsert(*t_type, "i", sizeof(double), H5T_NATIVE_DOUBLE);
+      }
+    else
+      {
+        Assert(false, ExcInternalError());
+      }
+
+    /*
+     * Create scalar attribute.
+     */
+    aid  = H5Screate(H5S_SCALAR);
+    attr = H5Acreate2(*hdf5_reference,
+                      attr_name.data(),
+                      *t_type,
+                      aid,
+                      H5P_DEFAULT,
+                      H5P_DEFAULT);
+
+    /*
+     * Write scalar attribute.
+     */
+    H5Awrite(attr, *t_type, &value);
+
+    H5Sclose(aid);
+    H5Aclose(attr);
+  }
+
+  template <typename T>
+  DataSet<T>::DataSet(const std::string    name,
+                      const hid_t &        parent_group_id,
+                      std::vector<hsize_t> dimensions,
+                      const bool           mpi,
+                      const Mode           mode)
+    : HDF5Object(name, mpi)
+    , rank(dimensions.size())
+    , dimensions(dimensions)
+  {
+    hdf5_reference = std::shared_ptr<hid_t>(new hid_t, [](auto pointer) {
+      // Relase the HDF5 resource
+      H5Dclose(*pointer);
+      delete pointer;
+    });
+    dataspace      = std::shared_ptr<hid_t>(new hid_t, [](auto pointer) {
+      // Relase the HDF5 resource
+      H5Sclose(*pointer);
+      delete pointer;
+    });
+
+    if (std::is_same<T, double>::value)
+      {
+        t_type  = std::shared_ptr<hid_t>(new hid_t);
+        *t_type = H5T_NATIVE_DOUBLE;
+      }
+    else if (std::is_same<T, int>::value)
+      {
+        t_type  = std::shared_ptr<hid_t>(new hid_t);
+        *t_type = H5T_NATIVE_INT;
+      }
+    else if (std::is_same<T, unsigned int>::value)
+      {
+        t_type  = std::shared_ptr<hid_t>(new hid_t);
+        *t_type = H5T_NATIVE_UINT;
+      }
+    else if (std::is_same<T, std::complex<double>>::value)
+      {
+        t_type  = std::shared_ptr<hid_t>(new hid_t, [](auto pointer) {
+          // Relase the HDF5 resource
+          H5Tclose(*pointer);
+          delete pointer;
+        });
+        *t_type = H5Tcreate(H5T_COMPOUND, sizeof(std::complex<double>));
+        //  The C++ standards committee agreed to mandate that the storage
+        //  format used for the std::complex type be binary-compatible with the
+        //  C99 type, i.e. an array T[2] with consecutive real [0] and imaginary
+        //  [1] parts.
+        H5Tinsert(*t_type, "r", 0, H5T_NATIVE_DOUBLE);
+        H5Tinsert(*t_type, "i", sizeof(double), H5T_NATIVE_DOUBLE);
+      }
+    else
+      {
+        Assert(false, ExcInternalError());
+      }
+
+    total_size = 1;
+    for (auto &&dimension : dimensions)
+      {
+        total_size *= dimension;
+      }
+
+    switch (mode)
+      {
+        case (Mode::create):
+          *dataspace = H5Screate_simple(rank, dimensions.data(), NULL);
+
+          *hdf5_reference = H5Dcreate2(parent_group_id,
+                                       name.data(),
+                                       *t_type,
+                                       *dataspace,
+                                       H5P_DEFAULT,
+                                       H5P_DEFAULT,
+                                       H5P_DEFAULT);
+          break;
+        case (Mode::open):
+          // Write the code for open
+          *hdf5_reference = H5Gopen2(parent_group_id, name.data(), H5P_DEFAULT);
+          break;
+        default:
+          Assert(false, ExcInternalError());
+          break;
+      }
+  }
+
+  template <typename T>
+  DataSet<T>::~DataSet()
+  {
+    // Close first the dataset
+    hdf5_reference.reset();
+    // Then close the dataspace
+    dataspace.reset();
+  }
+
+  template <typename T>
+  void
+  DataSet<T>::write_data(const std::vector<T> &data) const
+  {
+    AssertDimension(total_size, data.size());
+    hid_t plist;
+
+    if (mpi)
+      {
+        plist = H5Pcreate(H5P_DATASET_XFER);
+        H5Pset_dxpl_mpio(plist, H5FD_MPIO_COLLECTIVE);
+      }
+    else
+      {
+        plist = H5P_DEFAULT;
+      }
+
+    H5Dwrite(*hdf5_reference, *t_type, H5S_ALL, H5S_ALL, plist, data.data());
+
+    if (mpi)
+      {
+        H5Pclose(plist);
+      }
+  }
+
+  template <typename T>
+  void
+  DataSet<T>::write_data(const dealii::FullMatrix<T> &data) const
+  {
+    AssertDimension(total_size, data.m() * data.n());
+    hid_t plist;
+
+    if (mpi)
+      {
+        plist = H5Pcreate(H5P_DATASET_XFER);
+        H5Pset_dxpl_mpio(plist, H5FD_MPIO_COLLECTIVE);
+      }
+    else
+      {
+        plist = H5P_DEFAULT;
+      }
+
+    // The iterator of FullMatrix has to be converted to a pointer to the raw
+    // data
+    H5Dwrite(*hdf5_reference, *t_type, H5S_ALL, H5S_ALL, plist, &*data.begin());
+
+    if (mpi)
+      {
+        H5Pclose(plist);
+      }
+  }
+
+  template <typename T>
+  void
+  DataSet<T>::write_data_selection(const std::vector<T> &     data,
+                                   const std::vector<hsize_t> coordinates) const
+  {
+    AssertDimension(coordinates.size(), data.size() * rank);
+    hid_t                memory_dataspace;
+    hid_t                plist;
+    std::vector<hsize_t> data_dimensions = {data.size()};
+
+    memory_dataspace = H5Screate_simple(1, data_dimensions.data(), NULL);
+    H5Sselect_elements(*dataspace,
+                       H5S_SELECT_SET,
+                       data.size(),
+                       coordinates.data());
+
+    if (mpi)
+      {
+        plist = H5Pcreate(H5P_DATASET_XFER);
+        H5Pset_dxpl_mpio(plist, H5FD_MPIO_COLLECTIVE);
+      }
+    else
+      {
+        plist = H5P_DEFAULT;
+      }
+
+    H5Dwrite(*hdf5_reference,
+             *t_type,
+             memory_dataspace,
+             *dataspace,
+             plist,
+             data.data());
+
+    if (mpi)
+      {
+        H5Pclose(plist);
+      }
+    H5Sclose(memory_dataspace);
+  }
+
+  template <typename T>
+  void
+  DataSet<T>::write_data_hyperslab(const std::vector<T> &     data,
+                                   const std::vector<hsize_t> offset,
+                                   const std::vector<hsize_t> count) const
+  {
+    AssertDimension(std::accumulate(count.begin(),
+                                    count.end(),
+                                    1,
+                                    std::multiplies<unsigned int>()),
+                    data.size());
+    hid_t                memory_dataspace;
+    hid_t                plist;
+    std::vector<hsize_t> data_dimensions = {data.size()};
+
+    memory_dataspace = H5Screate_simple(1, data_dimensions.data(), NULL);
+    H5Sselect_hyperslab(
+      *dataspace, H5S_SELECT_SET, offset.data(), NULL, count.data(), NULL);
+
+    if (mpi)
+      {
+        plist = H5Pcreate(H5P_DATASET_XFER);
+        H5Pset_dxpl_mpio(plist, H5FD_MPIO_COLLECTIVE);
+      }
+    else
+      {
+        plist = H5P_DEFAULT;
+      }
+    H5Dwrite(*hdf5_reference,
+             *t_type,
+             memory_dataspace,
+             *dataspace,
+             plist,
+             &*data.begin());
+
+    if (mpi)
+      {
+        H5Pclose(plist);
+      }
+    H5Sclose(memory_dataspace);
+  }
+
+  template <typename T>
+  void
+  DataSet<T>::write_data_hyperslab(const dealii::FullMatrix<T> &data,
+                                   const std::vector<hsize_t>   offset,
+                                   const std::vector<hsize_t>   count) const
+  {
+    AssertDimension(std::accumulate(count.begin(),
+                                    count.end(),
+                                    1,
+                                    std::multiplies<unsigned int>()),
+                    data.m() * data.n());
+    hid_t                memory_dataspace;
+    hid_t                plist;
+    std::vector<hsize_t> data_dimensions = {data.m(), data.n()};
+
+    memory_dataspace = H5Screate_simple(2, data_dimensions.data(), NULL);
+    H5Sselect_hyperslab(
+      *dataspace, H5S_SELECT_SET, offset.data(), NULL, count.data(), NULL);
+
+    if (mpi)
+      {
+        plist = H5Pcreate(H5P_DATASET_XFER);
+        H5Pset_dxpl_mpio(plist, H5FD_MPIO_COLLECTIVE);
+      }
+    else
+      {
+        plist = H5P_DEFAULT;
+      }
+    H5Dwrite(*hdf5_reference,
+             *t_type,
+             memory_dataspace,
+             *dataspace,
+             plist,
+             &*data.begin());
+
+    if (mpi)
+      {
+        H5Pclose(plist);
+      }
+    H5Sclose(memory_dataspace);
+  }
+
+  template <typename T>
+  void
+  DataSet<T>::write_data_none() const
+  {
+    hid_t                memory_dataspace;
+    hid_t                plist;
+    std::vector<hsize_t> data_dimensions = {0};
+
+    memory_dataspace = H5Screate_simple(1, data_dimensions.data(), NULL);
+    H5Sselect_none(*dataspace);
+
+    if (mpi)
+      {
+        plist = H5Pcreate(H5P_DATASET_XFER);
+        H5Pset_dxpl_mpio(plist, H5FD_MPIO_COLLECTIVE);
+      }
+    else
+      {
+        plist = H5P_DEFAULT;
+      }
+
+    // The pointer of data can safely be NULL, see the discussion at the HDF5
+    // forum:
+    // https://forum.hdfgroup.org/t/parallel-i-o-does-not-support-filters-yet/884/17
+    H5Dwrite(
+      *hdf5_reference, *t_type, memory_dataspace, *dataspace, plist, NULL);
+
+    if (mpi)
+      {
+        H5Pclose(plist);
+      }
+    H5Sclose(memory_dataspace);
+  }
+
+  Group::Group(const std::string name,
+               const Group &     parentGroup,
+               const bool        mpi,
+               const Mode        mode)
+    : HDF5Object(name, mpi)
+  {
+    hdf5_reference = std::shared_ptr<hid_t>(new hid_t, [](auto pointer) {
+      // Relase the HDF5 resource
+      H5Gclose(*pointer);
+      delete pointer;
+    });
+    switch (mode)
+      {
+        case (Mode::create):
+          *hdf5_reference = H5Gcreate2(*(parentGroup.hdf5_reference),
+                                       name.data(),
+                                       H5P_DEFAULT,
+                                       H5P_DEFAULT,
+                                       H5P_DEFAULT);
+          break;
+        case (Mode::open):
+          *hdf5_reference =
+            H5Gopen2(*(parentGroup.hdf5_reference), name.data(), H5P_DEFAULT);
+          break;
+        default:
+          Assert(false, ExcInternalError());
+          break;
+      }
+  }
+
+  Group::Group(const std::string name, const bool mpi)
+    : HDF5Object(name, mpi)
+  {}
+
+  Group
+  Group::group(const std::string name)
+  {
+    return Group(name, *this, mpi, Mode::open);
+  }
+
+  Group
+  Group::create_group(const std::string name)
+  {
+    return Group(name, *this, mpi, Mode::create);
+  }
+
+  template <>
+  void
+  Group::write_dataset(const std::string          name,
+                       const std::vector<double> &data) const
+  {
+    std::vector<hsize_t> dimensions = {data.size()};
+    auto                 dataset    = create_dataset<double>(name, dimensions);
+    dataset.write_data(data);
+  }
+
+  template <>
+  void
+  Group::write_dataset(const std::string                 name,
+                       const dealii::FullMatrix<double> &data) const
+  {
+    std::vector<hsize_t> dimensions = {data.m(), data.n()};
+    auto                 dataset    = create_dataset<double>(name, dimensions);
+    dataset.write_data(data);
+  }
+
+  template <typename T>
+  DataSet<T>
+  Group::create_dataset(const std::string          name,
+                        const std::vector<hsize_t> dimensions) const
+  {
+    return DataSet<T>(
+      name, *hdf5_reference, dimensions, mpi, DataSet<T>::Mode::create);
+  }
+
+  File::File(const std::string name,
+             const bool        mpi,
+             const MPI_Comm    mpi_communicator,
+             const Mode        mode)
+    : Group(name, mpi)
+  {
+    hdf5_reference = std::shared_ptr<hid_t>(new hid_t, [](auto pointer) {
+      // Relase the HDF5 resource
+      H5Fclose(*pointer);
+      delete pointer;
+    });
+
+    hid_t          plist;
+    const MPI_Info info = MPI_INFO_NULL;
+
+    if (mpi)
+      {
+        plist = H5Pcreate(H5P_FILE_ACCESS);
+        H5Pset_fapl_mpio(plist, mpi_communicator, info);
+      }
+    else
+      {
+        plist = H5P_DEFAULT;
+      }
+
+    switch (mode)
+      {
+        case (Mode::create):
+          *hdf5_reference =
+            H5Fcreate(name.data(), H5F_ACC_TRUNC, H5P_DEFAULT, plist);
+          break;
+        case (Mode::open):
+          *hdf5_reference = H5Fopen(name.data(), H5F_ACC_RDWR, plist);
+          break;
+        default:
+          Assert(false, ExcInternalError());
+          break;
+      }
+
+    if (mpi)
+      {
+        // Relase the HDF5 resource
+        H5Pclose(plist);
+      }
+  }
+
+  File::File(const std::string name,
+             const MPI_Comm    mpi_communicator,
+             const Mode        mode)
+    : File(name, true, mpi_communicator, mode)
+  {}
+
+  File::File(const std::string name, const Mode mode)
+    : File(name, false, MPI_COMM_NULL, mode)
+  {}
+
+
+  // explicit instantiations of classes
+  template class DataSet<double>;
+  template class DataSet<std::complex<double>>;
+
+  // explicit instantiations of functions
+  template double
+  HDF5Object::attr<double>(const std::string attr_name) const;
+  template int
+  HDF5Object::attr<int>(const std::string attr_name) const;
+  template unsigned int
+  HDF5Object::attr<unsigned int>(const std::string attr_name) const;
+  template std::complex<double>
+  HDF5Object::attr<std::complex<double>>(const std::string attr_name) const;
+  // The specialization of HDF5Object::attr<std::string> has been defined above
+
+  template void
+  HDF5Object::write_attr<double>(const std::string attr_name,
+                                 double            value) const;
+  template void
+  HDF5Object::write_attr<int>(const std::string attr_name, int value) const;
+  template void
+  HDF5Object::write_attr<unsigned int>(const std::string attr_name,
+                                       unsigned int      value) const;
+  template void
+  HDF5Object::write_attr<std::complex<double>>(
+    const std::string    attr_name,
+    std::complex<double> value) const;
+
+  template DataSet<double>
+  Group::create_dataset<double>(const std::string          name,
+                                const std::vector<hsize_t> dimensions) const;
+  template DataSet<int>
+  Group::create_dataset<int>(const std::string          name,
+                             const std::vector<hsize_t> dimensions) const;
+  template DataSet<unsigned int>
+  Group::create_dataset<unsigned int>(
+    const std::string          name,
+    const std::vector<hsize_t> dimensions) const;
+  template DataSet<std::complex<double>>
+  Group::create_dataset<std::complex<double>>(
+    const std::string          name,
+    const std::vector<hsize_t> dimensions) const;
+
+
+} // namespace hdf5
+
+DEAL_II_NAMESPACE_CLOSE

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.