VTKHDF.h

// Copyright (C) 2024-2025 Chris Richardson, Jørgen S. Dokken
// Copyright (C) 2026 Garth N. Wells
//
// This file is part of DOLFINx (https://www.fenicsproject.org)
//
// SPDX-License-Identifier:    LGPL-3.0-or-later
 
#pragma once
 
#include "HDF5Interface.h"
#include <algorithm>
#include <array>
#include <basix/cell.h>
#include <basix/finite-element.h>
#include <basix/mdspan.hpp>
#include <cassert>
#include <concepts>
#include <cstddef>
#include <cstdint>
#include <dolfinx/common/IndexMap.h>
#include <dolfinx/common/MPI.h>
#include <dolfinx/common/log.h>
#include <dolfinx/fem/CoordinateElement.h>
#include <dolfinx/graph/partition.h>
#include <dolfinx/io/cells.h>
#include <dolfinx/mesh/Mesh.h>
#include <dolfinx/mesh/Topology.h>
#include <dolfinx/mesh/cell_types.h>
#include <dolfinx/mesh/utils.h>
#include <filesystem>
#include <format>
#include <iterator>
#include <limits>
#include <map>
#include <memory>
#include <numeric>
#include <optional>
#include <span>
#include <stdexcept>
#include <string>
#include <string_view>
#include <tuple>
#include <utility>
#include <vector>
 
namespace dolfinx::io::VTKHDF
{
namespace impl
{
// Construct a nodal coordinate element for a (VTK type, degree) group.
template <std::floating_point U>
fem::CoordinateElement<U>
coordinate_element(std::array<std::uint8_t, 2> type_degree)
{
  mesh::CellType cell_type = std::get<0>(cells::vtk_to_dolfinx(type_degree[0]));
 
  // VTK_QUADRATIC_WEDGE and VTK_QUADRATIC_PYRAMID. Both are serendipity
  // layouts (15 and 13 nodes) that Basix has no element for, and
  // vtk_to_dolfinx reports degree 2 for them, so they would otherwise
  // silently become 18- and 14-node Lagrange cells read from a 15- and
  // 13-node connectivity. Every rank sees the same cell types, so this
  // throws collectively.
  if (type_degree[0] == 26 or type_degree[0] == 27)
  {
    throw std::runtime_error(std::format(
        "VTK cell type {} is a serendipity layout that Basix cannot express.",
        static_cast<int>(type_degree[0])));
  }
 
  // VTK_QUADRATIC_QUAD and VTK_QUADRATIC_HEXAHEDRON, the two quadratic
  // serendipity layouts that Basix has an element for
  if (type_degree[0] == 23 or type_degree[0] == 25)
  {
    basix::cell::type cell = mesh::cell_type_to_basix_type(cell_type);
    basix::FiniteElement<U> element
        = basix::create_element<U>(basix::element::family::serendipity, cell, 2,
                                   basix::element::lagrange_variant::unset,
                                   basix::element::dpc_variant::unset, false);
    auto x = element.x();
    auto M = element.M();
    auto [geometry, geometry_shape] = basix::cell::geometry<U>(cell);
    const std::vector<std::vector<int>> edges = basix::cell::topology(cell)[1];
 
    // VTK stores edge midpoint values, not Basix's edge moments.
    for (std::size_t i = 0; i < edges.size(); ++i)
    {
      std::vector<U> midpoint(geometry_shape[1]);
      for (std::size_t j = 0; j < geometry_shape[1]; ++j)
      {
        midpoint[j] = (geometry[edges[i][0] * geometry_shape[1] + j]
                       + geometry[edges[i][1] * geometry_shape[1] + j])
                      / 2;
      }
      x[1][i] = {std::move(midpoint), {1, geometry_shape[1]}};
      M[1][i] = {{1}, {1, 1, 1, 1}};
    }
    std::array<std::vector<basix::impl::mdspan_t<const U, 2>>, 4> x_view;
    std::array<std::vector<basix::impl::mdspan_t<const U, 4>>, 4> M_view;
    for (std::size_t d = 0; d < 4; ++d)
    {
      for (const auto& [data, shape] : x[d])
        x_view[d].emplace_back(data.data(), shape);
      for (const auto& [data, shape] : M[d])
        M_view[d].emplace_back(data.data(), shape);
    }
    const auto& [wcoeffs, wshape] = element.wcoeffs();
    return fem::CoordinateElement<U>(std::make_shared<basix::FiniteElement<U>>(
        basix::create_custom_element<U>(
            cell, {}, {wcoeffs.data(), wshape}, x_view, M_view, 0,
            element.map_type(), element.sobolev_space(), false,
            element.embedded_subdegree(), element.embedded_superdegree(),
            element.polyset_type())));
  }
  int cell_degree = type_degree[1];
  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);
}
} // namespace impl

template <std::floating_point U>
void write_mesh(const std::filesystem::path& filename,
                const mesh::Mesh<U>& mesh)
{
  const int tdim = mesh.topology()->dim();
  std::vector<mesh::CellType> cell_types = mesh.topology()->entity_types(tdim);
 
  // Mapped before the file is opened so that an unsupported cell layout
  // throws on every rank before any collective IO has happened
  std::vector<std::int8_t> vtk_types;
  vtk_types.reserve(cell_types.size());
  for (std::size_t i = 0; i < cell_types.size(); ++i)
  {
    vtk_types.push_back(cells::get_vtk_cell_type(
        cell_types[i], mesh.geometry().cmaps().at(i).dim()));
  }
 
  hdf5::Handle h5file(hdf5::open_file(mesh.comm(), filename, "w", true),
                      H5Fclose);
 
  // Create VTKHDF group
  hdf5::add_group(h5file, "VTKHDF");
  {
    hdf5::Handle vtk_group(H5Gopen(h5file, "VTKHDF", H5P_DEFAULT), H5Gclose);
    if (vtk_group < 0)
      throw std::runtime_error("Failed to open HDF5 group \"VTKHDF\".");
    hdf5::set_attribute(vtk_group, "Version", std::vector{2, 2});
    hdf5::set_attribute(vtk_group, "Type", "UnstructuredGrid");
  }
 
  // Extract topology information for each cell type
  std::vector cell_index_maps = mesh.topology()->index_maps(tdim);
  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().dofmaps().at(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(),
                        vtk_types[i]);
  }
 
  // 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().at(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(), std::int64_t(0));
  std::int64_t offset_stop_position = std::accumulate(
      cell_stop_pos.begin(), cell_stop_pos.end(), std::int64_t(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(), std::int64_t(0));
 
  std::ranges::transform(topology_offsets, topology_offsets.begin(),
                         [topology_start](std::int64_t x)
                         { return x + topology_start; });
 
  std::int64_t num_all_cells_global = std::accumulate(
      num_cells_global.begin(), num_cells_global.end(), std::int64_t(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(), std::int64_t(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);
  h5file.close();
}

template <std::floating_point U>
void write_data(std::string_view point_or_cell,
                const std::filesystem::path& 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::invalid_argument("Selection must be Point or Cell.");
 
  const std::string poc(point_or_cell);
  std::string dataset_name = std::format("/VTKHDF/{}Data/u", poc);
  std::int32_t npoints = std::accumulate(
      index_maps.begin(), index_maps.end(), std::int32_t(0),
      [](std::int32_t a, auto& im) { return a + im->size_local(); });
  if (npoints == 0 and !data.empty())
  {
    throw std::invalid_argument(
        "Data supplied on a process with no local vertices/cells.");
  }
  if (npoints > 0 and data.size() % npoints != 0)
  {
    throw std::invalid_argument(
        "Data size mismatch with number of local vertices/cells.");
  }
 
  // A process with no local entities cannot determine the data width
  // from its (empty) data, but every process must build the same global
  // dataset shape. Take the width from the processes that do have data,
  // reducing (max width, -min width) so that an inconsistent width is
  // detected identically on all processes.
  const std::int32_t width_local = npoints == 0 ? 0 : data.size() / npoints;
  std::array<std::int32_t, 2> width_local_pair{
      width_local,
      npoints == 0 ? std::numeric_limits<std::int32_t>::min() : -width_local};
  std::array<std::int32_t, 2> width{0, 0};
  MPI_Allreduce(width_local_pair.data(), width.data(), 2, MPI_INT32_T, MPI_MAX,
                mesh.comm());
  const std::int32_t data_width = width[0];
  if (data_width == 0 or data_width != -width[1])
  {
    throw std::runtime_error(
        "Data width is zero, or differs between processes.");
  }
  logger().debug("Data vector width={}", data_width);
 
  hdf5::Handle h5file(hdf5::open_file(mesh.comm(), filename, "a", true),
                      H5Fclose);
  hdf5::add_group(h5file, "VTKHDF/Steps");
 
  std::int64_t point_data_offset = 0;
  {
    hdf5::Handle vtk_group(H5Gopen(h5file, "VTKHDF/Steps", H5P_DEFAULT),
                           H5Gclose);
    if (vtk_group < 0)
      throw std::runtime_error("Failed to open HDF5 group \"VTKHDF/Steps\".");
 
    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;
      {
        hdf5::Handle attr_id(H5Aopen(vtk_group, "NSteps", H5P_DEFAULT),
                             H5Aclose);
        if (attr_id < 0)
          throw std::runtime_error("Failed to open HDF5 attribute \"NSteps\".");
        if (H5Aread(attr_id, H5T_NATIVE_INT32, &nsteps) < 0)
          throw std::runtime_error("Failed to read HDF5 attribute \"NSteps\".");
      }
      hdf5::set_attribute(vtk_group, "NSteps", nsteps + 1);
 
      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];
    }
  }
 
  // 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, std::format("/VTKHDF/Steps/{}DataOffsets", poc));
  append_dataset(std::format("/VTKHDF/Steps/{}DataOffsets/u", poc),
                 point_data_offset);
 
  // Time values
  // FIXME: check these are increasing?
  append_dataset("/VTKHDF/Steps/Values", time);
 
  std::string group_name = std::format("/VTKHDF/{}Data", poc);
  hdf5::add_group(h5file, group_name);
 
  // Add point/cell data into dataset, extending each time by
  // global_size with each process writing its own part.
  std::int64_t range0 = std::accumulate(
      index_maps.begin(), index_maps.end(), std::int64_t(0),
      [](std::int64_t 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(), std::int64_t(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)
    {
      hdf5::Handle dset_id(hdf5::open_dataset(h5file, dataset_name), H5Dclose);
      hdf5::set_attribute(dset_id, "NumberOfComponents", data_width);
      hdf5::Handle vec_group(H5Gopen(h5file, group_name.c_str(), H5P_DEFAULT),
                             H5Gclose);
      if (vec_group < 0)
        throw std::runtime_error("Failed to open HDF5 group for vector data.");
      hdf5::set_attribute(vec_group, "Vectors", "u");
    }
  }
 
  h5file.close();
}

template <std::floating_point U>
mesh::Mesh<U> read_mesh(MPI_Comm comm, const std::filesystem::path& filename,
                        std::size_t gdim = 3,
                        std::optional<std::int32_t> max_facet_to_cell_links = 2,
                        int num_threads = 1)
{
  if (gdim > 3)
  {
    throw std::invalid_argument(
        "Geometric dimension must be less than or equal to 3.");
  }
 
  hdf5::Handle h5file(hdf5::open_file(comm, filename, "r", true), H5Fclose);
 
  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
      = common::local_range(rank, shape[0], mpi_size);
 
  std::vector<std::uint8_t> types = hdf5::read_dataset<std::uint8_t>(
      h5file, "/VTKHDF/Types", local_cell_range, true);
 
  // Read in offsets to determine the different cell-types in the mesh
  std::vector<std::int64_t> offsets = hdf5::read_dataset<std::int64_t>(
      h5file, "/VTKHDF/Offsets", {local_cell_range[0], local_cell_range[1] + 1},
      true);
 
  // The VTK type identifies the family; degree distinguishes Lagrange layouts.
  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];
    auto [cell_type, degree] = io::cells::vtk_to_dolfinx(types[i]);
    // If arbitrary order Lagrange VTK cell (indicated by -1), determine degree
    // from number of nodes
 
    std::uint8_t cell_degree
        = degree == -1 ? io::cells::cell_degree(cell_type, num_nodes)
                       : (std::uint8_t)degree;
    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 (VTK type, degree) to coordinate element index.
  std::map<std::array<std::uint8_t, 2>, std::int32_t> type_to_index;
  std::vector<mesh::CellType> dolfinx_cell_type;
  for (std::array<std::uint8_t, 2> ct : recv_types)
  {
    mesh::CellType cell_type = std::get<0>(io::cells::vtk_to_dolfinx(ct[0]));
    type_to_index.insert({ct, dolfinx_cell_type.size()});
    dolfinx_cell_type.push_back(cell_type);
  }
 
  std::vector npoints = hdf5::read_dataset<std::int64_t>(
      h5file, "/VTKHDF/NumberOfPoints", {0, 1}, true);
  logger().info("Mesh with {} points", npoints[0]);
  std::array<std::int64_t, 2> local_point_range
      = common::local_range(rank, npoints[0], mpi_size);
 
  std::vector<std::int64_t> x_shape
      = hdf5::get_dataset_shape(h5file, "/VTKHDF/Points");
  std::vector<U> points_local = hdf5::read_dataset<U>(h5file, "/VTKHDF/Points",
                                                      local_point_range, true);
 
  // Remove coordinates if 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);
  }
 
  std::vector<std::int64_t> topology = hdf5::read_dataset<std::int64_t>(
      h5file, "/VTKHDF/Connectivity", {offsets.front(), offsets.back()}, true);
  std::ranges::transform(offsets, offsets.begin(),
                         [offset = offsets.front()](std::int64_t x)
                         { return x - offset; });
  h5file.close();
 
  // Create cell topologies for each cell type.
  std::vector<std::vector<std::int64_t>> cells_local(recv_types.size());
  {
    std::vector<std::size_t> num_nodes_per_type(recv_types.size(), 0);
    for (std::size_t j = 0; j < types.size(); ++j)
    {
      std::int32_t type_index = type_to_index.at({types[j], cell_degrees[j]});
      num_nodes_per_type[type_index] += offsets[j + 1] - offsets[j];
    }
    for (std::size_t t = 0; t < cells_local.size(); ++t)
      cells_local[t].reserve(num_nodes_per_type[t]);
  }
  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::ranges::transform(recv_types, std::back_inserter(coordinate_elements),
                         impl::coordinate_element<U>);
 
  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},
                           graph::Partitioner{}, mesh::GhostMode::none,
                           max_facet_to_cell_links, num_threads);
}
} // namespace dolfinx::io::VTKHDF