DOLFINx 0.9.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 <complex>
15#include <concepts>
16#include <dolfinx/common/MPI.h>
17#include <dolfinx/common/types.h>
18#include <dolfinx/mesh/cell_types.h>
19#include <filesystem>
20#include <numeric>
21#include <pugixml.hpp>
22#include <span>
23#include <string>
24#include <utility>
25#include <vector>
26
27namespace dolfinx
28{
29
30namespace fem
31{
32template <dolfinx::scalar T, std::floating_point U>
33class Function;
34} // namespace fem
35
36namespace fem
37{
38template <std::floating_point T>
39class CoordinateElement;
40class ElementDofLayout;
41} // namespace fem
42
43namespace mesh
44{
45template <std::floating_point T>
46class Mesh;
47class Topology;
48} // namespace mesh
49
50namespace io::xdmf_utils
51{
52
55std::pair<std::string, int> get_cell_type(const pugi::xml_node& topology_node);
56
59std::array<std::string, 2> get_hdf5_paths(const pugi::xml_node& dataitem_node);
60
61std::filesystem::path
62get_hdf5_filename(const std::filesystem::path& xdmf_filename);
63
65std::vector<std::int64_t> get_dataset_shape(const pugi::xml_node& dataset_node);
66
68std::int64_t get_num_cells(const pugi::xml_node& topology_node);
69
71std::string vtk_cell_type_str(mesh::CellType cell_type, int num_nodes);
72
107template <typename T>
108std::pair<std::vector<std::int32_t>, std::vector<T>> distribute_entity_data(
109 const mesh::Topology& topology, std::span<const std::int64_t> nodes_g,
110 std::int64_t num_nodes_g, const fem::ElementDofLayout& cmap_dof_layout,
111 MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
112 const std::int32_t,
113 MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 2>>
114 xdofmap,
115 int entity_dim,
116 MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
117 const std::int64_t,
118 MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 2>>
119 entities,
120 std::span<const T> data);
121
123template <typename T>
124void add_data_item(pugi::xml_node& xml_node, hid_t h5_id,
125 const std::string& h5_path, std::span<const T> x,
126 std::int64_t offset, const std::vector<std::int64_t>& shape,
127 const std::string& number_type, bool use_mpi_io)
128{
129 // Add DataItem node
130 assert(xml_node);
131 pugi::xml_node data_item_node = xml_node.append_child("DataItem");
132 assert(data_item_node);
133
134 // Add dimensions attribute
135 std::string dims;
136 for (auto d : shape)
137 dims += std::to_string(d) + std::string(" ");
138 dims.pop_back();
139 data_item_node.append_attribute("Dimensions") = dims.c_str();
140
141 // Set type for topology data (needed by XDMF to prevent default to
142 // float)
143 if (!number_type.empty())
144 data_item_node.append_attribute("NumberType") = number_type.c_str();
145
146 // Add format attribute
147 if (h5_id < 0)
148 {
149 data_item_node.append_attribute("Format") = "XML";
150 assert(shape.size() == 2);
151 std::ostringstream s;
152 s.precision(16);
153 for (std::size_t i = 0; i < x.size(); ++i)
154 {
155 if ((i + 1) % shape[1] == 0 and shape[1] != 0)
156 s << x.data()[i] << std::endl;
157 else
158 s << x.data()[i] << " ";
159 }
160
161 data_item_node.append_child(pugi::node_pcdata).set_value(s.str().c_str());
162 }
163 else
164 {
165 data_item_node.append_attribute("Format") = "HDF";
166
167 // Get name of HDF5 file, including path
168 const std::filesystem::path p = io::hdf5::get_filename(h5_id);
169 const std::filesystem::path filename = p.filename().c_str();
170
171 // Add HDF5 filename and HDF5 internal path to XML file
172 const std::string xdmf_path
173 = filename.string() + std::string(":") + h5_path;
174 data_item_node.append_child(pugi::node_pcdata).set_value(xdmf_path.c_str());
175
176 // Compute data offset and range of values
177 std::int64_t local_shape0 = std::reduce(
178 std::next(shape.begin()), shape.end(), x.size(), std::divides{});
179
180 const std::array local_range{offset, offset + local_shape0};
181 io::hdf5::write_dataset(h5_id, h5_path, x.data(), local_range, shape,
182 use_mpi_io, false);
183
184 // Add partitioning attribute to dataset
185 // std::vector<std::size_t> partitions;
186 // std::vector<std::size_t> offset_tmp(1, offset);
187 // dolfinx::MPI::gather(comm, offset_tmp, partitions);
188 // dolfinx::MPI::broadcast(comm, partitions);
189 // io::hdf5::add_attribute(h5_id, h5_path, "partition", partitions);
190 }
191}
192
197template <typename T>
198std::vector<T> get_dataset(MPI_Comm comm, const pugi::xml_node& dataset_node,
199 hid_t h5_id,
200 std::array<std::int64_t, 2> range = {0, 0})
201{
202 // FIXME: Need to sort out dataset dimensions - can't depend on HDF5
203 // shape, and a Topology data item is not required to have a
204 // 'Dimensions' attribute since the dimensions can be determined from
205 // the number of cells and the cell type (for topology, one must
206 // supply cell type + (number of cells or dimensions)).
207 //
208 // A geometry data item must have 'Dimensions' attribute.
209
210 assert(dataset_node);
211 pugi::xml_attribute format_attr = dataset_node.attribute("Format");
212 assert(format_attr);
213
214 // Get data set shape from 'Dimensions' attribute (empty if not
215 // available)
216 const std::vector shape_xml = xdmf_utils::get_dataset_shape(dataset_node);
217
218 const std::string format = format_attr.as_string();
219 std::vector<T> data_vector;
220 // Only read ASCII on process 0
221 const int mpi_rank = dolfinx::MPI::rank(comm);
222 if (format == "XML")
223 {
224 if (mpi_rank == 0)
225 {
226 // Read data and trim any leading/trailing whitespace
227 pugi::xml_node data_node = dataset_node.first_child();
228 assert(data_node);
229 std::string data_str = data_node.value();
230
231 // Split data based on spaces and line breaks
232 std::vector<boost::iterator_range<std::string::iterator>> data_vector_str;
233 boost::split(data_vector_str, data_str, boost::is_any_of(" \n"));
234
235 // Add data to numerical vector
236 data_vector.reserve(data_vector_str.size());
237 for (auto& v : data_vector_str)
238 {
239 if (v.begin() != v.end())
240 data_vector.push_back(
241 boost::lexical_cast<T>(boost::copy_range<std::string>(v)));
242 }
243 }
244 }
245 else if (format == "HDF")
246 {
247 // Get file and data path
248 auto paths = xdmf_utils::get_hdf5_paths(dataset_node);
249
250 // Get data shape from HDF5 file
251 const std::vector shape_hdf5 = io::hdf5::get_dataset_shape(h5_id, paths[1]);
252
253 // FIXME: should we support empty data sets?
254 // Check that data set is not empty
255 assert(!shape_hdf5.empty());
256 assert(shape_hdf5[0] != 0);
257
258 // Determine range of data to read from HDF5 file. This is
259 // complicated by the XML Dimension attribute and the HDF5 storage
260 // possibly having different shapes, e.g. the HDF5 storage may be a
261 // flat array.
262
263 // If range = {0, 0} then no range is supplied and we must determine
264 // the range
265 if (range[0] == 0 and range[1] == 0)
266 {
267 if (shape_xml == shape_hdf5)
268 {
269 range = dolfinx::MPI::local_range(mpi_rank, shape_hdf5[0],
270 dolfinx::MPI::size(comm));
271 }
272 else if (!shape_xml.empty() and shape_hdf5.size() == 1)
273 {
274 // Size of dims > 0
275 std::int64_t d = std::reduce(shape_xml.begin(), shape_xml.end(),
276 std::int64_t(1), std::multiplies{});
277
278 // Check for data size consistency
279 if (d * shape_xml[0] != shape_hdf5[0])
280 {
281 throw std::runtime_error("Data size in XDMF/XML and size of HDF5 "
282 "dataset are inconsistent");
283 }
284
285 // Compute data range to read
286 range = dolfinx::MPI::local_range(mpi_rank, shape_xml[0],
287 dolfinx::MPI::rank(comm));
288 range[0] *= d;
289 range[1] *= d;
290 }
291 else
292 {
293 throw std::runtime_error("This combination of array shapes in XDMF and "
294 "HDF5 is not supported");
295 }
296 }
297
298 // Retrieve data
299 if (hid_t dset_id = io::hdf5::open_dataset(h5_id, paths[1]);
300 dset_id == H5I_INVALID_HID)
301 throw std::runtime_error("Failed to open HDF5 global dataset.");
302 else
303 {
304 data_vector = io::hdf5::read_dataset<T>(dset_id, range, true);
305 if (herr_t err = H5Dclose(dset_id); err < 0)
306 throw std::runtime_error("Failed to close HDF5 global dataset.");
307 }
308 }
309 else
310 throw std::runtime_error("Storage format \"" + format + "\" is unknown");
311
312 // Get dimensions for consistency (if available in DataItem node)
313 if (shape_xml.empty())
314 {
315 std::int64_t size = 1;
316 for (auto dim : shape_xml)
317 size *= dim;
318
319 std::int64_t size_global = 0;
320 const std::int64_t size_local = data_vector.size();
321 MPI_Allreduce(&size_local, &size_global, 1, MPI_INT64_T, MPI_SUM, comm);
322 if (size != size_global)
323 {
324 throw std::runtime_error(
325 "Data sizes in attribute and size of data read are inconsistent");
326 }
327 }
328
329 return data_vector;
330}
331
332} // namespace io::xdmf_utils
333} // namespace dolfinx
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:91
CellType
Cell type identifier.
Definition cell_types.h:22
Top-level namespace.
Definition defines.h:12