9#include "HDF5Interface.h"
11#include <basix/mdspan.hpp>
12#include <boost/algorithm/string.hpp>
13#include <boost/lexical_cast.hpp>
14#include <dolfinx/common/MPI.h>
15#include <dolfinx/common/types.h>
16#include <dolfinx/mesh/cell_types.h>
29template <dolfinx::scalar T, std::
floating_po
int U>
35template <std::
floating_po
int T>
42template <std::
floating_po
int T>
47namespace io::xdmf_utils
52std::pair<std::string, int> get_cell_type(
const pugi::xml_node& topology_node);
56std::array<std::string, 2> get_hdf5_paths(
const pugi::xml_node& dataitem_node);
59get_hdf5_filename(
const std::filesystem::path& xdmf_filename);
62std::vector<std::int64_t> get_dataset_shape(
const pugi::xml_node& dataset_node);
65std::int64_t get_num_cells(
const pugi::xml_node& topology_node);
68std::string vtk_cell_type_str(
mesh::CellType cell_type,
int num_nodes);
72void add_data_item(pugi::xml_node& xml_node, hid_t h5_id,
73 const std::string& h5_path, std::span<const T> x,
74 std::int64_t offset,
const std::vector<std::int64_t>& shape,
75 const std::string& number_type,
bool use_mpi_io)
79 pugi::xml_node data_item_node = xml_node.append_child(
"DataItem");
80 assert(data_item_node);
85 dims += std::to_string(d) + std::string(
" ");
87 data_item_node.append_attribute(
"Dimensions") = dims.c_str();
91 if (!number_type.empty())
92 data_item_node.append_attribute(
"NumberType") = number_type.c_str();
97 data_item_node.append_attribute(
"Format") =
"XML";
98 assert(shape.size() == 2);
101 for (std::size_t i = 0; i < x.size(); ++i)
103 if ((i + 1) % shape[1] == 0 and shape[1] != 0)
104 s << x.data()[i] << std::endl;
106 s << x.data()[i] <<
" ";
109 data_item_node.append_child(pugi::node_pcdata).set_value(s.str().c_str());
113 data_item_node.append_attribute(
"Format") =
"HDF";
116 const std::filesystem::path p = io::hdf5::get_filename(h5_id);
117 const std::filesystem::path filename = p.filename().c_str();
120 const std::string xdmf_path
121 = filename.string() + std::string(
":") + h5_path;
122 data_item_node.append_child(pugi::node_pcdata).set_value(xdmf_path.c_str());
125 std::int64_t local_shape0 = std::reduce(
126 std::next(shape.begin()), shape.end(), x.size(), std::divides{});
128 const std::array
local_range{offset, offset + local_shape0};
129 io::hdf5::write_dataset(h5_id, h5_path, x.data(), local_range, shape,
146std::vector<T> get_dataset(MPI_Comm comm,
const pugi::xml_node& dataset_node,
148 std::array<std::int64_t, 2> range = {0, 0})
158 assert(dataset_node);
159 pugi::xml_attribute format_attr = dataset_node.attribute(
"Format");
164 const std::vector shape_xml = xdmf_utils::get_dataset_shape(dataset_node);
166 const std::string format = format_attr.as_string();
167 std::vector<T> data_vector;
175 pugi::xml_node data_node = dataset_node.first_child();
177 std::string data_str = data_node.value();
180 std::vector<boost::iterator_range<std::string::iterator>> data_vector_str;
181 boost::split(data_vector_str, data_str, boost::is_any_of(
" \n"));
184 data_vector.reserve(data_vector_str.size());
185 for (
auto& v : data_vector_str)
187 if (v.begin() != v.end())
188 data_vector.push_back(
189 boost::lexical_cast<T>(boost::copy_range<std::string>(v)));
193 else if (format ==
"HDF")
196 auto paths = xdmf_utils::get_hdf5_paths(dataset_node);
199 const std::vector shape_hdf5 = io::hdf5::get_dataset_shape(h5_id, paths[1]);
203 assert(!shape_hdf5.empty());
204 assert(shape_hdf5[0] != 0);
213 if (range[0] == 0 and range[1] == 0)
215 if (shape_xml == shape_hdf5)
220 else if (!shape_xml.empty() and shape_hdf5.size() == 1)
223 std::int64_t d = std::reduce(shape_xml.begin(), shape_xml.end(),
224 std::int64_t(1), std::multiplies{});
227 if (d * shape_xml[0] != shape_hdf5[0])
229 throw std::runtime_error(
"Data size in XDMF/XML and size of HDF5 "
230 "dataset are inconsistent");
241 throw std::runtime_error(
"This combination of array shapes in XDMF and "
242 "HDF5 is not supported");
247 if (hid_t dset_id = io::hdf5::open_dataset(h5_id, paths[1]);
248 dset_id == H5I_INVALID_HID)
249 throw std::runtime_error(
"Failed to open HDF5 global dataset.");
252 data_vector = io::hdf5::read_dataset<T>(dset_id, range,
true);
253 if (herr_t err = H5Dclose(dset_id); err < 0)
254 throw std::runtime_error(
"Failed to close HDF5 global dataset.");
258 throw std::runtime_error(
"Storage format \"" + format +
"\" is unknown");
261 if (shape_xml.empty())
263 std::int64_t
size = 1;
264 for (
auto dim : shape_xml)
267 std::int64_t size_global = 0;
268 const std::int64_t size_local = data_vector.size();
269 MPI_Allreduce(&size_local, &size_global, 1, MPI_INT64_T, MPI_SUM, comm);
270 if (size != size_global)
272 throw std::runtime_error(
273 "Data sizes in attribute and size of data read are inconsistent");
Definition CoordinateElement.h:38
Definition ElementDofLayout.h:30
A Mesh consists of a set of connected and numbered mesh topological entities, and geometry data.
Definition Mesh.h:23
Topology stores the topology of a mesh, consisting of mesh entities and connectivity (incidence relat...
Definition Topology.h:46
int size(MPI_Comm comm)
Definition MPI.cpp:72
int rank(MPI_Comm comm)
Return process rank for the communicator.
Definition MPI.cpp:64
constexpr std::array< std::int64_t, 2 > local_range(int rank, std::int64_t N, int size)
Return local range for the calling process, partitioning the global [0, N - 1] range across all ranks...
Definition MPI.h:90
Finite element method functionality.
Definition assemble_expression_impl.h:23
Mesh data structures and algorithms on meshes.
Definition DofMap.h:32
CellType
Cell type identifier.
Definition cell_types.h:22
Top-level namespace.
Definition defines.h:12