The gmshio does currently not support prism elements, as there is no clear API in DOLFINx for reading facet tags.
Here is a modified version of the gmshio that can read the mesh and cell tags:
import gmsh, mpi4py
try:
from dolfinx.io import gmshio
except ImportError:
from dolfinx.io import gmsh as gmshio
import numpy as np
from dolfinx import default_real_type
import dolfinx
import numpy.typing as npt
import ufl
import basix.ufl
gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 0) # disable output messages
# Map from Gmsh cell type identifier (integer) to DOLFINx cell type and
# degree https://gmsh.info//doc/texinfo/gmsh.html#MSH-file-format
_gmsh_to_cells = {
1: ("interval", 1),
2: ("triangle", 1),
3: ("quadrilateral", 1),
4: ("tetrahedron", 1),
5: ("hexahedron", 1),
6: ("prism", 1),
8: ("interval", 2),
9: ("triangle", 2),
10: ("quadrilateral", 2),
11: ("tetrahedron", 2),
12: ("hexahedron", 2),
15: ("point", 0),
21: ("triangle", 3),
26: ("interval", 3),
29: ("tetrahedron", 3),
36: ("quadrilateral", 3),
92: ("hexahedron", 3),
}
def ufl_mesh(gmsh_cell: int, gdim: int, dtype: npt.DTypeLike) -> ufl.Mesh:
"""Create a UFL mesh from a Gmsh cell identifier and geometric
dimension.
See https://gmsh.info//doc/texinfo/gmsh.html#MSH-file-format.
Args:
gmsh_cell: Gmsh cell identifier.
gdim: Geometric dimension of the mesh.
Returns:
UFL Mesh using Lagrange elements (equispaced) of the
corresponding DOLFINx cell.
"""
try:
shape, degree = _gmsh_to_cells[gmsh_cell]
except KeyError as e:
print(f"Unknown cell type {gmsh_cell}.")
raise e
cell = ufl.Cell(shape)
element = basix.ufl.element(
basix.ElementFamily.P,
cell.cellname(),
degree,
basix.LagrangeVariant.equispaced,
shape=(gdim,),
dtype=dtype, # type: ignore[arg-type]
)
return ufl.Mesh(element)
def model_to_mesh(
model,
comm,
rank: int,
gdim: int = 3,
partitioner = None,
dtype=default_real_type,
) -> tuple[dolfinx.mesh.Mesh, dolfinx.mesh.MeshTags]:
"""Create a Mesh from a Gmsh model.
Creates a :class:`dolfinx.mesh.Mesh` from the physical entities of
the highest topological dimension in the Gmsh model. In parallel,
the gmsh model is processed on one MPI rank, and the
:class:`dolfinx.mesh.Mesh` is distributed across ranks.
Args:
model: Gmsh model.
comm: MPI communicator to use for mesh creation.
rank: MPI rank that the Gmsh model is initialized on.
gdim: Geometrical dimension of the mesh.
partitioner: Function that computes the parallel
distribution of cells across MPI ranks.
Returns:
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
large problems. For reuse, it is recommended to save the mesh
and corresponding tags using :class:`dolfinx.io.XDMFFile` after
creation for efficient access.
"""
valid_mesh = None
if comm.rank == rank:
assert model is not None, "Gmsh model is None on rank responsible for mesh creation."
# Get mesh geometry and mesh topology for each element
x = gmshio.extract_geometry(model)
try:
topologies, physical_groups = gmshio.extract_topology_and_markers(model)
valid_mesh = len(physical_groups) > 0
except ValueError as e:
topologies = gmshio.extract_topology_and_markers(model)
# 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
entity_tdim[i] = dim
num_nodes_per_element[i] = num_nodes
# Broadcast information to all other ranks
entity_tdim, element_ids, num_nodes_per_element = comm.bcast(
(entity_tdim, element_ids, num_nodes_per_element), root=rank
)
else:
entity_tdim, element_ids, num_nodes_per_element = comm.bcast((None, None, None), root=rank)
valid_mesh = comm.bcast(valid_mesh, root=rank)
if not valid_mesh:
raise RuntimeError("No 'physical groups' in gmsh mesh. Cannot continue.")
# Sort elements by descending dimension
perm_sort = np.argsort(entity_tdim)[::-1]
if len(perm_sort) > 1:
assert entity_tdim[perm_sort[0]] != entity_tdim[perm_sort[1]]
# Extract position of the highest topological entity and its
# topological dimension
cell_position = perm_sort[0]
tdim = int(entity_tdim[cell_position])
# Check that all cells are tagged once
error_msg = ""
if comm.rank == rank:
_elementTypes, _elementTags, _nodeTags = model.mesh.getElements(dim=tdim, tag=-1)
_elementType_dim = _elementTypes[0]
if _elementType_dim not in topologies.keys():
error_msg = "All cells are expected to be tagged once; none found"
num_cells_tagged = len(topologies[_elementType_dim]["entity_tags"])
if (num_cells := len(_elementTags[0])) != num_cells_tagged:
error_msg = (
"All cells are expected to be tagged once;"
+ f"found: {num_cells_tagged}, expected: {num_cells}"
)
num_cells_tagged_once = len(np.unique(topologies[_elementType_dim]["entity_tags"]))
if num_cells_tagged != num_cells_tagged_once:
error_msg = "All cells are expected to be tagged once, found duplicates"
error_msg = comm.bcast(error_msg, root=rank)
if error_msg != "":
raise RuntimeError(error_msg)
# Extract entity -> node connectivity for all cells and sub-entities
# marked in the GMSH model
meshtags: dict[int, tuple[npt.NDArray[np.int64], npt.NDArray[np.int32]]] = {}
for position in perm_sort:
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 a UFL Mesh object for the GMSH element with the highest
# topological 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 = gmshio.cell_perm_array(dolfinx.cpp.mesh.to_type(str(ufl_domain.ufl_cell())), num_nodes)
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 = dolfinx.mesh.create_mesh(
comm, cell_connectivity, ufl_domain, x[:, :gdim].astype(dtype, copy=False), partitioner
)
assert tdim == mesh.topology.dim, (
f"{mesh.topology.dim=} does not match Gmsh model dimension {tdim}"
)
# Create MeshTags for all sub entities
topology = mesh.topology
codim_to_name = {0: "cell", 1: "facet", 2: "ridge", 3: "peak"}
dolfinx_meshtags: dict[str, dolfinx.mesh.MeshTags | None] = {}
for codim in [0]:
key = f"{codim_to_name[codim]}_tags"
if (
codim == 1 and topology.cell_type == dolfinx.mesh.CellType.prism
) or topology.cell_type == dolfinx.mesh.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
# 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
local_entities, local_values = dolfinx.io.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 = dolfinx.graph.adjacencylist(local_entities)
et = dolfinx.mesh.meshtags_from_entities(
mesh, tdim - codim, adj, local_values.astype(np.int32, copy=False)
)
et.name = key
return mesh, et
occ = gmsh.model.occ
gdim = 2
rect1 = occ.addRectangle(0, 0, 0, 1, 1)
extruded = occ.extrude(
[(gdim, rect1)],
0, 0, 1, # extrusion vector
numElements=[1], # only one element in extrusion direction
recombine=True
)
occ.synchronize()
# add_sub physical tags
all_doms = occ.getEntities(gdim+1)
for j, dom in enumerate(all_doms):
gmsh.model.addPhysicalGroup(dom[0], [dom[1]], j + 1) # create the main group/node
# number all boundaries
all_edges = gmsh.model.getEntities(gdim)
for j, edge in enumerate(all_edges):
gmsh.model.addPhysicalGroup(edge[0], [edge[1]], edge[1]) # create the main group/node
gmsh.model.mesh.generate(gdim+1)
mpi_rank = 0
mesh, ct = model_to_mesh(gmsh.model, mpi4py.MPI.COMM_WORLD, mpi_rank, gdim+1)
with dolfinx.io.XDMFFile(mpi4py.MPI.COMM_WORLD, "mwe.xdmf", "w") as xdmf:
xdmf.write_mesh(mesh)
xdmf.write_meshtags(ct, mesh.geometry)
