Adding read_meshtags_from_legacy_h5

src/adios4dolfinx/__init__.py

@@ -18,7 +18,7 @@
     write_mesh,
     write_meshtags,
 )
-from .legacy_readers import read_function_from_legacy_h5, read_mesh_from_legacy_h5
+from .legacy_readers import read_function_from_legacy_h5, read_mesh_from_legacy_h5, read_meshtags_from_legacy_h5
 from .original_checkpoint import write_function_on_input_mesh, write_mesh_input_order
 from .snapshot import snapshot_checkpoint
 

src/adios4dolfinx/legacy_readers.py

@@ -24,7 +24,7 @@
     read_array,
     resolve_adios_scope,
 )
-from .comm_helpers import send_dofs_and_recv_values
+from .comm_helpers import send_dofs_and_recv_values, compute_local_range
 from .utils import (
     compute_dofmap_pos,
     compute_insert_position,
@@ -36,6 +36,7 @@
 
 __all__ = [
     "read_mesh_from_legacy_h5",
+    "read_meshtags_from_legacy_h5",
     "read_function_from_legacy_h5",
 ]
 
@@ -153,6 +154,7 @@ def read_dofmap_legacy(
         global_dofs = np.zeros_like(cells, dtype=np.int64)
         input_cell_positions = cells - local_cell_range[0]
         read_pos = in_offsets[input_cell_positions].astype(np.int32) + dof_pos - in_offsets[0]
+        read_pos = read_pos.astype(np.int32)
         global_dofs = mapped_dofmap[read_pos]
         del input_cell_positions, read_pos
 
@@ -325,6 +327,69 @@ def read_mesh_from_legacy_h5(
     return dolfinx.mesh.create_mesh(MPI.COMM_WORLD, mesh_topology, mesh_geometry, domain)
 
 
+def read_meshtags_from_legacy_h5(
+    filename: pathlib.Path,
+    mesh: dolfinx.mesh.Mesh,
+    dim: int,
+    group: str,
+) -> dolfinx.mesh.MeshTags:
+    """
+    Read MeshFuncion from a `h5`-file generated by legacy DOLFIN `HDF5File.write`
+    or `XDMF.write_checkpoint`.
+
+    Args:
+        filename : Path to `h5` or `xdmf` file
+        mesh : The mesh associated with the MeshFunction / MeshTags
+        dim : Dimension of the MeshFunction / MeshTags
+        group : Group within the `h5` file where the function is stored
+    """
+
+    # Make sure we use the HDF5File and check that the file is present
+    filename = pathlib.Path(filename).with_suffix(".h5")
+    if not filename.is_file():
+        raise FileNotFoundError(f"File {filename} does not exist")
+
+    # Create ADIOS2 reader
+    adios = adios2.ADIOS(mesh.comm)
+    with ADIOSFile(
+        adios=adios,
+        filename=filename,
+        mode=adios2.Mode.Read,
+        io_name="MeshFunction reader",
+        engine="HDF5",
+    ) as adios_file:
+
+        # Get topology
+        if f"{group}/topology" not in adios_file.io.AvailableVariables().keys():
+            raise KeyError(f"MeshFunction topology not found at '{group}/topology'")
+
+        in_topology = adios_file.io.InquireVariable(f"{group}/topology")
+        shape = in_topology.Shape()
+        local_cell_range = compute_local_range(mesh.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"),)
+        adios_file.file.Get(in_topology, topology, adios2.Mode.Sync)
+
+        # Get values
+        if f"{group}/values" not in adios_file.io.AvailableVariables().keys():
+            raise KeyError(f"MeshFunction values not found at '{group}/values'")
+
+        in_values = adios_file.io.InquireVariable(f"{group}/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]], [local_cell_range[1] - local_cell_range[0]]])
+        adios_file.file.Get(in_values, values, adios2.Mode.Sync)
+        # Convert to int32
+        values = values.astype(np.int32)
+
+        # Create meshtags
+        local_entities, local_values = dolfinx.io.utils.distribute_entity_data(mesh, dim, topology, values)
+        adj = dolfinx.cpp.graph.AdjacencyList_int32(local_entities)
+
+        mt = dolfinx.mesh.meshtags_from_entities(mesh, dim, adj, local_values.astype(np.int32, copy=False))
+        return mt
+
+
 def read_function_from_legacy_h5(
     filename: pathlib.Path,
     comm: MPI.Intracomm,