DOLFINx 0.10.0.0
DOLFINx C++ interface
Loading...
Searching...
No Matches
xdmf_utils.h
1// Copyright (C) 2012-2024 Chris N. Richardson, Jørgen S. Dokken
2//
3// This file is part of DOLFINx (https://www.fenicsproject.org)
4//
5// SPDX-License-Identifier: LGPL-3.0-or-later
6
7#pragma once
8
9#include "HDF5Interface.h"
10#include <array>
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>
17#include <filesystem>
18#include <numeric>
19#include <pugixml.hpp>
20#include <span>
21#include <string>
22#include <utility>
23#include <vector>
24
25namespace dolfinx
26{
27namespace fem
28{
29template <dolfinx::scalar T, std::floating_point U>
30class Function;
31} // namespace fem
32
33namespace fem
34{
35template <std::floating_point T>
38} // namespace fem
39
40namespace mesh
41{
42template <std::floating_point T>
43class Mesh;
44class Topology;
45} // namespace mesh
46
47namespace io::xdmf_utils
48{
49
52std::pair<std::string, int> get_cell_type(const pugi::xml_node& topology_node);
53
56std::array<std::string, 2> get_hdf5_paths(const pugi::xml_node& dataitem_node);
57
58std::filesystem::path
59get_hdf5_filename(const std::filesystem::path& xdmf_filename);
60
62std::vector<std::int64_t> get_dataset_shape(const pugi::xml_node& dataset_node);
63
65std::int64_t get_num_cells(const pugi::xml_node& topology_node);
66
68std::string vtk_cell_type_str(mesh::CellType cell_type, int num_nodes);
69
71template <typename T>
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)
76{
77 // Add DataItem node
78 assert(xml_node);
79 pugi::xml_node data_item_node = xml_node.append_child("DataItem");
80 assert(data_item_node);
81
82 // Add dimensions attribute
83 std::string dims;
84 for (auto d : shape)
85 dims += std::to_string(d) + std::string(" ");
86 dims.pop_back();
87 data_item_node.append_attribute("Dimensions") = dims.c_str();
88
89 // Set type for topology data (needed by XDMF to prevent default to
90 // float)
91 if (!number_type.empty())
92 data_item_node.append_attribute("NumberType") = number_type.c_str();
93
94 // Add format attribute
95 if (h5_id < 0)
96 {
97 data_item_node.append_attribute("Format") = "XML";
98 assert(shape.size() == 2);
99 std::ostringstream s;
100 s.precision(16);
101 for (std::size_t i = 0; i < x.size(); ++i)
102 {
103 if ((i + 1) % shape[1] == 0 and shape[1] != 0)
104 s << x.data()[i] << std::endl;
105 else
106 s << x.data()[i] << " ";
107 }
108
109 data_item_node.append_child(pugi::node_pcdata).set_value(s.str().c_str());
110 }
111 else
112 {
113 data_item_node.append_attribute("Format") = "HDF";
114
115 // Get name of HDF5 file, including path
116 const std::filesystem::path p = io::hdf5::get_filename(h5_id);
117 const std::filesystem::path filename = p.filename().c_str();
118
119 // Add HDF5 filename and HDF5 internal path to XML file
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());
123
124 // Compute data offset and range of values
125 std::int64_t local_shape0 = std::reduce(
126 std::next(shape.begin()), shape.end(), x.size(), std::divides{});
127
128 const std::array local_range{offset, offset + local_shape0};
129 io::hdf5::write_dataset(h5_id, h5_path, x.data(), local_range, shape,
130 use_mpi_io, false);
131
132 // Add partitioning attribute to dataset
133 // std::vector<std::size_t> partitions;
134 // std::vector<std::size_t> offset_tmp(1, offset);
135 // dolfinx::MPI::gather(comm, offset_tmp, partitions);
136 // dolfinx::MPI::broadcast(comm, partitions);
137 // io::hdf5::add_attribute(h5_id, h5_path, "partition", partitions);
138 }
139}
140
145template <typename T>
146std::vector<T> get_dataset(MPI_Comm comm, const pugi::xml_node& dataset_node,
147 hid_t h5_id,
148 std::array<std::int64_t, 2> range = {0, 0})
149{
150 // FIXME: Need to sort out dataset dimensions - can't depend on HDF5
151 // shape, and a Topology data item is not required to have a
152 // 'Dimensions' attribute since the dimensions can be determined from
153 // the number of cells and the cell type (for topology, one must
154 // supply cell type + (number of cells or dimensions)).
155 //
156 // A geometry data item must have 'Dimensions' attribute.
157
158 assert(dataset_node);
159 pugi::xml_attribute format_attr = dataset_node.attribute("Format");
160 assert(format_attr);
161
162 // Get data set shape from 'Dimensions' attribute (empty if not
163 // available)
164 const std::vector shape_xml = xdmf_utils::get_dataset_shape(dataset_node);
165
166 const std::string format = format_attr.as_string();
167 std::vector<T> data_vector;
168 // Only read ASCII on process 0
169 const int mpi_rank = dolfinx::MPI::rank(comm);
170 if (format == "XML")
171 {
172 if (mpi_rank == 0)
173 {
174 // Read data and trim any leading/trailing whitespace
175 pugi::xml_node data_node = dataset_node.first_child();
176 assert(data_node);
177 std::string data_str = data_node.value();
178
179 // Split data based on spaces and line breaks
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"));
182
183 // Add data to numerical vector
184 data_vector.reserve(data_vector_str.size());
185 for (auto& v : data_vector_str)
186 {
187 if (v.begin() != v.end())
188 data_vector.push_back(
189 boost::lexical_cast<T>(boost::copy_range<std::string>(v)));
190 }
191 }
192 }
193 else if (format == "HDF")
194 {
195 // Get file and data path
196 auto paths = xdmf_utils::get_hdf5_paths(dataset_node);
197
198 // Get data shape from HDF5 file
199 const std::vector shape_hdf5 = io::hdf5::get_dataset_shape(h5_id, paths[1]);
200
201 // FIXME: should we support empty data sets?
202 // Check that data set is not empty
203 assert(!shape_hdf5.empty());
204 assert(shape_hdf5[0] != 0);
205
206 // Determine range of data to read from HDF5 file. This is
207 // complicated by the XML Dimension attribute and the HDF5 storage
208 // possibly having different shapes, e.g. the HDF5 storage may be a
209 // flat array.
210
211 // If range = {0, 0} then no range is supplied and we must determine
212 // the range
213 if (range[0] == 0 and range[1] == 0)
214 {
215 if (shape_xml == shape_hdf5)
216 {
217 range = dolfinx::MPI::local_range(mpi_rank, shape_hdf5[0],
218 dolfinx::MPI::size(comm));
219 }
220 else if (!shape_xml.empty() and shape_hdf5.size() == 1)
221 {
222 // Size of dims > 0
223 std::int64_t d = std::reduce(shape_xml.begin(), shape_xml.end(),
224 std::int64_t(1), std::multiplies{});
225
226 // Check for data size consistency
227 if (d * shape_xml[0] != shape_hdf5[0])
228 {
229 throw std::runtime_error("Data size in XDMF/XML and size of HDF5 "
230 "dataset are inconsistent");
231 }
232
233 // Compute data range to read
234 range = dolfinx::MPI::local_range(mpi_rank, shape_xml[0],
235 dolfinx::MPI::rank(comm));
236 range[0] *= d;
237 range[1] *= d;
238 }
239 else
240 {
241 throw std::runtime_error("This combination of array shapes in XDMF and "
242 "HDF5 is not supported");
243 }
244 }
245
246 // Retrieve data
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.");
250 else
251 {
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.");
255 }
256 }
257 else
258 throw std::runtime_error("Storage format \"" + format + "\" is unknown");
259
260 // Get dimensions for consistency (if available in DataItem node)
261 if (shape_xml.empty())
262 {
263 std::int64_t size = 1;
264 for (auto dim : shape_xml)
265 size *= dim;
266
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)
271 {
272 throw std::runtime_error(
273 "Data sizes in attribute and size of data read are inconsistent");
274 }
275 }
276
277 return data_vector;
278}
279
280} // namespace io::xdmf_utils
281} // namespace dolfinx
Definition CoordinateElement.h:38
Definition ElementDofLayout.h:30
Definition Function.h:47
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