DOLFINx 0.10.0.0
DOLFINx C++ interface
Loading...
Searching...
No Matches
VTKHDF.h
1// Copyright (C) 2024-2025 Chris 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#include "HDF5Interface.h"
8#include <algorithm>
9#include <concepts>
10#include <dolfinx/common/IndexMap.h>
11#include <dolfinx/io/cells.h>
12#include <dolfinx/mesh/Mesh.h>
13#include <dolfinx/mesh/Topology.h>
14#include <dolfinx/mesh/utils.h>
15#include <map>
16#include <vector>
17
18namespace dolfinx::io::VTKHDF
19{
27template <std::floating_point U>
28void write_mesh(std::string filename, const mesh::Mesh<U>& mesh)
29{
30 hid_t h5file = hdf5::open_file(mesh.comm(), filename, "w", true);
31
32 // Create VTKHDF group
33 hdf5::add_group(h5file, "VTKHDF");
34 hid_t vtk_group = H5Gopen(h5file, "VTKHDF", H5P_DEFAULT);
35
36 // Create "Version" attribute
37 hsize_t dims = 2;
38 hid_t space_id = H5Screate_simple(1, &dims, NULL);
39 hid_t attr_id = H5Acreate(vtk_group, "Version", H5T_NATIVE_INT32, space_id,
40 H5P_DEFAULT, H5P_DEFAULT);
41 std::array<std::int32_t, 2> version = {2, 2};
42 H5Awrite(attr_id, H5T_NATIVE_INT32, version.data());
43 H5Aclose(attr_id);
44 H5Sclose(space_id);
45
46 // Create "Type" attribute
47 space_id = H5Screate(H5S_SCALAR);
48 hid_t atype = H5Tcopy(H5T_C_S1);
49 H5Tset_size(atype, 16);
50 H5Tset_strpad(atype, H5T_STR_NULLTERM);
51 attr_id
52 = H5Acreate(vtk_group, "Type", atype, space_id, H5P_DEFAULT, H5P_DEFAULT);
53 H5Awrite(attr_id, atype, "UnstructuredGrid");
54 H5Aclose(attr_id);
55 H5Sclose(space_id);
56 H5Gclose(vtk_group);
57
58 // Extract topology information for each cell type
59 std::vector<mesh::CellType> cell_types
60 = mesh.topology()->entity_types(mesh.topology()->dim());
61
62 std::vector cell_index_maps
63 = mesh.topology()->index_maps(mesh.topology()->dim());
64 std::vector<std::int32_t> num_cells;
65 std::vector<std::int64_t> num_cells_global;
66 for (auto& im : cell_index_maps)
67 {
68 num_cells.push_back(im->size_local());
69 num_cells_global.push_back(im->size_global());
70 }
71
72 // Geometry dofmap and points
73 std::shared_ptr<const common::IndexMap> geom_imap
74 = mesh.geometry().index_map();
75 std::int64_t size_global = geom_imap->size_global();
76 std::vector<std::int64_t> geom_global_shape = {size_global, 3};
77 std::array<std::int64_t, 2> geom_irange = geom_imap->local_range();
78 hdf5::write_dataset(h5file, "/VTKHDF/Points", mesh.geometry().x().data(),
79 geom_irange, geom_global_shape, true, false);
80 hdf5::write_dataset(h5file, "VTKHDF/NumberOfPoints", &size_global, {0, 1},
81 {1}, true, false);
82
83 // Note: VTKHDF stores the cells as an adjacency list, where cell
84 // types might be jumbled up
85 std::vector<std::int64_t> topology_flattened;
86 std::vector<std::int64_t> topology_offsets;
87 std::vector<std::uint8_t> vtkcelltypes;
88 for (std::size_t i = 0; i < cell_index_maps.size(); ++i)
89 {
90 md::mdspan<const std::int32_t, md::dextents<std::size_t, 2>> g_dofmap
91 = mesh.geometry().dofmap(i);
92
93 std::vector<std::uint16_t> perm
94 = cells::perm_vtk(cell_types[i], g_dofmap.extent(1));
95 std::vector<std::uint16_t> inverse_perm = cells::transpose(perm);
96 std::vector<std::int32_t> local_dm;
97 local_dm.reserve(g_dofmap.extent(1) * num_cells[i]);
98 for (int j = 0; j < num_cells[i]; ++j)
99 for (std::size_t k = 0; k < g_dofmap.extent(1); ++k)
100 local_dm.push_back(g_dofmap(j, inverse_perm[k]));
101
102 std::vector<std::int64_t> global_dm(local_dm.size());
103 geom_imap->local_to_global(local_dm, global_dm);
104
105 topology_flattened.insert(topology_flattened.end(), global_dm.begin(),
106 global_dm.end());
107 topology_offsets.insert(topology_offsets.end(), g_dofmap.extent(0),
108 g_dofmap.extent(1));
109 vtkcelltypes.insert(
110 vtkcelltypes.end(), cell_index_maps[i]->size_local(),
111 cells::get_vtk_cell_type(cell_types[i], mesh.topology()->dim()));
112 }
113
114 // Create topo_offsets
115 std::partial_sum(topology_offsets.cbegin(), topology_offsets.cend(),
116 topology_offsets.begin());
117
118 std::vector<int> num_nodes_per_cell;
119 std::vector<std::int64_t> cell_start_pos;
120 std::vector<std::int64_t> cell_stop_pos;
121 for (std::size_t i = 0; i < cell_index_maps.size(); ++i)
122 {
123 num_nodes_per_cell.push_back(mesh.geometry().cmaps()[i].dim());
124 std::array<std::int64_t, 2> r = cell_index_maps[i]->local_range();
125 cell_start_pos.push_back(r[0]);
126 cell_stop_pos.push_back(r[1]);
127 }
128
129 // Compute overall cell offset from offsets for each cell type
130 std::int64_t offset_start_position
131 = std::accumulate(cell_start_pos.begin(), cell_start_pos.end(), 0);
132 std::int64_t offset_stop_position
133 = std::accumulate(cell_stop_pos.begin(), cell_stop_pos.end(), 0);
134
135 // Compute overall topology offset from offsets for each cell type
136 std::int64_t topology_start
137 = std::inner_product(num_nodes_per_cell.begin(), num_nodes_per_cell.end(),
138 cell_start_pos.begin(), 0);
139
140 std::transform(topology_offsets.cbegin(), topology_offsets.cend(),
141 topology_offsets.begin(),
142 [topology_start](auto x) { return x + topology_start; });
143
144 std::int64_t num_all_cells_global
145 = std::accumulate(num_cells_global.begin(), num_cells_global.end(), 0);
146 hdf5::write_dataset(h5file, "/VTKHDF/Offsets", topology_offsets.data(),
147 {offset_start_position + 1, offset_stop_position + 1},
148 {num_all_cells_global + 1}, true, false);
149
150 // Store global mesh connectivity
151 std::int64_t topology_size_global
152 = std::inner_product(num_nodes_per_cell.begin(), num_nodes_per_cell.end(),
153 num_cells_global.begin(), 0);
154
155 std::int64_t topology_stop = topology_start + topology_flattened.size();
156 hdf5::write_dataset(h5file, "/VTKHDF/Connectivity", topology_flattened.data(),
157 {topology_start, topology_stop}, {topology_size_global},
158 true, false);
159
160 // Store cell types
161 hdf5::write_dataset(h5file, "/VTKHDF/Types", vtkcelltypes.data(),
162 {offset_start_position, offset_stop_position},
163 {num_all_cells_global}, true, false);
164 hdf5::write_dataset(h5file, "/VTKHDF/NumberOfConnectivityIds",
165 &topology_size_global, {0, 1}, {1}, true, false);
166 hdf5::write_dataset(h5file, "/VTKHDF/NumberOfCells", &num_all_cells_global,
167 {0, 1}, {1}, true, false);
168 hdf5::close_file(h5file);
169}
170
180template <std::floating_point U>
181mesh::Mesh<U> read_mesh(MPI_Comm comm, std::string filename,
182 std::size_t gdim = 3)
183{
184 hid_t h5file = hdf5::open_file(comm, filename, "r", true);
185
186 std::vector<std::int64_t> shape
187 = hdf5::get_dataset_shape(h5file, "/VTKHDF/Types");
188 int rank = dolfinx::MPI::rank(comm);
189 int mpi_size = dolfinx::MPI::size(comm);
190 std::array<std::int64_t, 2> local_cell_range
191 = dolfinx::MPI::local_range(rank, shape[0], mpi_size);
192
193 hid_t dset_id = hdf5::open_dataset(h5file, "/VTKHDF/Types");
194 std::vector<std::uint8_t> types
195 = hdf5::read_dataset<std::uint8_t>(dset_id, local_cell_range, true);
196 H5Dclose(dset_id);
197
198 // Create reverse map (VTK -> DOLFINx cell type)
199 std::map<std::uint8_t, mesh::CellType> vtk_to_dolfinx;
200 {
201 for (auto type : {mesh::CellType::point, mesh::CellType::interval,
202 mesh::CellType::triangle, mesh::CellType::quadrilateral,
203 mesh::CellType::tetrahedron, mesh::CellType::prism,
204 mesh::CellType::pyramid, mesh::CellType::hexahedron})
205 {
206 vtk_to_dolfinx.insert(
207 {cells::get_vtk_cell_type(type, mesh::cell_dim(type)), type});
208 }
209 }
210
211 // Read in offsets to determine the different cell-types in the mesh
212 dset_id = hdf5::open_dataset(h5file, "/VTKHDF/Offsets");
213 std::vector<std::int64_t> offsets = hdf5::read_dataset<std::int64_t>(
214 dset_id, {local_cell_range[0], local_cell_range[1] + 1}, true);
215 H5Dclose(dset_id);
216
217 // Convert cell offsets to cell type and cell degree tuples
218 std::vector<std::array<std::uint8_t, 2>> types_unique;
219 std::vector<std::uint8_t> cell_degrees;
220 for (std::size_t i = 0; i < types.size(); ++i)
221 {
222 std::int64_t num_nodes = offsets[i + 1] - offsets[i];
223 std::uint8_t cell_degree
224 = cells::cell_degree(vtk_to_dolfinx.at(types[i]), num_nodes);
225 types_unique.push_back({types[i], cell_degree});
226 cell_degrees.push_back(cell_degree);
227 }
228 {
229 std::ranges::sort(types_unique);
230 auto [unique_end, range_end] = std::ranges::unique(types_unique);
231 types_unique.erase(unique_end, range_end);
232 }
233
234 // Share cell types with all processes to make global list of cell
235 // types
236 // FIXME: amount of data is small, but number of connections does not
237 // scale
238 int count = 2 * types_unique.size();
239 std::vector<std::int32_t> recv_count(mpi_size);
240 MPI_Allgather(&count, 1, MPI_INT32_T, recv_count.data(), 1, MPI_INT32_T,
241 comm);
242 std::vector<std::int32_t> recv_offsets(mpi_size + 1, 0);
243 std::partial_sum(recv_count.begin(), recv_count.end(),
244 recv_offsets.begin() + 1);
245
246 std::vector<std::array<std::uint8_t, 2>> recv_types;
247 {
248 std::vector<std::uint8_t> send_types;
249 for (std::array<std::uint8_t, 2> t : types_unique)
250 send_types.insert(send_types.end(), t.begin(), t.end());
251
252 std::vector<std::uint8_t> recv_types_buffer(recv_offsets.back());
253 MPI_Allgatherv(send_types.data(), send_types.size(), MPI_UINT8_T,
254 recv_types_buffer.data(), recv_count.data(),
255 recv_offsets.data(), MPI_UINT8_T, comm);
256
257 for (std::size_t i = 0; i < recv_types_buffer.size(); i += 2)
258 recv_types.push_back({recv_types_buffer[i], recv_types_buffer[i + 1]});
259
260 std::ranges::sort(recv_types);
261 auto [unique_end, range_end] = std::ranges::unique(recv_types);
262 recv_types.erase(unique_end, range_end);
263 }
264
265 // Map from VTKCellType to index in list of (cell types, degree)
266 std::map<std::array<std::uint8_t, 2>, std::int32_t> type_to_index;
267 std::vector<mesh::CellType> dolfinx_cell_type;
268 std::vector<std::uint8_t> dolfinx_cell_degree;
269 for (std::array<std::uint8_t, 2> ct : recv_types)
270 {
271 mesh::CellType cell_type = vtk_to_dolfinx.at(ct[0]);
272 type_to_index.insert({ct, dolfinx_cell_degree.size()});
273 dolfinx_cell_degree.push_back(ct[1]);
274 dolfinx_cell_type.push_back(cell_type);
275 }
276
277 dset_id = hdf5::open_dataset(h5file, "/VTKHDF/NumberOfPoints");
278 std::vector npoints = hdf5::read_dataset<std::int64_t>(dset_id, {0, 1}, true);
279 H5Dclose(dset_id);
280 spdlog::info("Mesh with {} points", npoints[0]);
281 std::array<std::int64_t, 2> local_point_range
282 = dolfinx::MPI::local_range(rank, npoints[0], mpi_size);
283
284 std::vector<std::int64_t> x_shape
285 = hdf5::get_dataset_shape(h5file, "/VTKHDF/Points");
286 dset_id = hdf5::open_dataset(h5file, "/VTKHDF/Points");
287 std::vector<U> points_local
288 = hdf5::read_dataset<U>(dset_id, local_point_range, true);
289 H5Dclose(dset_id);
290
291 // Remove coordinates if gdim != 3
292 assert(gdim <= 3);
293 std::vector<U> points_pruned((local_point_range[1] - local_point_range[0])
294 * gdim);
295 for (std::int64_t i = 0; i < local_point_range[1] - local_point_range[0]; ++i)
296 {
297 std::copy_n(points_local.begin() + i * 3, gdim,
298 points_pruned.begin() + i * gdim);
299 }
300
301 dset_id = hdf5::open_dataset(h5file, "/VTKHDF/Connectivity");
302 std::vector<std::int64_t> topology = hdf5::read_dataset<std::int64_t>(
303 dset_id, {offsets.front(), offsets.back()}, true);
304 H5Dclose(dset_id);
305 std::transform(offsets.cbegin(), offsets.cend(), offsets.begin(),
306 [offset = offsets.front()](auto x) { return x - offset; });
307 hdf5::close_file(h5file);
308
309 // Create cell topologies for each celltype in mesh
310 std::vector<std::vector<std::int64_t>> cells_local(recv_types.size());
311 for (std::size_t j = 0; j < types.size(); ++j)
312 {
313 std::int32_t type_index = type_to_index.at({types[j], cell_degrees[j]});
314 mesh::CellType cell_type = dolfinx_cell_type[type_index];
315 std::vector<std::uint16_t> perm
316 = cells::perm_vtk(cell_type, offsets[j + 1] - offsets[j]);
317 for (std::int64_t k = 0; k < offsets[j + 1] - offsets[j]; ++k)
318 cells_local[type_index].push_back(topology[perm[k] + offsets[j]]);
319 }
320
321 // Make coordinate elements
322 std::vector<fem::CoordinateElement<U>> coordinate_elements;
323 std::transform(
324 dolfinx_cell_type.cbegin(), dolfinx_cell_type.cend(),
325 dolfinx_cell_degree.cbegin(), std::back_inserter(coordinate_elements),
326 [](auto cell_type, auto cell_degree)
327 {
328 basix::element::lagrange_variant variant
329 = (cell_degree > 2) ? basix::element::lagrange_variant::equispaced
330 : basix::element::lagrange_variant::unset;
331 return fem::CoordinateElement<U>(cell_type, cell_degree, variant);
332 });
333
334 auto part = create_cell_partitioner(mesh::GhostMode::none);
335 std::vector<std::span<const std::int64_t>> cells_span(cells_local.begin(),
336 cells_local.end());
337 return mesh::create_mesh(comm, comm, cells_span, coordinate_elements, comm,
338 points_pruned, {(std::size_t)x_shape[0], gdim},
339 part);
340}
341} // namespace dolfinx::io::VTKHDF
Functions supporting mesh operations.
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
std::int8_t get_vtk_cell_type(mesh::CellType cell, int dim)
Get VTK cell identifier.
Definition cells.cpp:698
std::vector< std::uint16_t > perm_vtk(mesh::CellType type, int num_nodes)
Permutation array to map from VTK to DOLFINx node ordering.
Definition cells.cpp:516
std::vector< std::uint16_t > transpose(std::span< const std::uint16_t > map)
Compute the transpose of a re-ordering map.
Definition cells.cpp:670
int cell_degree(mesh::CellType type, int num_nodes)
Get the Lagrange order of a given cell with a given number of nodes.
Definition cells.cpp:592
int cell_dim(CellType type)
Return topological dimension of cell type.
Definition cell_types.cpp:134
CellType
Cell type identifier.
Definition cell_types.h:22
Mesh< typename std::remove_reference_t< typename U::value_type > > create_mesh(MPI_Comm comm, MPI_Comm commt, std::vector< std::span< const std::int64_t > > cells, const std::vector< fem::CoordinateElement< typename std::remove_reference_t< typename U::value_type > > > &elements, MPI_Comm commg, const U &x, std::array< std::size_t, 2 > xshape, const CellPartitionFunction &partitioner, const CellReorderFunction &reorder_fn=graph::reorder_gps)
Create a distributed mesh::Mesh from mesh data and using the provided graph partitioning function for...
Definition utils.h:848
CellPartitionFunction create_cell_partitioner(mesh::GhostMode ghost_mode=mesh::GhostMode::none, const graph::partition_fn &partfn=&graph::partition_graph)
Create a function that computes destination rank for mesh cells on this rank by applying the default ...
Definition utils.cpp:100
std::string version()
Return DOLFINx version string.
Definition defines.cpp:10