Loading...
Searching...
No Matches
VTKHDF.h
// Copyright (C) 2024-2025 Chris Richardson, Jørgen S. Dokken
//
// This file is part of DOLFINx (https://www.fenicsproject.org)
//
// SPDX-License-Identifier: LGPL-3.0-or-later
#include "HDF5Interface.h"
#include <algorithm>
#include <concepts>
#include <dolfinx/common/IndexMap.h>
#include <dolfinx/io/cells.h>
#include <dolfinx/mesh/Mesh.h>
#include <dolfinx/mesh/Topology.h>
#include <dolfinx/mesh/utils.h>
#include <map>
#include <vector>
namespace dolfinx::io::VTKHDF
{
template <std::floating_point U>
void write_mesh(std::string filename, const mesh::Mesh<U>& mesh)
{
hid_t h5file = hdf5::open_file(mesh.comm(), filename, "w", true);
// Create VTKHDF group
hdf5::add_group(h5file, "VTKHDF");
hid_t vtk_group = H5Gopen(h5file, "VTKHDF", H5P_DEFAULT);
hdf5::set_attribute(vtk_group, "Version", std::vector{2, 2});
hdf5::set_attribute(vtk_group, "Type", "UnstructuredGrid");
H5Gclose(vtk_group);
// Extract topology information for each cell type
std::vector<mesh::CellType> cell_types
= mesh.topology()->entity_types(mesh.topology()->dim());
std::vector cell_index_maps
= mesh.topology()->index_maps(mesh.topology()->dim());
std::vector<std::int32_t> num_cells;
std::vector<std::int64_t> num_cells_global;
for (auto& im : cell_index_maps)
{
num_cells.push_back(im->size_local());
num_cells_global.push_back(im->size_global());
}
// Geometry dofmap and points
std::shared_ptr<const common::IndexMap> geom_imap
= mesh.geometry().index_map();
std::int64_t size_global = geom_imap->size_global();
std::vector<std::int64_t> geom_global_shape = {size_global, 3};
std::array<std::int64_t, 2> geom_irange = geom_imap->local_range();
hdf5::write_dataset(h5file, "/VTKHDF/Points", mesh.geometry().x().data(),
geom_irange, geom_global_shape, true, false);
hdf5::write_dataset(h5file, "VTKHDF/NumberOfPoints", &size_global, {0, 1},
{1}, true, false);
// Note: VTKHDF stores the cells as an adjacency list, where cell
// types might be jumbled up
std::vector<std::int64_t> topology_flattened;
std::vector<std::int64_t> topology_offsets;
std::vector<std::uint8_t> vtkcelltypes;
for (std::size_t i = 0; i < cell_index_maps.size(); ++i)
{
md::mdspan<const std::int32_t, md::dextents<std::size_t, 2>> g_dofmap
= mesh.geometry().dofmap(i);
std::vector<std::uint16_t> perm
= cells::perm_vtk(cell_types[i], g_dofmap.extent(1));
std::vector<std::uint16_t> inverse_perm = cells::transpose(perm);
std::vector<std::int32_t> local_dm;
local_dm.reserve(g_dofmap.extent(1) * num_cells[i]);
for (int j = 0; j < num_cells[i]; ++j)
for (std::size_t k = 0; k < g_dofmap.extent(1); ++k)
local_dm.push_back(g_dofmap(j, inverse_perm[k]));
std::vector<std::int64_t> global_dm(local_dm.size());
geom_imap->local_to_global(local_dm, global_dm);
topology_flattened.insert(topology_flattened.end(), global_dm.begin(),
global_dm.end());
topology_offsets.insert(topology_offsets.end(), g_dofmap.extent(0),
g_dofmap.extent(1));
vtkcelltypes.insert(
vtkcelltypes.end(), cell_index_maps[i]->size_local(),
cells::get_vtk_cell_type(cell_types[i], mesh.topology()->dim()));
}
// Create topo_offsets
std::partial_sum(topology_offsets.cbegin(), topology_offsets.cend(),
topology_offsets.begin());
std::vector<int> num_nodes_per_cell;
std::vector<std::int64_t> cell_start_pos;
std::vector<std::int64_t> cell_stop_pos;
for (std::size_t i = 0; i < cell_index_maps.size(); ++i)
{
num_nodes_per_cell.push_back(mesh.geometry().cmaps()[i].dim());
std::array<std::int64_t, 2> r = cell_index_maps[i]->local_range();
cell_start_pos.push_back(r[0]);
cell_stop_pos.push_back(r[1]);
}
// Compute overall cell offset from offsets for each cell type
std::int64_t offset_start_position
= std::accumulate(cell_start_pos.begin(), cell_start_pos.end(), 0);
std::int64_t offset_stop_position
= std::accumulate(cell_stop_pos.begin(), cell_stop_pos.end(), 0);
// Compute overall topology offset from offsets for each cell type
std::int64_t topology_start
= std::inner_product(num_nodes_per_cell.begin(), num_nodes_per_cell.end(),
cell_start_pos.begin(), 0);
std::transform(topology_offsets.cbegin(), topology_offsets.cend(),
topology_offsets.begin(),
[topology_start](auto x) { return x + topology_start; });
std::int64_t num_all_cells_global
= std::accumulate(num_cells_global.begin(), num_cells_global.end(), 0);
hdf5::write_dataset(h5file, "/VTKHDF/Offsets", topology_offsets.data(),
{offset_start_position + 1, offset_stop_position + 1},
{num_all_cells_global + 1}, true, false);
// Store global mesh connectivity
std::int64_t topology_size_global
= std::inner_product(num_nodes_per_cell.begin(), num_nodes_per_cell.end(),
num_cells_global.begin(), 0);
std::int64_t topology_stop = topology_start + topology_flattened.size();
hdf5::write_dataset(h5file, "/VTKHDF/Connectivity", topology_flattened.data(),
{topology_start, topology_stop}, {topology_size_global},
true, false);
// Store cell types
hdf5::write_dataset(h5file, "/VTKHDF/Types", vtkcelltypes.data(),
{offset_start_position, offset_stop_position},
{num_all_cells_global}, true, false);
hdf5::write_dataset(h5file, "/VTKHDF/NumberOfConnectivityIds",
&topology_size_global, {0, 1}, {1}, true, false);
hdf5::write_dataset(h5file, "/VTKHDF/NumberOfCells", &num_all_cells_global,
{0, 1}, {1}, true, false);
hdf5::close_file(h5file);
}
template <std::floating_point U>
void write_data(std::string point_or_cell, std::string filename,
const mesh::Mesh<U>& mesh, const std::vector<U>& data,
double time)
{
std::vector<std::shared_ptr<const common::IndexMap>> index_maps;
if (point_or_cell == "Point")
index_maps = {mesh.geometry().index_map()};
else if (point_or_cell == "Cell")
index_maps = mesh.topology()->index_maps(mesh.topology()->dim());
else
throw std::runtime_error("Selection must be Point or Cell");
std::string dataset_name = "/VTKHDF/" + point_or_cell + "Data/u";
int npoints
= std::accumulate(index_maps.begin(), index_maps.end(), 0,
[](int a, auto im) { return a + im->size_local(); });
int data_width = data.size() / npoints;
if (data.size() % npoints != 0)
{
throw std::runtime_error(
"Data size mismatch with number of local vertices/cells");
}
spdlog::debug("Data vector width={}", data_width);
hid_t h5file = hdf5::open_file(mesh.comm(), filename, "a", true);
hdf5::add_group(h5file, "VTKHDF/Steps");
hid_t vtk_group = H5Gopen(h5file, "VTKHDF/Steps", H5P_DEFAULT);
std::int64_t point_data_offset = 0;
if (htri_t attr_exists = H5Aexists(vtk_group, "NSteps"); attr_exists < 0)
throw std::runtime_error("Error checking attribute");
else if (attr_exists == 0)
hdf5::set_attribute(vtk_group, "NSteps", 1);
else
{
// Read and increment attribute
std::int32_t nsteps = 0;
hid_t attr_id = H5Aopen(vtk_group, "NSteps", H5P_DEFAULT);
H5Aread(attr_id, H5T_NATIVE_INT32, &nsteps);
nsteps++;
H5Awrite(attr_id, H5T_NATIVE_INT32, &nsteps);
H5Aclose(attr_id);
std::vector<std::int64_t> data_shape
= hdf5::get_dataset_shape(h5file, dataset_name);
assert(data_shape.size() == 2);
point_data_offset = data_shape[0];
}
H5Gclose(vtk_group);
// Add a single value to end of a 1D dataset
auto append_dataset
= [&h5file]<typename T>(const std::string& dset_name, T value)
{
std::int32_t s = 0;
if (hdf5::has_dataset(h5file, dset_name))
{
std::vector<std::int64_t> shape
= hdf5::get_dataset_shape(h5file, dset_name);
assert(shape.size() == 1);
s = shape[0];
}
hdf5::write_dataset(h5file, dset_name, &value, {s, s + 1}, {s + 1}, true,
true);
};
// Mesh remains the same, so these values are the same for each time step
append_dataset("/VTKHDF/Steps/CellOffsets", 0);
append_dataset("/VTKHDF/Steps/ConnectivityIdOffsets", 0);
append_dataset("/VTKHDF/Steps/NumberOfParts", 1);
append_dataset("/VTKHDF/Steps/PartOffsets", 0);
append_dataset("/VTKHDF/Steps/PointOffsets", 0);
// Add the current data size to the end of the offset array
hdf5::add_group(h5file, "/VTKHDF/Steps/" + point_or_cell + "DataOffsets");
append_dataset("/VTKHDF/Steps/" + point_or_cell + "DataOffsets/u",
point_data_offset);
// Time values
// FIXME: check these are increasing?
append_dataset("/VTKHDF/Steps/Values", time);
std::string group_name = "/VTKHDF/" + point_or_cell + "Data";
hdf5::add_group(h5file, group_name);
// Add point/cell data into dataset, extending each time by
// size_global with each process writing its own part.
std::int64_t range0 = std::accumulate(index_maps.begin(), index_maps.end(), 0,
[](int a, auto im)
{ return a + im->local_range()[0]; });
std::array<std::int64_t, 2> range{range0, range0 + npoints};
std::int64_t global_size = std::accumulate(
index_maps.begin(), index_maps.end(), 0,
[](std::int64_t a, auto im) { return a + im->size_global(); });
std::vector<std::int64_t> shape0 = {global_size, data_width};
if (hdf5::has_dataset(h5file, dataset_name))
{
std::vector<std::int64_t> shape
= hdf5::get_dataset_shape(h5file, dataset_name);
assert(shape.size() == 2);
std::int64_t offset = shape[0];
range[0] += offset;
range[1] += offset;
shape0[0] += offset;
hdf5::write_dataset(h5file, dataset_name, data.data(), range, shape0, true,
true);
}
else
{
hdf5::write_dataset(h5file, dataset_name, data.data(), range, shape0, true,
true);
if (data_width > 1)
{
hid_t dset_id = hdf5::open_dataset(h5file, dataset_name);
hdf5::set_attribute(dset_id, "NumberOfComponents", data_width);
H5Dclose(dset_id);
hid_t vtk_group = H5Gopen(h5file, group_name.c_str(), H5P_DEFAULT);
hdf5::set_attribute(vtk_group, "Vectors", "u");
H5Gclose(vtk_group);
}
}
hdf5::close_file(h5file);
}
template <std::floating_point U>
mesh::Mesh<U> read_mesh(MPI_Comm comm, std::string filename,
std::size_t gdim = 3,
std::optional<std::int32_t> max_facet_to_cell_links = 2)
{
hid_t h5file = hdf5::open_file(comm, filename, "r", true);
std::vector<std::int64_t> shape
= hdf5::get_dataset_shape(h5file, "/VTKHDF/Types");
int rank = dolfinx::MPI::rank(comm);
int mpi_size = dolfinx::MPI::size(comm);
std::array<std::int64_t, 2> local_cell_range
= dolfinx::MPI::local_range(rank, shape[0], mpi_size);
hid_t dset_id = hdf5::open_dataset(h5file, "/VTKHDF/Types");
std::vector<std::uint8_t> types
= hdf5::read_dataset<std::uint8_t>(dset_id, local_cell_range, true);
H5Dclose(dset_id);
// Create reverse map (VTK -> DOLFINx cell type)
std::map<std::uint8_t, mesh::CellType> vtk_to_dolfinx;
{
for (auto type : {mesh::CellType::point, mesh::CellType::interval,
mesh::CellType::triangle, mesh::CellType::quadrilateral,
mesh::CellType::tetrahedron, mesh::CellType::prism,
mesh::CellType::pyramid, mesh::CellType::hexahedron})
{
vtk_to_dolfinx.insert(
{cells::get_vtk_cell_type(type, mesh::cell_dim(type)), type});
}
}
// Read in offsets to determine the different cell-types in the mesh
dset_id = hdf5::open_dataset(h5file, "/VTKHDF/Offsets");
std::vector<std::int64_t> offsets = hdf5::read_dataset<std::int64_t>(
dset_id, {local_cell_range[0], local_cell_range[1] + 1}, true);
H5Dclose(dset_id);
// Convert cell offsets to cell type and cell degree tuples
std::vector<std::array<std::uint8_t, 2>> types_unique;
std::vector<std::uint8_t> cell_degrees;
for (std::size_t i = 0; i < types.size(); ++i)
{
std::int64_t num_nodes = offsets[i + 1] - offsets[i];
std::uint8_t cell_degree
= cells::cell_degree(vtk_to_dolfinx.at(types[i]), num_nodes);
types_unique.push_back({types[i], cell_degree});
cell_degrees.push_back(cell_degree);
}
{
std::ranges::sort(types_unique);
auto [unique_end, range_end] = std::ranges::unique(types_unique);
types_unique.erase(unique_end, range_end);
}
// Share cell types with all processes to make global list of cell
// types
// FIXME: amount of data is small, but number of connections does not
// scale
int count = 2 * types_unique.size();
std::vector<std::int32_t> recv_count(mpi_size);
MPI_Allgather(&count, 1, MPI_INT32_T, recv_count.data(), 1, MPI_INT32_T,
comm);
std::vector<std::int32_t> recv_offsets(mpi_size + 1, 0);
std::partial_sum(recv_count.begin(), recv_count.end(),
recv_offsets.begin() + 1);
std::vector<std::array<std::uint8_t, 2>> recv_types;
{
std::vector<std::uint8_t> send_types;
for (std::array<std::uint8_t, 2> t : types_unique)
send_types.insert(send_types.end(), t.begin(), t.end());
std::vector<std::uint8_t> recv_types_buffer(recv_offsets.back());
MPI_Allgatherv(send_types.data(), send_types.size(), MPI_UINT8_T,
recv_types_buffer.data(), recv_count.data(),
recv_offsets.data(), MPI_UINT8_T, comm);
for (std::size_t i = 0; i < recv_types_buffer.size(); i += 2)
recv_types.push_back({recv_types_buffer[i], recv_types_buffer[i + 1]});
std::ranges::sort(recv_types);
auto [unique_end, range_end] = std::ranges::unique(recv_types);
recv_types.erase(unique_end, range_end);
}
// Map from VTKCellType to index in list of (cell types, degree)
std::map<std::array<std::uint8_t, 2>, std::int32_t> type_to_index;
std::vector<mesh::CellType> dolfinx_cell_type;
std::vector<std::uint8_t> dolfinx_cell_degree;
for (std::array<std::uint8_t, 2> ct : recv_types)
{
mesh::CellType cell_type = vtk_to_dolfinx.at(ct[0]);
type_to_index.insert({ct, dolfinx_cell_degree.size()});
dolfinx_cell_degree.push_back(ct[1]);
dolfinx_cell_type.push_back(cell_type);
}
dset_id = hdf5::open_dataset(h5file, "/VTKHDF/NumberOfPoints");
std::vector npoints = hdf5::read_dataset<std::int64_t>(dset_id, {0, 1}, true);
H5Dclose(dset_id);
spdlog::info("Mesh with {} points", npoints[0]);
std::array<std::int64_t, 2> local_point_range
= dolfinx::MPI::local_range(rank, npoints[0], mpi_size);
std::vector<std::int64_t> x_shape
= hdf5::get_dataset_shape(h5file, "/VTKHDF/Points");
dset_id = hdf5::open_dataset(h5file, "/VTKHDF/Points");
std::vector<U> points_local
= hdf5::read_dataset<U>(dset_id, local_point_range, true);
H5Dclose(dset_id);
// Remove coordinates if gdim != 3
assert(gdim <= 3);
std::vector<U> points_pruned((local_point_range[1] - local_point_range[0])
* gdim);
for (std::int64_t i = 0; i < local_point_range[1] - local_point_range[0]; ++i)
{
std::copy_n(points_local.begin() + i * 3, gdim,
points_pruned.begin() + i * gdim);
}
dset_id = hdf5::open_dataset(h5file, "/VTKHDF/Connectivity");
std::vector<std::int64_t> topology = hdf5::read_dataset<std::int64_t>(
dset_id, {offsets.front(), offsets.back()}, true);
H5Dclose(dset_id);
std::transform(offsets.cbegin(), offsets.cend(), offsets.begin(),
[offset = offsets.front()](auto x) { return x - offset; });
hdf5::close_file(h5file);
// Create cell topologies for each celltype in mesh
std::vector<std::vector<std::int64_t>> cells_local(recv_types.size());
for (std::size_t j = 0; j < types.size(); ++j)
{
std::int32_t type_index = type_to_index.at({types[j], cell_degrees[j]});
mesh::CellType cell_type = dolfinx_cell_type[type_index];
std::vector<std::uint16_t> perm
= cells::perm_vtk(cell_type, offsets[j + 1] - offsets[j]);
for (std::int64_t k = 0; k < offsets[j + 1] - offsets[j]; ++k)
cells_local[type_index].push_back(topology[perm[k] + offsets[j]]);
}
// Make coordinate elements
std::vector<fem::CoordinateElement<U>> coordinate_elements;
std::transform(
dolfinx_cell_type.cbegin(), dolfinx_cell_type.cend(),
dolfinx_cell_degree.cbegin(), std::back_inserter(coordinate_elements),
[](auto cell_type, auto cell_degree)
{
basix::element::lagrange_variant variant
= (cell_degree > 2) ? basix::element::lagrange_variant::equispaced
: basix::element::lagrange_variant::unset;
return fem::CoordinateElement<U>(cell_type, cell_degree, variant);
});
auto part = create_cell_partitioner(mesh::GhostMode::none,
dolfinx::graph::partition_graph,
max_facet_to_cell_links);
std::vector<std::span<const std::int64_t>> cells_span(cells_local.begin(),
cells_local.end());
return mesh::create_mesh(comm, comm, cells_span, coordinate_elements, comm,
points_pruned, {(std::size_t)x_shape[0], gdim}, part,
max_facet_to_cell_links);
}
} // namespace dolfinx::io::VTKHDF
utils.h
Functions supporting mesh operations.
dolfinx::MPI::size
int size(MPI_Comm comm)
Definition MPI.cpp:72
dolfinx::MPI::rank
int rank(MPI_Comm comm)
Return process rank for the communicator.
Definition MPI.cpp:64
dolfinx::MPI::local_range
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:89
dolfinx::graph::partition_graph
AdjacencyList< std::int32_t > partition_graph(MPI_Comm comm, int nparts, const AdjacencyList< std::int64_t > &local_graph, bool ghosting)
Partition graph across processes using the default graph partitioner.
Definition partition.cpp:21
dolfinx::io::cells::get_vtk_cell_type
std::int8_t get_vtk_cell_type(mesh::CellType cell, int dim)
Get VTK cell identifier.
Definition cells.cpp:707
dolfinx::io::cells::perm_vtk
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:525
dolfinx::io::cells::transpose
std::vector< std::uint16_t > transpose(std::span< const std::uint16_t > map)
Compute the transpose of a re-ordering map.
Definition cells.cpp:679
dolfinx::io::cells::cell_degree
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:601
dolfinx::mesh::create_cell_partitioner
CellPartitionFunction create_cell_partitioner(mesh::GhostMode ghost_mode=mesh::GhostMode::none, const graph::partition_fn &partfn=&graph::partition_graph, std::optional< std::int32_t > max_facet_to_cell_links=2)
Create a function that computes destination rank for mesh cells on this rank by applying the default ...
Definition utils.cpp:100
dolfinx::mesh::cell_dim
int cell_dim(CellType type)
Return topological dimension of cell type.
Definition cell_types.cpp:134
dolfinx::mesh::CellType
CellType
Cell type identifier.
Definition cell_types.h:22
dolfinx::mesh::create_mesh
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, std::optional< std::int32_t > max_facet_to_cell_links, 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:1025