Support for tagged domains from legacy h5 files

[Reidmen]

Question:

Is it planned to allow reading meshes that include marked boundaries?

I would guess something along the lines:

def read_mesh_from_legacy_checkpoint(
    filename: str, cell_type: str = "tetrahedron"
) -> Tuple[dolfinx.mesh.Mesh, list[NDArray[np.float64]]]:
    """
    Returns: 
        mesh
        list of array of boundary values
    """
    pass

Context: My group's workflow uses tagged meshes (coming from a segmentation procedure from MRI images). We are in the process of updating all our codebase to dolfinx, but still many meshes are in the legacy h5 format.

[jorgensd]

Are the meshes and markers in h5-format or xdmf format? If they are in xdmffiles They can already be read by dolfinx.

[Reidmen]

Sadly, the whole pipeline ends up with meshes that are in h5 format. But thanks for the info about xdmf! I'll look in more detail!

[finsberg]

The following code can be used to read a facet tags from a tagged mesh (also works in parallel)

import adios4dolfinx
import adios2
import numpy as np

comm = MPI.COMM_WORLD

# Mesh file in XDMF
with dolfinx.io.XDMFFile(comm, mesh_file, "r") as xdmf:
    mesh = xdmf.read_mesh(name="mesh")

mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim-1)

adios = adios2.ADIOS(mesh.comm)
io_name = "ReadMeshFunction"
io = adios.DeclareIO(io_name)
io.SetEngine("HDF5")
# Read the corresponding H5File
file = io.Open(str(mesh_file.with_suffix(".h5")), adios2.Mode.Read)
file.BeginStep()
in_topology = io.InquireVariable("/MeshFunction/0/mesh/topology")
shape = in_topology.Shape()
local_cell_range = adios4dolfinx.comm_helpers.compute_local_range(comm, shape[0])

in_topology.SetSelection(
    [[local_cell_range[0], 0], [local_cell_range[1] - local_cell_range[0], shape[1]]]
)
topology = np.empty(
    (local_cell_range[1] - local_cell_range[0], shape[1]),
    dtype=in_topology.Type().strip("_t"),
)
file.Get(in_topology, topology, adios2.Mode.Sync)

in_values = io.InquireVariable("/MeshFunction/0/values")        
values = np.empty(local_cell_range[1] - local_cell_range[0], dtype=in_values.Type().strip("_t"))
in_values.SetSelection(
    [[local_cell_range[0], 0], [local_cell_range[1] - local_cell_range[0], 1]]
)

file.Get(in_values, values, adios2.Mode.Sync)

file.EndStep()
file.Close()
adios.RemoveIO(io_name)

# Convert to int32
values = values.astype(np.int32)
local_entities, local_values = dolfinx.io.utils.distribute_entity_data(mesh, mesh.topology.dim - 1, topology, values)
adj = dolfinx.cpp.graph.AdjacencyList_int32(local_entities)

ct = dolfinx.mesh.meshtags_from_entities(mesh, mesh.topology.dim - 1, adj, local_values.astype(np.int32, copy=False))
ct.name = "Facet tags"

# Save for visualization in paraview
mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim)
with dolfinx.io.XDMFFile(mesh.comm, "ffun.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_meshtags(ct, mesh.geometry)

[Reidmen]

Thanks a lot @finsberg !