]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Support to HDF5 export of data on simplex meshes 10852/head
authorPasquale Africa <pasqualeclaudio.africa@polimi.it>
Thu, 27 Aug 2020 20:13:17 +0000 (20:13 +0000)
committerPasquale Africa <pasqualeclaudio.africa@polimi.it>
Tue, 1 Sep 2020 13:34:14 +0000 (13:34 +0000)
doc/news/changes/minor/20200827PasqualeAfrica [new file with mode: 0644]
include/deal.II/base/data_out_base.h
source/base/data_out_base.cc
tests/simplex/data_out_write_hdf5_01.cc [new file with mode: 0644]
tests/simplex/data_out_write_hdf5_01.mpirun=4.with_simplex_support=on.with_hdf5=on.output [new file with mode: 0644]

diff --git a/doc/news/changes/minor/20200827PasqualeAfrica b/doc/news/changes/minor/20200827PasqualeAfrica
new file mode 100644 (file)
index 0000000..064ba9f
--- /dev/null
@@ -0,0 +1,3 @@
+Added: DataOut now supports HDF5 file format with simplex meshes.
+<br>
+(Pasquale Claudio Africa, 2020/08/27)
index 03fcfeff1741e875605bc09c0bbf5c2cb12e7ce5..db588abba1e3881c5c67ee8e397002051bb5b0e8 100644 (file)
@@ -1496,12 +1496,9 @@ namespace DataOutBase
     unsigned int node_dim;
 
     /**
-     * The number of vertices per cell. Equal to
-     * GeometryInfo<node_dim>::vertices_per_cell. We need to store
-     * it as a run-time variable here because the dimension
-     * node_dim is also a run-time variable.
+     * The number of cells stored in @ref filtered_cells.
      */
-    unsigned int vertices_per_cell;
+    unsigned int num_cells;
 
     /**
      * Map of points to an internal index.
@@ -3331,10 +3328,21 @@ public:
   /**
    * Get the XDMF content associated with this entry.
    * If the entry is not valid, this returns an empty string.
+   *
+   * @deprecated Use @ref get_xdmf_content(const unsigned int, const ReferenceCell::Type &) instead.
    */
+  DEAL_II_DEPRECATED
   std::string
   get_xdmf_content(const unsigned int indent_level) const;
 
+  /**
+   * Get the XDMF content associated with this entry.
+   * If the entry is not valid, this returns an empty string.
+   */
+  std::string
+  get_xdmf_content(const unsigned int         indent_level,
+                   const ReferenceCell::Type &reference_cell_type) const;
+
 private:
   /**
    * Whether this entry is valid and contains data to be written.
index 79faad781ee917e3c51cf39742de1eb836177d87..86fbe61c6c52d95a30fa04bc7642e1d8f47bc8cc 100644 (file)
@@ -315,7 +315,7 @@ namespace DataOutBase
   DataOutFilter::DataOutFilter()
     : flags(false, true)
     , node_dim(numbers::invalid_unsigned_int)
-    , vertices_per_cell(numbers::invalid_unsigned_int)
+    , num_cells(numbers::invalid_unsigned_int)
   {}
 
 
@@ -323,7 +323,7 @@ namespace DataOutBase
   DataOutFilter::DataOutFilter(const DataOutBase::DataOutFilterFlags &flags)
     : flags(flags)
     , node_dim(numbers::invalid_unsigned_int)
-    , vertices_per_cell(numbers::invalid_unsigned_int)
+    , num_cells(numbers::invalid_unsigned_int)
   {}
 
 
@@ -363,6 +363,10 @@ namespace DataOutBase
                                    const unsigned int pt_index)
   {
     filtered_cells[cell_index] = filtered_points[pt_index];
+
+    // (Re)-initialize counter at any first call to this method.
+    if (cell_index == 0)
+      num_cells = 1;
   }
 
 
@@ -432,7 +436,7 @@ namespace DataOutBase
   unsigned int
   DataOutFilter::n_cells() const
   {
-    return filtered_cells.size() / vertices_per_cell;
+    return num_cells;
   }
 
 
@@ -465,9 +469,11 @@ namespace DataOutBase
                             const unsigned int d2,
                             const unsigned int d3)
   {
+    ++num_cells;
+
     const unsigned int base_entry =
       index * GeometryInfo<dim>::vertices_per_cell;
-    vertices_per_cell = GeometryInfo<dim>::vertices_per_cell;
+
     internal_add_cell(base_entry + 0, start);
     if (dim >= 1)
       {
@@ -494,11 +500,14 @@ namespace DataOutBase
                                    const unsigned int start,
                                    const unsigned int n_points)
   {
-    (void)index;
-    (void)start;
-    (void)n_points;
+    ++num_cells;
 
-    Assert(false, ExcNotImplemented());
+    const unsigned int base_entry = index * n_points;
+
+    for (unsigned int i = 0; i < n_points; ++i)
+      {
+        internal_add_cell(base_entry + i, start + i);
+      }
   }
 
 
@@ -802,7 +811,7 @@ namespace
     n_cells = 0;
     for (const auto &patch : patches)
       {
-        // The following formula doesn't work for non-tensor products
+        // The following formula doesn't hold for non-tensor products.
         if (patch.reference_cell_type == ReferenceCell::get_hypercube(dim))
           {
             n_nodes += Utilities::fixed_power<dim>(patch.n_subdivisions + 1);
@@ -1493,6 +1502,7 @@ namespace
                         unsigned int d3)
   {
     stream << GeometryInfo<dim>::vertices_per_cell << '\t' << start;
+
     if (dim >= 1)
       stream << '\t' << start + d1;
     {
@@ -2502,7 +2512,7 @@ namespace DataOutBase
           {
             const unsigned int n_subdivisions = patch.n_subdivisions;
             const unsigned int n              = n_subdivisions + 1;
-            // Length of loops in all dimensons
+            // Length of loops in all dimensions
             const unsigned int n1 = (dim > 0) ? n_subdivisions : 1;
             const unsigned int n2 = (dim > 1) ? n_subdivisions : 1;
             const unsigned int n3 = (dim > 2) ? n_subdivisions : 1;
@@ -3584,7 +3594,7 @@ namespace DataOutBase
     }
 
     // max. and min. height of solution
-    Assert(patches.size() > 0, ExcInternalError());
+    Assert(patches.size() > 0, ExcNoPatches());
     double hmin = patches[0].data(0, 0);
     double hmax = patches[0].data(0, 0);
 
@@ -4948,14 +4958,8 @@ namespace DataOutBase
 
     for (const auto &patch : patches)
       {
-        // special treatment of simplices since they are not subdivided
-        if (patch.reference_cell_type != ReferenceCell::get_hypercube(dim))
-          {
-            n_nodes += patch.data.n_cols();
-            n_cells += 1;
-            n_points_an_n_cell += patch.data.n_cols() + 1;
-          }
-        else
+        // The following formulas don't hold for non-tensor products.
+        if (patch.reference_cell_type == ReferenceCell::get_hypercube(dim))
           {
             n_nodes += Utilities::fixed_power<dim>(patch.n_subdivisions + 1);
 
@@ -4973,6 +4977,12 @@ namespace DataOutBase
                   (1 + GeometryInfo<dim>::vertices_per_cell);
               }
           }
+        else
+          {
+            n_nodes += patch.data.n_cols();
+            n_cells += 1;
+            n_points_an_n_cell += patch.data.n_cols() + 1;
+          }
       }
 
     // in gmv format the vertex coordinates and the data have an order that is a
@@ -7352,8 +7362,13 @@ DataOutInterface<dim, spacedim>::write_xdmf_file(
         << "    <Grid Name=\"CellTime\" GridType=\"Collection\" CollectionType=\"Temporal\">\n";
 
       // Write out all the entries indented
+      const auto &patches = get_patches();
+      Assert(patches.size() > 0, DataOutBase::ExcNoPatches());
+
       for (it = entries.begin(); it != entries.end(); ++it)
-        xdmf_file << it->get_xdmf_content(3);
+        {
+          xdmf_file << it->get_xdmf_content(3, patches[0].reference_cell_type);
+        }
 
       xdmf_file << "    </Grid>\n";
       xdmf_file << "  </Domain>\n";
@@ -7544,12 +7559,12 @@ DataOutBase::write_hdf5_parallel(
 template <int dim, int spacedim>
 void
 DataOutBase::write_hdf5_parallel(
-  const std::vector<Patch<dim, spacedim>> & /*patches*/,
-  const DataOutBase::DataOutFilter &data_filter,
-  const bool                        write_mesh_file,
-  const std::string &               mesh_filename,
-  const std::string &               solution_filename,
-  MPI_Comm                          comm)
+  const std::vector<Patch<dim, spacedim>> &patches,
+  const DataOutBase::DataOutFilter &       data_filter,
+  const bool                               write_mesh_file,
+  const std::string &                      mesh_filename,
+  const std::string &                      solution_filename,
+  MPI_Comm                                 comm)
 {
   AssertThrow(
     spacedim >= 2,
@@ -7562,6 +7577,7 @@ DataOutBase::write_hdf5_parallel(
 #ifndef DEAL_II_WITH_HDF5
   // throw an exception, but first make sure the compiler does not warn about
   // the now unused function arguments
+  (void)patches;
   (void)data_filter;
   (void)write_mesh_file;
   (void)mesh_filename;
@@ -7570,15 +7586,19 @@ DataOutBase::write_hdf5_parallel(
   AssertThrow(false, ExcMessage("HDF5 support is disabled."));
 #else
 #  ifndef DEAL_II_WITH_MPI
+  (void)comm;
+#  endif
+
   // verify that there are indeed patches to be written out. most of the times,
   // people just forget to call build_patches when there are no patches, so a
   // warning is in order. that said, the assertion is disabled if we support MPI
   // since then it can happen that on the coarsest mesh, a processor simply has
   // no cells it actually owns, and in that case it is legit if there are no
   // patches
-  Assert(data_filter.n_nodes() > 0, ExcNoPatches());
-  (void)comm;
-#  endif
+  Assert(patches.size() > 0, ExcNoPatches());
+
+  const auto &cell_info =
+    ReferenceCell::internal::Info::get_cell(patches[0].reference_cell_type);
 
   hid_t h5_mesh_file_id = -1, h5_solution_file_id, file_plist_id, plist_id;
   hid_t node_dataspace, node_dataset, node_file_dataspace,
@@ -7673,7 +7693,7 @@ DataOutBase::write_hdf5_parallel(
       AssertThrow(node_dataspace >= 0, ExcIO());
 
       cell_ds_dim[0] = global_node_cell_count[1];
-      cell_ds_dim[1] = GeometryInfo<dim>::vertices_per_cell;
+      cell_ds_dim[1] = cell_info.n_vertices();
       cell_dataspace = H5Screate_simple(2, cell_ds_dim, nullptr);
       AssertThrow(cell_dataspace >= 0, ExcIO());
 
@@ -7734,7 +7754,7 @@ DataOutBase::write_hdf5_parallel(
 
       // And repeat for cells
       count[0] = local_node_cell_count[1];
-      count[1] = GeometryInfo<dim>::vertices_per_cell;
+      count[1] = cell_info.n_vertices();
       offset[0] = global_node_cell_offsets[1];
       offset[1] = 0;
       cell_memory_dataspace = H5Screate_simple(2, count, nullptr);
@@ -8501,6 +8521,17 @@ namespace
 
 std::string
 XDMFEntry::get_xdmf_content(const unsigned int indent_level) const
+{
+  return get_xdmf_content(indent_level,
+                          ReferenceCell::get_hypercube(dimension));
+}
+
+
+
+std::string
+XDMFEntry::get_xdmf_content(
+  const unsigned int         indent_level,
+  const ReferenceCell::Type &reference_cell_type) const
 {
   if (!valid)
     return "";
@@ -8531,17 +8562,51 @@ XDMFEntry::get_xdmf_content(const unsigned int indent_level) const
            << "\" NumberOfElements=\"" << num_cells
            << "\" NodesPerElement=\"2\">\n";
       else if (dimension == 2)
-        ss << indent(indent_level + 1) << "<Topology TopologyType=\""
-           << "Quadrilateral"
-           << "\" NumberOfElements=\"" << num_cells << "\">\n";
+        {
+          Assert(reference_cell_type == ReferenceCell::Type::Quad ||
+                   reference_cell_type == ReferenceCell::Type::Tri,
+                 ExcNotImplemented());
+
+          ss << indent(indent_level + 1) << "<Topology TopologyType=\"";
+          if (reference_cell_type == ReferenceCell::Type::Quad)
+            {
+              ss << "Quadrilateral"
+                 << "\" NumberOfElements=\"" << num_cells << "\">\n"
+                 << indent(indent_level + 2) << "<DataItem Dimensions=\""
+                 << num_cells << " " << (1 << dimension);
+            }
+          else // if (reference_cell_type == ReferenceCell::Type::Tri)
+            {
+              ss << "Triangle"
+                 << "\" NumberOfElements=\"" << num_cells << "\">\n"
+                 << indent(indent_level + 2) << "<DataItem Dimensions=\""
+                 << num_cells << " " << 3;
+            }
+        }
       else if (dimension == 3)
-        ss << indent(indent_level + 1) << "<Topology TopologyType=\""
-           << "Hexahedron"
-           << "\" NumberOfElements=\"" << num_cells << "\">\n";
+        {
+          Assert(reference_cell_type == ReferenceCell::Type::Hex ||
+                   reference_cell_type == ReferenceCell::Type::Tet,
+                 ExcNotImplemented());
+
+          ss << indent(indent_level + 1) << "<Topology TopologyType=\"";
+          if (reference_cell_type == ReferenceCell::Type::Hex)
+            {
+              ss << "Hexahedron"
+                 << "\" NumberOfElements=\"" << num_cells << "\">\n"
+                 << indent(indent_level + 2) << "<DataItem Dimensions=\""
+                 << num_cells << " " << (1 << dimension);
+            }
+          else // if (reference_cell_type == ReferenceCell::Type::Tet)
+            {
+              ss << "Tetrahedron"
+                 << "\" NumberOfElements=\"" << num_cells << "\">\n"
+                 << indent(indent_level + 2) << "<DataItem Dimensions=\""
+                 << num_cells << " " << 4;
+            }
+        }
 
-      ss << indent(indent_level + 2) << "<DataItem Dimensions=\"" << num_cells
-         << " " << (1 << dimension)
-         << "\" NumberType=\"UInt\" Format=\"HDF\">\n";
+      ss << "\" NumberType=\"UInt\" Format=\"HDF\">\n";
       ss << indent(indent_level + 3) << h5_mesh_filename << ":/cells\n";
       ss << indent(indent_level + 2) << "</DataItem>\n";
       ss << indent(indent_level + 1) << "</Topology>\n";
diff --git a/tests/simplex/data_out_write_hdf5_01.cc b/tests/simplex/data_out_write_hdf5_01.cc
new file mode 100644 (file)
index 0000000..48f0397
--- /dev/null
@@ -0,0 +1,146 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// Test DataOut with HDF5 for simplex meshes.
+
+#include <deal.II/dofs/dof_handler.h>
+
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/mapping_fe.h>
+
+#include <deal.II/grid/grid_in.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <deal.II/simplex/fe_lib.h>
+#include <deal.II/simplex/grid_generator.h>
+
+#include "../tests.h"
+
+using namespace dealii;
+
+template <int dim>
+class RightHandSideFunction : public Function<dim>
+{
+public:
+  RightHandSideFunction(const unsigned int n_components)
+    : Function<dim>(n_components)
+  {}
+
+  virtual double
+  value(const Point<dim> &p, const unsigned int component = 0) const
+  {
+    return p[component % dim] * p[component % dim];
+  }
+};
+
+template <int dim, int spacedim = dim>
+void
+test(const FiniteElement<dim, spacedim> &fe, const unsigned int n_components)
+{
+  Triangulation<dim, spacedim> tria;
+  GridGenerator::subdivided_hyper_cube_with_simplices(tria, dim == 2 ? 4 : 2);
+
+  DoFHandler<dim> dof_handler(tria);
+
+  dof_handler.distribute_dofs(fe);
+
+  Vector<double> solution(dof_handler.n_dofs());
+
+  MappingFE<dim> mapping(Simplex::FE_P<dim>(1));
+
+  VectorTools::interpolate(mapping,
+                           dof_handler,
+                           RightHandSideFunction<dim>(n_components),
+                           solution);
+
+  static unsigned int counter = 0;
+
+  for (unsigned int n_subdivisions = 1; n_subdivisions <= 2; ++n_subdivisions)
+    {
+      DataOut<dim> data_out;
+
+      data_out.attach_dof_handler(dof_handler);
+      data_out.add_data_vector(solution, "solution");
+
+      data_out.build_patches(mapping, n_subdivisions);
+
+      const std::string output_basename("test." + std::to_string(dim) + "." +
+                                        std::to_string(counter++));
+
+      DataOutBase::DataOutFilter data_filter(
+        DataOutBase::DataOutFilterFlags(true, true));
+      data_out.write_filtered_data(data_filter);
+      data_out.write_hdf5_parallel(data_filter,
+                                   output_basename + ".h5",
+                                   MPI_COMM_SELF);
+
+      std::vector<XDMFEntry> xdmf_entries({data_out.create_xdmf_entry(
+        data_filter, output_basename + ".h5", 0, MPI_COMM_SELF)});
+
+      data_out.write_xdmf_file(xdmf_entries,
+                               output_basename + ".xdmf",
+                               MPI_COMM_SELF);
+
+      data_out.clear();
+
+      // Sadly hdf5 is binary and we can not use hd5dump because it might
+      // not be in the path. At least we can make sure that both the xdmf and
+      // the h5 file are created.
+      std::ifstream h5((output_basename + ".h5").c_str());
+      AssertThrow(h5.good(), ExcIO());
+
+      std::ifstream xdmf((output_basename + ".xdmf").c_str());
+      AssertThrow(h5.good(), ExcIO());
+
+      deallog << "Files " << output_basename + ".h5"
+              << " and " << output_basename + ".xdmf"
+              << " created succesfully!" << std::endl;
+    }
+}
+
+int
+main(int argc, char **argv)
+{
+  initlog();
+
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+
+  {
+    const unsigned int dim = 2;
+    test<dim>(Simplex::FE_P<dim>(2) /*=degree*/, 1);
+    test<dim>(FESystem<dim>(Simplex::FE_P<dim>(2 /*=degree*/), dim), dim);
+    test<dim>(FESystem<dim>(Simplex::FE_P<dim>(2 /*=degree*/),
+                            dim,
+                            Simplex::FE_P<dim>(1 /*=degree*/),
+                            1),
+              dim + 1);
+  }
+  {
+    const unsigned int dim = 3;
+    test<dim>(Simplex::FE_P<dim>(2) /*=degree*/, 1);
+    test<dim>(FESystem<dim>(Simplex::FE_P<dim>(2 /*=degree*/), dim), dim);
+    test<dim>(FESystem<dim>(Simplex::FE_P<dim>(2 /*=degree*/),
+                            dim,
+                            Simplex::FE_P<dim>(1 /*=degree*/),
+                            1),
+              dim + 1);
+  }
+}
diff --git a/tests/simplex/data_out_write_hdf5_01.mpirun=4.with_simplex_support=on.with_hdf5=on.output b/tests/simplex/data_out_write_hdf5_01.mpirun=4.with_simplex_support=on.with_hdf5=on.output
new file mode 100644 (file)
index 0000000..553a7f5
--- /dev/null
@@ -0,0 +1,13 @@
+
+DEAL::Files test.2.0.h5 and test.2.0.xdmf created succesfully!
+DEAL::Files test.2.1.h5 and test.2.1.xdmf created succesfully!
+DEAL::Files test.2.2.h5 and test.2.2.xdmf created succesfully!
+DEAL::Files test.2.3.h5 and test.2.3.xdmf created succesfully!
+DEAL::Files test.2.4.h5 and test.2.4.xdmf created succesfully!
+DEAL::Files test.2.5.h5 and test.2.5.xdmf created succesfully!
+DEAL::Files test.3.0.h5 and test.3.0.xdmf created succesfully!
+DEAL::Files test.3.1.h5 and test.3.1.xdmf created succesfully!
+DEAL::Files test.3.2.h5 and test.3.2.xdmf created succesfully!
+DEAL::Files test.3.3.h5 and test.3.3.xdmf created succesfully!
+DEAL::Files test.3.4.h5 and test.3.4.xdmf created succesfully!
+DEAL::Files test.3.5.h5 and test.3.5.xdmf created succesfully!
\ No newline at end of file

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.