From: Pasquale Africa Date: Thu, 27 Aug 2020 20:13:17 +0000 (+0000) Subject: Support to HDF5 export of data on simplex meshes X-Git-Tag: v9.3.0-rc1~1146^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F10852%2Fhead;p=dealii.git Support to HDF5 export of data on simplex meshes --- diff --git a/doc/news/changes/minor/20200827PasqualeAfrica b/doc/news/changes/minor/20200827PasqualeAfrica new file mode 100644 index 0000000000..064ba9f46c --- /dev/null +++ b/doc/news/changes/minor/20200827PasqualeAfrica @@ -0,0 +1,3 @@ +Added: DataOut now supports HDF5 file format with simplex meshes. +
+(Pasquale Claudio Africa, 2020/08/27) diff --git a/include/deal.II/base/data_out_base.h b/include/deal.II/base/data_out_base.h index 03fcfeff17..db588abba1 100644 --- a/include/deal.II/base/data_out_base.h +++ b/include/deal.II/base/data_out_base.h @@ -1496,12 +1496,9 @@ namespace DataOutBase unsigned int node_dim; /** - * The number of vertices per cell. Equal to - * GeometryInfo::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. diff --git a/source/base/data_out_base.cc b/source/base/data_out_base.cc index 79faad781e..86fbe61c6c 100644 --- a/source/base/data_out_base.cc +++ b/source/base/data_out_base.cc @@ -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::vertices_per_cell; - vertices_per_cell = GeometryInfo::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(patch.n_subdivisions + 1); @@ -1493,6 +1502,7 @@ namespace unsigned int d3) { stream << GeometryInfo::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(patch.n_subdivisions + 1); @@ -4973,6 +4977,12 @@ namespace DataOutBase (1 + GeometryInfo::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::write_xdmf_file( << " \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 << " \n"; xdmf_file << " \n"; @@ -7544,12 +7559,12 @@ DataOutBase::write_hdf5_parallel( template void DataOutBase::write_hdf5_parallel( - const std::vector> & /*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> &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::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::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) << "\n"; + { + Assert(reference_cell_type == ReferenceCell::Type::Quad || + reference_cell_type == ReferenceCell::Type::Tri, + ExcNotImplemented()); + + ss << indent(indent_level + 1) << "\n" + << indent(indent_level + 2) << "\n" + << indent(indent_level + 2) << "\n"; + { + Assert(reference_cell_type == ReferenceCell::Type::Hex || + reference_cell_type == ReferenceCell::Type::Tet, + ExcNotImplemented()); + + ss << indent(indent_level + 1) << "\n" + << indent(indent_level + 2) << "\n" + << indent(indent_level + 2) << "\n"; + ss << "\" NumberType=\"UInt\" Format=\"HDF\">\n"; ss << indent(indent_level + 3) << h5_mesh_filename << ":/cells\n"; ss << indent(indent_level + 2) << "\n"; ss << indent(indent_level + 1) << "\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 index 0000000000..48f03976fb --- /dev/null +++ b/tests/simplex/data_out_write_hdf5_01.cc @@ -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 + +#include +#include + +#include +#include +#include + +#include +#include + +#include +#include + +#include "../tests.h" + +using namespace dealii; + +template +class RightHandSideFunction : public Function +{ +public: + RightHandSideFunction(const unsigned int n_components) + : Function(n_components) + {} + + virtual double + value(const Point &p, const unsigned int component = 0) const + { + return p[component % dim] * p[component % dim]; + } +}; + +template +void +test(const FiniteElement &fe, const unsigned int n_components) +{ + Triangulation tria; + GridGenerator::subdivided_hyper_cube_with_simplices(tria, dim == 2 ? 4 : 2); + + DoFHandler dof_handler(tria); + + dof_handler.distribute_dofs(fe); + + Vector solution(dof_handler.n_dofs()); + + MappingFE mapping(Simplex::FE_P(1)); + + VectorTools::interpolate(mapping, + dof_handler, + RightHandSideFunction(n_components), + solution); + + static unsigned int counter = 0; + + for (unsigned int n_subdivisions = 1; n_subdivisions <= 2; ++n_subdivisions) + { + DataOut 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 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(Simplex::FE_P(2) /*=degree*/, 1); + test(FESystem(Simplex::FE_P(2 /*=degree*/), dim), dim); + test(FESystem(Simplex::FE_P(2 /*=degree*/), + dim, + Simplex::FE_P(1 /*=degree*/), + 1), + dim + 1); + } + { + const unsigned int dim = 3; + test(Simplex::FE_P(2) /*=degree*/, 1); + test(FESystem(Simplex::FE_P(2 /*=degree*/), dim), dim); + test(FESystem(Simplex::FE_P(2 /*=degree*/), + dim, + Simplex::FE_P(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 index 0000000000..553a7f5878 --- /dev/null +++ b/tests/simplex/data_out_write_hdf5_01.mpirun=4.with_simplex_support=on.with_hdf5=on.output @@ -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