From 89739ce8121dbc0b326875c2191cf3a8723aa199 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= <dokken@simula.no>
Date: Wed, 14 May 2025 13:55:46 +0000
Subject: [PATCH 01/12] Refactor gmshio to use loops and 'correct' names for
 codimensional entities. Ref M. Scroggs paper on permutations, Table 1

---
 python/demo/demo_gmsh.py    |  16 ++--
 python/dolfinx/io/gmshio.py | 181 ++++++++++++------------------------
 2 files changed, 65 insertions(+), 132 deletions(-)

diff --git a/python/demo/demo_gmsh.py b/python/demo/demo_gmsh.py
index 1f23887418c..f6630c56b47 100644
--- a/python/demo/demo_gmsh.py
+++ b/python/demo/demo_gmsh.py
@@ -167,10 +167,10 @@ def create_mesh(comm: MPI.Comm, model: gmsh.model, name: str, filename: str, mod
         mesh_data.cell_tags.name = f"{name}_cells"
     if mesh_data.facet_tags is not None:
         mesh_data.facet_tags.name = f"{name}_facets"
-    if mesh_data.edge_tags is not None:
-        mesh_data.edge_tags.name = f"{name}_edges"
-    if mesh_data.vertex_tags is not None:
-        mesh_data.vertex_tags.name = f"{name}_vertices"
+    if mesh_data.ridge_tags is not None:
+        mesh_data.ridge_tags.name = f"{name}_edges"
+    if mesh_data.peak_tags is not None:
+        mesh_data.peak_tags.name = f"{name}_vertices"
     with XDMFFile(mesh_data.mesh.comm, filename, mode) as file:
         mesh_data.mesh.topology.create_connectivity(2, 3)
         mesh_data.mesh.topology.create_connectivity(1, 3)
@@ -188,15 +188,15 @@ def create_mesh(comm: MPI.Comm, model: gmsh.model, name: str, filename: str, mod
                 mesh_data.mesh.geometry,
                 geometry_xpath=f"/Xdmf/Domain/Grid[@Name='{name}']/Geometry",
             )
-        if mesh_data.edge_tags is not None:
+        if mesh_data.ridge_tags is not None:
             file.write_meshtags(
-                mesh_data.edge_tags,
+                mesh_data.ridge_tags,
                 mesh_data.mesh.geometry,
                 geometry_xpath=f"/Xdmf/Domain/Grid[@Name='{name}']/Geometry",
             )
-        if mesh_data.vertex_tags is not None:
+        if mesh_data.peak_tags is not None:
             file.write_meshtags(
-                mesh_data.vertex_tags,
+                mesh_data.peak_tags,
                 mesh_data.mesh.geometry,
                 geometry_xpath=f"/Xdmf/Domain/Grid[@Name='{name}']/Geometry",
             )
diff --git a/python/dolfinx/io/gmshio.py b/python/dolfinx/io/gmshio.py
index f17e8002102..ba5b3eaace2 100644
--- a/python/dolfinx/io/gmshio.py
+++ b/python/dolfinx/io/gmshio.py
@@ -1,4 +1,4 @@
-# Copyright (C) 2022 Jørgen S. Dokken
+# Copyright (C) 2022-2025 Jørgen S. Dokken, Henrik N. T. Finsberg and Paul T. Kühner
 #
 # This file is part of DOLFINx (https://www.fenicsproject.org)
 #
@@ -78,9 +78,9 @@ class MeshData(typing.NamedTuple):
     Args:
         mesh: Mesh.
         cell_tags: MeshTags for cells.
-        facet_tags: MeshTags for facets.
-        edge_tags: MeshTags for edges.
-        vertex_tags: MeshTags for vertices.
+        facet_tags: MeshTags for facets (codim 1).
+        ridge_tags: MeshTags for edges (codim 2).
+        peak_tags: MeshTags for vertices (codim 3).
         physical_groups: Physical groups in the mesh, where the key
             is the physical name and the value is a tuple with the
             dimension and tag.
@@ -89,8 +89,8 @@ class MeshData(typing.NamedTuple):
     mesh: Mesh
     cell_tags: typing.Optional[MeshTags]
     facet_tags: typing.Optional[MeshTags]
-    edge_tags: typing.Optional[MeshTags]
-    vertex_tags: typing.Optional[MeshTags]
+    ridge_tags: typing.Optional[MeshTags]
+    peak_tags: typing.Optional[MeshTags]
     physical_groups: dict[str, tuple[int, int]]
 
 
@@ -277,8 +277,10 @@ def model_to_mesh(
             distribution of cells across MPI ranks.
 
     Returns:
-        MeshData with mesh, cell tags, facet tags, edge tags,
-        vertex tags and physical groups.
+        MeshData with mesh and tags of corresponding entities by codimension.
+        Codimension 0 is the cell tags, codimension 1 is the facet tags,
+        codimension 2 is the ridge tags and codimension 3 is the peak tags
+        as well as a lookup table from the physical groups by name to integer.
 
     Note:
         For performance, this function should only be called once for
@@ -311,139 +313,70 @@ def model_to_mesh(
         num_nodes = cell_information[perm_sort[-1]]["num_nodes"]
         cell_id, num_nodes = comm.bcast([cell_id, num_nodes], root=rank)
 
-        # Check for facet, edge and vertex data and broadcast relevant info if True
-        has_facet_data = (tdim - 1) in cell_dimensions
-        has_edge_data = (tdim - 2) in cell_dimensions
-        has_vertex_data = (tdim - 3) in cell_dimensions
-
-        has_facet_data = comm.bcast(has_facet_data, root=rank)
-        if has_facet_data:
-            num_facet_nodes = comm.bcast(cell_information[perm_sort[-2]]["num_nodes"], root=rank)
-            gmsh_facet_id = cell_information[perm_sort[-2]]["id"]
-            marked_facets = np.asarray(topologies[gmsh_facet_id]["topology"], dtype=np.int64)
-            facet_values = np.asarray(topologies[gmsh_facet_id]["cell_data"], dtype=np.int32)
-
-        has_edge_data = comm.bcast(has_edge_data, root=rank)
-        if has_edge_data:
-            num_edge_nodes = comm.bcast(cell_information[perm_sort[-3]]["num_nodes"], root=rank)
-            gmsh_edge_id = cell_information[perm_sort[-3]]["id"]
-            marked_edges = np.asarray(topologies[gmsh_edge_id]["topology"], dtype=np.int64)
-            edge_values = np.asarray(topologies[gmsh_edge_id]["cell_data"], dtype=np.int32)
-
-        has_vertex_data = comm.bcast(has_vertex_data, root=rank)
-        if has_vertex_data:
-            num_vertex_nodes = comm.bcast(cell_information[perm_sort[-4]]["num_nodes"], root=rank)
-            gmsh_vertex_id = cell_information[perm_sort[-4]]["id"]
-            marked_vertices = np.asarray(topologies[gmsh_vertex_id]["topology"], dtype=np.int64)
-            vertex_values = np.asarray(topologies[gmsh_vertex_id]["cell_data"], dtype=np.int32)
-
-        cells = np.asarray(topologies[cell_id]["topology"], dtype=np.int64)
-        cell_values = np.asarray(topologies[cell_id]["cell_data"], dtype=np.int32)
         physical_groups = comm.bcast(physical_groups, root=rank)
     else:
-        cell_id, num_nodes = comm.bcast([None, None], root=rank)
-        cells, x = np.empty([0, num_nodes], dtype=np.int32), np.empty([0, gdim], dtype=dtype)
-        cell_values = np.empty((0,), dtype=np.int32)
-
-        has_facet_data = comm.bcast(None, root=rank)
-        if has_facet_data:
-            num_facet_nodes = comm.bcast(None, root=rank)
-            marked_facets = np.empty((0, num_facet_nodes), dtype=np.int32)
-            facet_values = np.empty((0,), dtype=np.int32)
-
-        has_edge_data = comm.bcast(None, root=rank)
-        if has_edge_data:
-            num_edge_nodes = comm.bcast(None, root=rank)
-            marked_edges = np.empty((0, num_edge_nodes), dtype=np.int32)
-            edge_values = np.empty((0,), dtype=np.int32)
-
-        has_vertex_data = comm.bcast(None, root=rank)
-        if has_vertex_data:
-            num_vertex_nodes = comm.bcast(None, root=rank)
-            marked_vertices = np.empty((0, num_vertex_nodes), dtype=np.int32)
-            vertex_values = np.empty((0,), dtype=np.int32)
-
         physical_groups = comm.bcast(None, root=rank)
 
+    # Check for facet, edge and vertex data and broadcast relevant info if True
+    meshtags: dict[int, tuple[bool, npt.NDArray[np.int32], npt.NDArray[np.int32]]] = {}
+    for codim in [0, 1, 2, 3]:
+        # Check if we have any tagged entities on root process and distribute
+        if comm.bcast(tdim - codim in cell_dimensions, root=rank):
+            if comm.rank == rank:
+                position = -1 - codim
+                gmsh_entity_id = cell_information[perm_sort[position]]["id"]
+                marked_entities = np.asarray(topologies[gmsh_entity_id]["topology"], dtype=np.int64)
+                entity_values = np.asarray(topologies[gmsh_entity_id]["cell_data"], dtype=np.int32)
+                num_entity_nodes = comm.bcast(
+                    cell_information[perm_sort[position]]["num_nodes"], root=rank
+                )
+                meshtags[codim] = (marked_entities, entity_values)
+            else:
+                num_entity_nodes = comm.bcast(None, root=rank)
+                marked_entities = np.empty((0, num_entity_nodes), dtype=np.int32)
+                entity_values = np.empty((0,), dtype=np.int32)
+                meshtags[codim] = (marked_entities, entity_values)
+
     # Create distributed mesh
     ufl_domain = ufl_mesh(cell_id, gdim, dtype=dtype)
     gmsh_cell_perm = cell_perm_array(_cpp.mesh.to_type(str(ufl_domain.ufl_cell())), num_nodes)
-    cells = cells[:, gmsh_cell_perm].copy()
+    cells = meshtags[0][0][:, gmsh_cell_perm].copy()
     mesh = create_mesh(comm, cells, x[:, :gdim].astype(dtype, copy=False), ufl_domain, partitioner)
 
-    # Create MeshTags for cells
-    local_entities, local_values = distribute_entity_data(
-        mesh, mesh.topology.dim, cells, cell_values
-    )
-    mesh.topology.create_connectivity(mesh.topology.dim, 0)
-    adj = adjacencylist(local_entities)
-    ct = meshtags_from_entities(
-        mesh, mesh.topology.dim, adj, local_values.astype(np.int32, copy=False)
-    )
-    ct.name = "Cell tags"
-
+    codim_to_name = {0: "cell", 1: "facet", 2: "ridge", 3: "peak"}
+    dolfinx_meshtags: dict[int, MeshTags] = {}
     # Create MeshTags for facets
     topology = mesh.topology
     tdim = topology.dim
-    if has_facet_data:
-        # Permute facets from MSH to DOLFINx ordering
-        # FIXME: This does not work for prism meshes
-        if topology.cell_type == CellType.prism or topology.cell_type == CellType.pyramid:
-            raise RuntimeError(f"Unsupported cell type {topology.cell_type}")
-
-        facet_type = _cpp.mesh.cell_entity_type(
-            _cpp.mesh.to_type(str(ufl_domain.ufl_cell())), tdim - 1, 0
-        )
-        gmsh_facet_perm = cell_perm_array(facet_type, num_facet_nodes)
-        marked_facets = marked_facets[:, gmsh_facet_perm]
-
-        local_entities, local_values = distribute_entity_data(
-            mesh, tdim - 1, marked_facets, facet_values
-        )
-        mesh.topology.create_connectivity(topology.dim - 1, tdim)
-        adj = adjacencylist(local_entities)
-        ft = meshtags_from_entities(mesh, tdim - 1, adj, local_values.astype(np.int32, copy=False))
-        ft.name = "Facet tags"
-    else:
-        ft = None
-
-    if has_edge_data:
-        # Permute edges from MSH to DOLFINx ordering
-        edge_type = _cpp.mesh.cell_entity_type(
-            _cpp.mesh.to_type(str(ufl_domain.ufl_cell())), tdim - 2, 0
-        )
-        gmsh_edge_perm = cell_perm_array(edge_type, num_edge_nodes)
-        marked_edges = marked_edges[:, gmsh_edge_perm]
-
+    for codim in [0, 1, 2, 3]:
+        key = f"{codim_to_name[codim]}_tags"
+        if (
+            (codim == 1
+            and topology.cell_type == CellType.prism)
+            or topology.cell_type == CellType.pyramid
+        ):
+            raise RuntimeError(f"Unsupported facet tag for type {topology.cell_type}")
+
+        meshtag_data = meshtags.get(codim, None)
+        if meshtag_data is None:
+            dolfinx_meshtags[key] = None
+            continue
+
+        (marked_entities, entity_values) = meshtag_data
+
+        # Create MeshTags for cells
         local_entities, local_values = distribute_entity_data(
-            mesh, tdim - 2, marked_edges, edge_values
+            mesh, tdim - codim, marked_entities, entity_values
         )
-        mesh.topology.create_connectivity(topology.dim - 2, tdim)
+        mesh.topology.create_connectivity(tdim - codim, 0)
         adj = adjacencylist(local_entities)
-        et = meshtags_from_entities(mesh, tdim - 2, adj, local_values.astype(np.int32, copy=False))
-        et.name = "Edge tags"
-    else:
-        et = None
-
-    if has_vertex_data:
-        # Permute vertices from MSH to DOLFINx ordering
-        vertex_type = _cpp.mesh.cell_entity_type(
-            _cpp.mesh.to_type(str(ufl_domain.ufl_cell())), tdim - 3, 0
-        )
-        gmsh_vertex_perm = cell_perm_array(vertex_type, num_vertex_nodes)
-        marked_vertices = marked_vertices[:, gmsh_vertex_perm]
-
-        local_entities, local_values = distribute_entity_data(
-            mesh, tdim - 3, marked_vertices, vertex_values
+        et = meshtags_from_entities(
+            mesh, tdim - codim, adj, local_values.astype(np.int32, copy=False)
         )
-        mesh.topology.create_connectivity(topology.dim - 3, tdim)
-        adj = adjacencylist(local_entities)
-        vt = meshtags_from_entities(mesh, tdim - 3, adj, local_values.astype(np.int32, copy=False))
-        vt.name = "Vertex tags"
-    else:
-        vt = None
+        et.name = key
+        dolfinx_meshtags[key] = et
 
-    return MeshData(mesh, ct, ft, et, vt, physical_groups)
+    return MeshData(mesh, **dolfinx_meshtags, physical_groups=physical_groups)
 
 
 def read_from_msh(

From f37f36f5ea592e048a334a6763613f626b6cf0d8 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= <dokken@simula.no>
Date: Wed, 14 May 2025 13:57:48 +0000
Subject: [PATCH 02/12] Ruff format

---
 python/dolfinx/io/gmshio.py | 6 ++----
 1 file changed, 2 insertions(+), 4 deletions(-)

diff --git a/python/dolfinx/io/gmshio.py b/python/dolfinx/io/gmshio.py
index ba5b3eaace2..518b96be7f7 100644
--- a/python/dolfinx/io/gmshio.py
+++ b/python/dolfinx/io/gmshio.py
@@ -351,10 +351,8 @@ def model_to_mesh(
     for codim in [0, 1, 2, 3]:
         key = f"{codim_to_name[codim]}_tags"
         if (
-            (codim == 1
-            and topology.cell_type == CellType.prism)
-            or topology.cell_type == CellType.pyramid
-        ):
+            codim == 1 and topology.cell_type == CellType.prism
+        ) or topology.cell_type == CellType.pyramid:
             raise RuntimeError(f"Unsupported facet tag for type {topology.cell_type}")
 
         meshtag_data = meshtags.get(codim, None)

From 2af54fdeef7e0d8ca9901792ea8c5930565ff4b8 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= <dokken@simula.no>
Date: Wed, 14 May 2025 14:01:48 +0000
Subject: [PATCH 03/12] Come on mypy!

---
 python/dolfinx/io/gmshio.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/python/dolfinx/io/gmshio.py b/python/dolfinx/io/gmshio.py
index 518b96be7f7..f4991e311d1 100644
--- a/python/dolfinx/io/gmshio.py
+++ b/python/dolfinx/io/gmshio.py
@@ -318,7 +318,7 @@ def model_to_mesh(
         physical_groups = comm.bcast(None, root=rank)
 
     # Check for facet, edge and vertex data and broadcast relevant info if True
-    meshtags: dict[int, tuple[bool, npt.NDArray[np.int32], npt.NDArray[np.int32]]] = {}
+    meshtags: dict[int, tuple[npt.NDArray[np.int32], npt.NDArray[np.int32]]] = {}
     for codim in [0, 1, 2, 3]:
         # Check if we have any tagged entities on root process and distribute
         if comm.bcast(tdim - codim in cell_dimensions, root=rank):
@@ -344,7 +344,7 @@ def model_to_mesh(
     mesh = create_mesh(comm, cells, x[:, :gdim].astype(dtype, copy=False), ufl_domain, partitioner)
 
     codim_to_name = {0: "cell", 1: "facet", 2: "ridge", 3: "peak"}
-    dolfinx_meshtags: dict[int, MeshTags] = {}
+    dolfinx_meshtags: dict[str, MeshTags] = {}
     # Create MeshTags for facets
     topology = mesh.topology
     tdim = topology.dim
@@ -374,7 +374,7 @@ def model_to_mesh(
         et.name = key
         dolfinx_meshtags[key] = et
 
-    return MeshData(mesh, **dolfinx_meshtags, physical_groups=physical_groups)
+    return MeshData(mesh, physical_groups=physical_groups, **dolfinx_meshtags)
 
 
 def read_from_msh(

From df1035f3dc40b8dc3174f326bddb225735c2c662 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= <dokken@simula.no>
Date: Wed, 14 May 2025 14:03:04 +0000
Subject: [PATCH 04/12] Optional

---
 python/dolfinx/io/gmshio.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/python/dolfinx/io/gmshio.py b/python/dolfinx/io/gmshio.py
index f4991e311d1..9025f23d975 100644
--- a/python/dolfinx/io/gmshio.py
+++ b/python/dolfinx/io/gmshio.py
@@ -344,7 +344,7 @@ def model_to_mesh(
     mesh = create_mesh(comm, cells, x[:, :gdim].astype(dtype, copy=False), ufl_domain, partitioner)
 
     codim_to_name = {0: "cell", 1: "facet", 2: "ridge", 3: "peak"}
-    dolfinx_meshtags: dict[str, MeshTags] = {}
+    dolfinx_meshtags: dict[str, typing.Optional[MeshTags]] = {}
     # Create MeshTags for facets
     topology = mesh.topology
     tdim = topology.dim

From b10e6f2e88933ecbf08be8e248732fa28c964d9c Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= <dokken@simula.no>
Date: Wed, 14 May 2025 14:22:37 +0000
Subject: [PATCH 05/12] Add missing create connectivity

---
 python/demo/demo_axis.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/python/demo/demo_axis.py b/python/demo/demo_axis.py
index f1a8385fb8f..e0f76db6d7e 100644
--- a/python/demo/demo_axis.py
+++ b/python/demo/demo_axis.py
@@ -590,6 +590,7 @@ def create_eps_mu(pml, rho, eps_bkg, mu_bkg):
 # Marker functions for the scattering efficiency integral
 marker = fem.Function(D)
 scatt_facets = mesh_data.facet_tags.find(scatt_tag)
+mesh_data.mesh.topology.create_connectivity(tdim - 1, tdim)
 incident_cells = mesh.compute_incident_entities(
     mesh_data.mesh.topology, scatt_facets, tdim - 1, tdim
 )

From b56653a6807b53d4dbe70c0fd780a7985016dc2c Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?J=C3=B8rgen=20Schartum=20Dokken?= <dokken92@gmail.com>
Date: Wed, 14 May 2025 21:50:13 +0200
Subject: [PATCH 06/12] Change connectivity creation

---
 python/dolfinx/io/gmshio.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/python/dolfinx/io/gmshio.py b/python/dolfinx/io/gmshio.py
index 9025f23d975..f388f01e3d8 100644
--- a/python/dolfinx/io/gmshio.py
+++ b/python/dolfinx/io/gmshio.py
@@ -366,7 +366,7 @@ def model_to_mesh(
         local_entities, local_values = distribute_entity_data(
             mesh, tdim - codim, marked_entities, entity_values
         )
-        mesh.topology.create_connectivity(tdim - codim, 0)
+        mesh.topology.create_connectivity(tdim - codim, tdim)
         adj = adjacencylist(local_entities)
         et = meshtags_from_entities(
             mesh, tdim - codim, adj, local_values.astype(np.int32, copy=False)

From 03b183eabd85be4d36867947627aed87f0555da7 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= <dokken@simula.no>
Date: Wed, 14 May 2025 20:22:50 +0000
Subject: [PATCH 07/12] Create demo that tags peaks, ridges, facets and cells
 in a grid

---
 python/demo/demo_gmsh.py    | 19 +++++++++++++------
 python/dolfinx/io/gmshio.py |  8 +++++++-
 2 files changed, 20 insertions(+), 7 deletions(-)

diff --git a/python/demo/demo_gmsh.py b/python/demo/demo_gmsh.py
index f6630c56b47..2c4bc829a0d 100644
--- a/python/demo/demo_gmsh.py
+++ b/python/demo/demo_gmsh.py
@@ -36,7 +36,8 @@
 
 # +
 def gmsh_sphere(model: gmsh.model, name: str) -> gmsh.model:
-    """Create a Gmsh model of a sphere.
+    """Create a Gmsh model of a sphere and tag sub entitites
+    from all co-dimensions (peaks, ridges, facets and cells).
 
     Args:
         model: Gmsh model to add the mesh to.
@@ -53,9 +54,15 @@ def gmsh_sphere(model: gmsh.model, name: str) -> gmsh.model:
     # Synchronize OpenCascade representation with gmsh model
     model.occ.synchronize()
 
-    # Add physical marker for cells. It is important to call this
-    # function after OpenCascade synchronization
-    model.add_physical_group(dim=3, tags=[sphere])
+    # Add physical tag for sphere
+    model.add_physical_group(dim=3, tags=[sphere], tag=1)
+
+    # Embed all sub-entities from the GMSH model into the sphere and tag them
+    for dim in [0, 1, 2]:
+        entities = model.getEntities(dim)
+        entity_ids = [entity[1] for entity in entities]
+        model.mesh.embed(dim, entity_ids, 3, sphere)
+        model.add_physical_group(dim=dim, tags=entity_ids, tag=dim)
 
     # Generate the mesh
     model.mesh.generate(dim=3)
@@ -168,9 +175,9 @@ def create_mesh(comm: MPI.Comm, model: gmsh.model, name: str, filename: str, mod
     if mesh_data.facet_tags is not None:
         mesh_data.facet_tags.name = f"{name}_facets"
     if mesh_data.ridge_tags is not None:
-        mesh_data.ridge_tags.name = f"{name}_edges"
+        mesh_data.ridge_tags.name = f"{name}_ridges"
     if mesh_data.peak_tags is not None:
-        mesh_data.peak_tags.name = f"{name}_vertices"
+        mesh_data.peak_tags.name = f"{name}_peaks"
     with XDMFFile(mesh_data.mesh.comm, filename, mode) as file:
         mesh_data.mesh.topology.create_connectivity(2, 3)
         mesh_data.mesh.topology.create_connectivity(1, 3)
diff --git a/python/dolfinx/io/gmshio.py b/python/dolfinx/io/gmshio.py
index f388f01e3d8..bfcd973e2f4 100644
--- a/python/dolfinx/io/gmshio.py
+++ b/python/dolfinx/io/gmshio.py
@@ -184,7 +184,13 @@ def extract_topology_and_markers(
             # one cell-type,
             # i.e. facets of prisms and pyramid meshes are not supported
             (entity_types, entity_tags, entity_topologies) = model.mesh.getElements(dim, tag=entity)
-            assert len(entity_types) == 1
+
+            if len(entity_types) == 2:
+                raise RuntimeError(
+                    f"Unsupported mesh with multiple cell types {entity_types} for entity {entity}"
+                )
+            elif len(entity_types) == 0:
+                continue
 
             # Determine number of local nodes per element to create the
             # topology of the elements

From 61fa6dbbb3e07d1b6640f63266ad42fc998d5ea2 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= <dokken@simula.no>
Date: Wed, 14 May 2025 21:18:19 +0000
Subject: [PATCH 08/12] Add documentation

---
 python/dolfinx/io/gmshio.py | 91 ++++++++++++++++++++-----------------
 1 file changed, 50 insertions(+), 41 deletions(-)

diff --git a/python/dolfinx/io/gmshio.py b/python/dolfinx/io/gmshio.py
index bfcd973e2f4..6e41436977d 100644
--- a/python/dolfinx/io/gmshio.py
+++ b/python/dolfinx/io/gmshio.py
@@ -300,60 +300,61 @@ def model_to_mesh(
         x = extract_geometry(model)
         topologies, physical_groups = extract_topology_and_markers(model)
 
-        # Extract Gmsh cell id, dimension of cell and number of nodes to
-        # cell for each
+        # Extract Gmsh cell id, dimension of cell and number of nodes in cell
         num_cell_types = len(topologies.keys())
-        cell_information = dict()
-        cell_dimensions = np.zeros(num_cell_types, dtype=np.int32)
+        element_ids = np.zeros(num_cell_types, dtype=np.int32)
+        element_dimensions = np.zeros(num_cell_types, dtype=np.int32)
+        num_nodes_per_element = np.zeros(num_cell_types, dtype=np.int32)
         for i, element in enumerate(topologies.keys()):
             _, dim, _, num_nodes, _, _ = model.mesh.getElementProperties(element)
-            cell_information[i] = {"id": element, "dim": dim, "num_nodes": num_nodes}
-            cell_dimensions[i] = dim
+            element_ids[i] = element
+            element_dimensions[i] = dim
+            num_nodes_per_element[i] = num_nodes
 
-        # Sort elements by ascending dimension
-        perm_sort = np.argsort(cell_dimensions)
-
-        # Broadcast cell type data and geometric dimension
-        cell_id = cell_information[perm_sort[-1]]["id"]
-        tdim = cell_information[perm_sort[-1]]["dim"]
-        num_nodes = cell_information[perm_sort[-1]]["num_nodes"]
-        cell_id, num_nodes = comm.bcast([cell_id, num_nodes], root=rank)
-
-        physical_groups = comm.bcast(physical_groups, root=rank)
+        # Broadcast information to all other ranks
+        element_dimensions, element_ids, num_nodes_per_element = comm.bcast(
+            (element_dimensions, element_ids, num_nodes_per_element), root=rank
+        )
     else:
-        physical_groups = comm.bcast(None, root=rank)
+        element_dimensions, element_ids, num_nodes_per_element = comm.bcast(
+            (None, None, None), root=rank
+        )
+
+    # Sort elements by descending dimension
+    assert len(np.unique(element_dimensions)) == len(element_dimensions)
+    perm_sort = np.argsort(element_dimensions)[::-1]
+    tdim = int(element_dimensions[perm_sort[0]])
 
-    # Check for facet, edge and vertex data and broadcast relevant info if True
     meshtags: dict[int, tuple[npt.NDArray[np.int32], npt.NDArray[np.int32]]] = {}
-    for codim in [0, 1, 2, 3]:
-        # Check if we have any tagged entities on root process and distribute
-        if comm.bcast(tdim - codim in cell_dimensions, root=rank):
-            if comm.rank == rank:
-                position = -1 - codim
-                gmsh_entity_id = cell_information[perm_sort[position]]["id"]
-                marked_entities = np.asarray(topologies[gmsh_entity_id]["topology"], dtype=np.int64)
-                entity_values = np.asarray(topologies[gmsh_entity_id]["cell_data"], dtype=np.int32)
-                num_entity_nodes = comm.bcast(
-                    cell_information[perm_sort[position]]["num_nodes"], root=rank
-                )
-                meshtags[codim] = (marked_entities, entity_values)
-            else:
-                num_entity_nodes = comm.bcast(None, root=rank)
-                marked_entities = np.empty((0, num_entity_nodes), dtype=np.int32)
-                entity_values = np.empty((0,), dtype=np.int32)
-                meshtags[codim] = (marked_entities, entity_values)
+    for position in perm_sort:
+        # Extract entity data for tagged codim entities
+        codim = tdim - element_dimensions[position]
+        if comm.rank == rank:
+            gmsh_entity_id = element_ids[position]
+            marked_entities = np.asarray(topologies[gmsh_entity_id]["topology"], dtype=np.int64)
+            entity_values = np.asarray(topologies[gmsh_entity_id]["cell_data"], dtype=np.int32)
+            meshtags[codim] = (marked_entities, entity_values)
+        else:
+            marked_entities = np.empty((0, num_nodes_per_element[position]), dtype=np.int32)
+            entity_values = np.empty((0,), dtype=np.int32)
+            meshtags[codim] = (marked_entities, entity_values)
 
     # Create distributed mesh
-    ufl_domain = ufl_mesh(cell_id, gdim, dtype=dtype)
+    ufl_domain = ufl_mesh(element_ids[perm_sort[0]], gdim, dtype=dtype)
+    num_nodes = num_nodes_per_element[perm_sort[0]]
     gmsh_cell_perm = cell_perm_array(_cpp.mesh.to_type(str(ufl_domain.ufl_cell())), num_nodes)
     cells = meshtags[0][0][:, gmsh_cell_perm].copy()
+    if comm.rank != rank:
+        x = np.empty([0, gdim], dtype=dtype)  # No nodes on other than root rank
     mesh = create_mesh(comm, cells, x[:, :gdim].astype(dtype, copy=False), ufl_domain, partitioner)
+    assert tdim == mesh.topology.dim, (
+        f"{mesh.topology.dim=} does not match Gmsh model dimension {tdim}"
+    )
 
+    # Create MeshTags for entity data extracted from GMSH
+    topology = mesh.topology
     codim_to_name = {0: "cell", 1: "facet", 2: "ridge", 3: "peak"}
     dolfinx_meshtags: dict[str, typing.Optional[MeshTags]] = {}
-    # Create MeshTags for facets
-    topology = mesh.topology
-    tdim = topology.dim
     for codim in [0, 1, 2, 3]:
         key = f"{codim_to_name[codim]}_tags"
         if (
@@ -366,12 +367,14 @@ def model_to_mesh(
             dolfinx_meshtags[key] = None
             continue
 
+        # Distribute entity data [[e0_v0, e0_v1, ...], [e1_v0, e1_v1, ...], ...]
+        # which is made in global input indices to local indices on the
+        # owning process
         (marked_entities, entity_values) = meshtag_data
-
-        # Create MeshTags for cells
         local_entities, local_values = distribute_entity_data(
             mesh, tdim - codim, marked_entities, entity_values
         )
+        # Create MeshTags object from the local entities
         mesh.topology.create_connectivity(tdim - codim, tdim)
         adj = adjacencylist(local_entities)
         et = meshtags_from_entities(
@@ -380,6 +383,12 @@ def model_to_mesh(
         et.name = key
         dolfinx_meshtags[key] = et
 
+    # Broadcast physical groups (string to integer mapping) to all ranks
+    if comm.rank == rank:
+        physical_groups = comm.bcast(physical_groups, root=rank)
+    else:
+        physical_groups = comm.bcast(None, root=rank)
+
     return MeshData(mesh, physical_groups=physical_groups, **dolfinx_meshtags)
 
 

From c7b84f7ee5f50d1ee3c37debd74de1e74b2b8a04 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?J=C3=B8rgen=20Schartum=20Dokken?= <dokken92@gmail.com>
Date: Wed, 14 May 2025 23:20:24 +0200
Subject: [PATCH 09/12] Revert new connectivity

---
 python/demo/demo_axis.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/python/demo/demo_axis.py b/python/demo/demo_axis.py
index e0f76db6d7e..f1a8385fb8f 100644
--- a/python/demo/demo_axis.py
+++ b/python/demo/demo_axis.py
@@ -590,7 +590,6 @@ def create_eps_mu(pml, rho, eps_bkg, mu_bkg):
 # Marker functions for the scattering efficiency integral
 marker = fem.Function(D)
 scatt_facets = mesh_data.facet_tags.find(scatt_tag)
-mesh_data.mesh.topology.create_connectivity(tdim - 1, tdim)
 incident_cells = mesh.compute_incident_entities(
     mesh_data.mesh.topology, scatt_facets, tdim - 1, tdim
 )

From 13cca1f91661038f83f7e0fea069ef0f2f3d0e81 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?J=C3=B8rgen=20Schartum=20Dokken?= <dokken92@gmail.com>
Date: Fri, 16 May 2025 13:30:44 +0200
Subject: [PATCH 10/12] Generalize check

---
 python/dolfinx/io/gmshio.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/python/dolfinx/io/gmshio.py b/python/dolfinx/io/gmshio.py
index 6e41436977d..ca211101186 100644
--- a/python/dolfinx/io/gmshio.py
+++ b/python/dolfinx/io/gmshio.py
@@ -185,7 +185,7 @@ def extract_topology_and_markers(
             # i.e. facets of prisms and pyramid meshes are not supported
             (entity_types, entity_tags, entity_topologies) = model.mesh.getElements(dim, tag=entity)
 
-            if len(entity_types) == 2:
+            if len(entity_types) > 1:
                 raise RuntimeError(
                     f"Unsupported mesh with multiple cell types {entity_types} for entity {entity}"
                 )

From 047b76975cffa05fc168e997e501d4b520e34165 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?J=C3=B8rgen=20Schartum=20Dokken?= <dokken92@gmail.com>
Date: Fri, 16 May 2025 13:31:13 +0200
Subject: [PATCH 11/12] Renaming as suggested by M

Co-authored-by: Michal Habera <michal.habera@gmail.com>
---
 python/dolfinx/io/gmshio.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/python/dolfinx/io/gmshio.py b/python/dolfinx/io/gmshio.py
index ca211101186..15bdacb5017 100644
--- a/python/dolfinx/io/gmshio.py
+++ b/python/dolfinx/io/gmshio.py
@@ -79,8 +79,8 @@ class MeshData(typing.NamedTuple):
         mesh: Mesh.
         cell_tags: MeshTags for cells.
         facet_tags: MeshTags for facets (codim 1).
-        ridge_tags: MeshTags for edges (codim 2).
-        peak_tags: MeshTags for vertices (codim 3).
+        ridge_tags: MeshTags for ridges (codim 2).
+        peak_tags: MeshTags for peaks (codim 3).
         physical_groups: Physical groups in the mesh, where the key
             is the physical name and the value is a tuple with the
             dimension and tag.

From 0edee6b539784f0a165d6d97fa3112c35d7e0cfd Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= <dokken@simula.no>
Date: Thu, 22 May 2025 14:12:44 +0000
Subject: [PATCH 12/12] Improve documentation and naming

---
 python/dolfinx/io/gmshio.py | 55 ++++++++++++++++++++++---------------
 1 file changed, 33 insertions(+), 22 deletions(-)

diff --git a/python/dolfinx/io/gmshio.py b/python/dolfinx/io/gmshio.py
index 15bdacb5017..1e78309895b 100644
--- a/python/dolfinx/io/gmshio.py
+++ b/python/dolfinx/io/gmshio.py
@@ -300,58 +300,69 @@ def model_to_mesh(
         x = extract_geometry(model)
         topologies, physical_groups = extract_topology_and_markers(model)
 
-        # Extract Gmsh cell id, dimension of cell and number of nodes in cell
-        num_cell_types = len(topologies.keys())
-        element_ids = np.zeros(num_cell_types, dtype=np.int32)
-        element_dimensions = np.zeros(num_cell_types, dtype=np.int32)
-        num_nodes_per_element = np.zeros(num_cell_types, dtype=np.int32)
+        # Extract Gmsh entity (cell) id, topological dimension and number of nodes
+        # which is used to create an appropriate coordinate element, and seperate
+        # higher topological entities from lower topological entities (e.g. facets,
+        # ridges and peaks).
+        num_unique_entities = len(topologies.keys())
+        element_ids = np.zeros(num_unique_entities, dtype=np.int32)
+        entity_tdim = np.zeros(num_unique_entities, dtype=np.int32)
+        num_nodes_per_element = np.zeros(num_unique_entities, dtype=np.int32)
         for i, element in enumerate(topologies.keys()):
             _, dim, _, num_nodes, _, _ = model.mesh.getElementProperties(element)
             element_ids[i] = element
-            element_dimensions[i] = dim
+            entity_tdim[i] = dim
             num_nodes_per_element[i] = num_nodes
 
         # Broadcast information to all other ranks
-        element_dimensions, element_ids, num_nodes_per_element = comm.bcast(
-            (element_dimensions, element_ids, num_nodes_per_element), root=rank
+        entity_tdim, element_ids, num_nodes_per_element = comm.bcast(
+            (entity_tdim, element_ids, num_nodes_per_element), root=rank
         )
     else:
-        element_dimensions, element_ids, num_nodes_per_element = comm.bcast(
-            (None, None, None), root=rank
-        )
+        entity_tdim, element_ids, num_nodes_per_element = comm.bcast((None, None, None), root=rank)
 
     # Sort elements by descending dimension
-    assert len(np.unique(element_dimensions)) == len(element_dimensions)
-    perm_sort = np.argsort(element_dimensions)[::-1]
-    tdim = int(element_dimensions[perm_sort[0]])
+    assert len(np.unique(entity_tdim)) == len(entity_tdim)
+    perm_sort = np.argsort(entity_tdim)[::-1]
+
+    # Extract position of the highest topological entity and its topological dimension
+    cell_position = perm_sort[0]
+    tdim = int(entity_tdim[cell_position])
 
+    # Extract entity -> node connectivity for all cells and sub-entities marked in the GMSH model
     meshtags: dict[int, tuple[npt.NDArray[np.int32], npt.NDArray[np.int32]]] = {}
     for position in perm_sort:
-        # Extract entity data for tagged codim entities
-        codim = tdim - element_dimensions[position]
+        codim = tdim - entity_tdim[position]
         if comm.rank == rank:
             gmsh_entity_id = element_ids[position]
             marked_entities = np.asarray(topologies[gmsh_entity_id]["topology"], dtype=np.int64)
             entity_values = np.asarray(topologies[gmsh_entity_id]["cell_data"], dtype=np.int32)
             meshtags[codim] = (marked_entities, entity_values)
         else:
+            # Any other process than input rank does not have any entities
             marked_entities = np.empty((0, num_nodes_per_element[position]), dtype=np.int32)
             entity_values = np.empty((0,), dtype=np.int32)
             meshtags[codim] = (marked_entities, entity_values)
 
-    # Create distributed mesh
-    ufl_domain = ufl_mesh(element_ids[perm_sort[0]], gdim, dtype=dtype)
-    num_nodes = num_nodes_per_element[perm_sort[0]]
+    # Create a UFL Mesh object for the GMSH element with the highest topoligcal dimension
+    ufl_domain = ufl_mesh(element_ids[cell_position], gdim, dtype=dtype)
+
+    # Get cell->node connectivity and  permute to FEniCS ordering
+    num_nodes = num_nodes_per_element[cell_position]
     gmsh_cell_perm = cell_perm_array(_cpp.mesh.to_type(str(ufl_domain.ufl_cell())), num_nodes)
-    cells = meshtags[0][0][:, gmsh_cell_perm].copy()
+    cell_connectivity = meshtags[0][0][:, gmsh_cell_perm].copy()
+
+    # Create a distributed mesh, where mesh nodes are only destributed from the input rank
     if comm.rank != rank:
         x = np.empty([0, gdim], dtype=dtype)  # No nodes on other than root rank
-    mesh = create_mesh(comm, cells, x[:, :gdim].astype(dtype, copy=False), ufl_domain, partitioner)
+    mesh = create_mesh(
+        comm, cell_connectivity, x[:, :gdim].astype(dtype, copy=False), ufl_domain, partitioner
+    )
     assert tdim == mesh.topology.dim, (
         f"{mesh.topology.dim=} does not match Gmsh model dimension {tdim}"
     )
 
-    # Create MeshTags for entity data extracted from GMSH
+    # Create MeshTags for all sub entities
     topology = mesh.topology
     codim_to_name = {0: "cell", 1: "facet", 2: "ridge", 3: "peak"}
     dolfinx_meshtags: dict[str, typing.Optional[MeshTags]] = {}