Coverage for /dolfinx-env/lib/python3.14/site-packages/io4dolfinx/backends/gmsh/backend.py: 91%
379 statements
« prev ^ index » next coverage.py v7.17.0, created at 2026-10-23 13:26 +0000
1"""
2Gmsh interface to io4dolfinx.
4Reads and writes Gmsh ``.msh`` files (MSH 2.2 and 4.1, ASCII and binary) through the
5Gmsh Python API, reusing :mod:`dolfinx.io.gmsh` for the geometry, the topology and the
6node ordering. All file access happens on rank 0. Only the element types of
7``dolfinx.io.gmsh`` are supported, and DOLFINx 0.11 or later is required.
9In DOLFINx terms:
11- As in :func:`dolfinx.io.gmsh.model_to_mesh`, the cells are the elements of the
12 highest dimension in a physical group, and every cell is in exactly one physical
13 group. Node ``i`` has node tag ``i + 1``. An MSH file holds a single mesh, so meshes
14 are not time dependent.
15- Point and cell data are post-processing views (``$NodeData``/``$ElementData``) keyed
16 by node and element tags. The steps of a view are its time series. Complex data is
17 split into the views ``real_{name}`` and ``imag_{name}``, as in
18 :class:`dolfinx.io.XDMFFile`.
19- Mesh tags are physical groups, with the physical tag as the tag value.
21SPDX License identifier: MIT
23Copyright: Jørgen S. Dokken and Simula Research Laboratory
24"""
26from pathlib import Path
27from typing import Any, Callable, TypeVar
29from mpi4py import MPI
31import basix
32import dolfinx
33import fenicsx_compat
34import numpy as np
35import numpy.typing as npt
36from dolfinx.io.gmsh import (
37 _gmsh_to_cells,
38 cell_perm_array,
39 extract_geometry,
40 extract_topology_and_markers,
41)
43try:
44 import gmsh
45except ImportError:
46 raise ModuleNotFoundError("The gmsh-backend requires gmsh")
48from ...structures import ArrayData, FunctionData, MeshData, MeshTagsData, ReadMeshData
49from ...utils import check_file_exists, codim_names
50from .. import FileMode, ReadMode
52read_mode = ReadMode.serial
54# Gmsh options changed by this backend, restored if Gmsh was initialized by the caller
55_options = ("General.Terminal", "Mesh.Binary", "Mesh.SaveAll", "PostProcessing.SaveMesh")
57Local = TypeVar("Local")
58Shared = TypeVar("Shared")
61def _run_on_root(
62 comm: MPI.Comm, func: Callable[[list[int]], tuple[Local, Shared]]
63) -> tuple[Local | None, Shared]:
64 """Run a function in a Gmsh session on rank 0 and share its outcome with all ranks.
66 If Gmsh is already initialized, e.g. by the caller, the function runs in a new model,
67 which is removed afterwards together with the views it created, and the previous
68 model and options are restored. Otherwise Gmsh is initialized and finalized around it.
70 Args:
71 comm: MPI communicator. Collective on all ranks.
72 func: Function to run on rank 0. It takes the tags of the views that existed
73 before the session, which do not belong to the file, and returns a tuple of
74 data kept on rank 0 and data broadcast to all ranks.
76 Returns:
77 Tuple containing:
78 - The first output of ``func`` on rank 0, ``None`` on other ranks
79 - The second output of ``func`` on all ranks
81 Raises:
82 NotImplementedError: If DOLFINx is older than 0.11.
83 Exception: An exception raised by ``func`` is raised on all ranks.
84 """
85 if fenicsx_compat.before("0.11.0.dev0"):
86 raise NotImplementedError("The gmsh backend requires DOLFINx 0.11 or later.")
88 result: Local | None = None
89 shared: Shared | None = None
90 error: Exception | None = None
91 if comm.rank == 0:
92 initialized = gmsh.isInitialized()
93 if initialized:
94 previous_model = gmsh.model.getCurrent()
95 previous_options = {option: gmsh.option.getNumber(option) for option in _options}
96 # Views are global in gmsh, not per model, so existing ones must be skipped and kept
97 preexisting_views = list(gmsh.view.getTags())
98 gmsh.model.add("io4dolfinx")
99 else:
100 gmsh.initialize(interruptible=False)
101 preexisting_views = []
102 try:
103 gmsh.option.setNumber("General.Terminal", 0)
104 result, shared = func(preexisting_views)
105 except Exception as e:
106 error = e
107 finally:
108 if initialized:
109 for view in gmsh.view.getTags():
110 if view not in preexisting_views:
111 gmsh.view.remove(view)
112 gmsh.model.remove()
113 gmsh.model.setCurrent(previous_model)
114 for option, value in previous_options.items():
115 gmsh.option.setNumber(option, value)
116 else:
117 gmsh.finalize()
118 error, shared = comm.bcast((error, shared), root=0)
119 if error is not None:
120 raise error
121 return result, shared # type: ignore[return-value]
124def _find_view(
125 name: str, preexisting_views: list[int]
126) -> tuple[list[int], npt.NDArray[np.float64]] | None:
127 """Find the views of the current Gmsh session holding the data with a given name.
129 A gmsh view holds all time steps of one field, each step storing its time.
130 Complex data is stored as two views, ``real_{name}`` and ``imag_{name}``, as
131 :class:`dolfinx.io.XDMFFile` does.
133 Args:
134 name: Name of the data
135 preexisting_views: Tags of the views that do not belong to the file, which are skipped
137 Returns:
138 ``None`` if there is no such data, otherwise a tuple containing:
139 - The view tags, ``[view]`` for real data and ``[real, imag]`` for complex data
140 - The time of each step of the views
141 """
142 views = {
143 gmsh.view.option.getString(view, "Name"): view
144 for view in gmsh.view.getTags()
145 if view not in preexisting_views
146 }
147 real, imag = f"real_{name}", f"imag_{name}"
148 if name in views:
149 found = [views[name]]
150 elif real in views and imag in views:
151 found = [views[real], views[imag]]
152 elif real in views or imag in views:
153 raise RuntimeError(
154 f"Found only one part of complex data '{name}', expected '{real}' and '{imag}'."
155 )
156 else:
157 return None
158 num_steps = int(gmsh.view.option.getNumber(found[0], "NbTimeStep"))
159 times = [gmsh.view.getModelData(found[0], step)[3] for step in range(num_steps)]
160 return found, np.array(times, dtype=np.float64)
163def _tagged_entities() -> dict[
164 int,
165 tuple[npt.NDArray[np.int64], npt.NDArray[np.int64], npt.NDArray[np.int32], str, int] | None,
166]:
167 """Get the elements in physical groups of the current Gmsh model, by dimension.
169 Uses :func:`dolfinx.io.gmsh.extract_topology_and_markers`, so as in
170 :func:`dolfinx.io.gmsh.model_to_mesh` the highest dimension holds the cells, and an
171 element in several physical groups appears once per group.
173 Returns:
174 Map from topological dimension to ``None`` if the elements of that dimension are
175 of several types or of a type not supported by :mod:`dolfinx.io.gmsh`, otherwise
176 to a tuple containing:
177 - The element tags, sorted
178 - The element nodes as zero-based node indices in DOLFINx ordering, of shape
179 ``(num_elements, num_nodes_per_element)``, one row per element tag
180 - The physical tag of each element
181 - The DOLFINx cell type of the elements
182 - The degree of the elements
183 """
184 topologies, _ = extract_topology_and_markers(gmsh.model)
185 if len(topologies) == 0:
186 raise ValueError("No physical groups in the Gmsh model.")
187 types_by_dim: dict[int, list[int]] = {}
188 for gmsh_type in topologies:
189 if gmsh_type in _gmsh_to_cells:
190 cell_type = dolfinx.mesh.to_type(_gmsh_to_cells[gmsh_type][0])
191 dim = dolfinx.mesh.cell_dim(cell_type)
192 else:
193 dim = gmsh.model.mesh.getElementProperties(gmsh_type)[1]
194 types_by_dim.setdefault(dim, []).append(gmsh_type)
196 entities: dict[int, Any] = {}
197 for dim, types in types_by_dim.items():
198 if len(types) > 1 or types[0] not in _gmsh_to_cells:
199 entities[dim] = None
200 continue
201 cell_name, degree = _gmsh_to_cells[types[0]]
202 data = topologies[types[0]]
203 tags = data["entity_tags"].astype(np.int64)
204 order = np.argsort(tags, kind="stable")
205 nodes = data["topology"].astype(np.int64)
206 perm = cell_perm_array(dolfinx.mesh.to_type(cell_name), nodes.shape[1])
207 values = data["cell_data"].astype(np.int32)
208 entities[dim] = (tags[order], nodes[order][:, perm], values[order], cell_name, degree)
209 return entities
212def get_default_backend_args(arguments: dict[str, Any] | None) -> dict[str, Any]:
213 """Get default backend arguments given a set of input arguments.
215 Args:
216 arguments: Input backend arguments
218 Returns:
219 The arguments, with defaults for:
220 - ``"gdim"``: Geometric dimension of a mesh read from file, ``3`` as in
221 :func:`dolfinx.io.gmsh.read_from_msh`
222 - ``"dtype"``: Floating point type of the coordinates and data read from file,
223 ``np.float64``
224 - ``"binary"``: Write binary instead of ASCII MSH. ``None`` keeps the encoding
225 of an existing file and writes a new one as ASCII.
226 """
227 args: dict[str, Any] = {"gdim": 3, "dtype": np.float64, "binary": None}
228 args.update(arguments or {})
229 return args
232def write_attributes(
233 filename: Path | str,
234 comm: MPI.Comm,
235 name: str,
236 attributes: dict[str, np.ndarray],
237 backend_args: dict[str, Any] | None,
238):
239 """Write attributes to file.
241 Args:
242 filename: Path to file to write to
243 comm: MPI communicator used in storage
244 name: Name of the attribute group
245 attributes: Dictionary of attributes to write
246 backend_args: Arguments to backend
247 """
248 raise NotImplementedError("The gmsh backend cannot write attributes.")
251def read_attributes(
252 filename: Path | str,
253 comm: MPI.Comm,
254 name: str,
255 backend_args: dict[str, Any] | None,
256) -> dict[str, Any]:
257 """Read attributes from file.
259 Args:
260 filename: Path to file to read from
261 comm: MPI communicator used in storage
262 name: Name of the attribute group
263 backend_args: Arguments to backend
265 Returns:
266 Dictionary of attributes read from file
267 """
268 raise NotImplementedError("The gmsh backend cannot read attributes.")
271def read_timestamps(
272 filename: Path | str,
273 comm: MPI.Comm,
274 function_name: str,
275 backend_args: dict[str, Any] | None,
276) -> npt.NDArray[np.float64 | str]: # type: ignore[type-var]
277 """Read the times of the steps of point or cell data.
279 Args:
280 filename: Path to file to read from
281 comm: MPI communicator used in storage
282 function_name: Name of the point or cell data
283 backend_args: Arguments to backend
285 Returns:
286 The time of each step, on all ranks
287 """
288 check_file_exists(filename)
290 def read(preexisting_views: list[int]) -> tuple[None, npt.NDArray[np.float64]]:
291 gmsh.merge(str(filename))
292 if (found := _find_view(function_name, preexisting_views)) is None:
293 raise KeyError(f"Could not find data '{function_name}' in {filename}.")
294 return None, found[1]
296 return _run_on_root(comm, read)[1]
299def read_function_names(
300 filename: Path | str, comm: MPI.Comm, backend_args: dict[str, Any] | None
301) -> list[str]:
302 """Read the names of all point and cell data in a file.
304 Args:
305 filename: Path to file
306 comm: MPI communicator to launch IO on.
307 backend_args: Arguments to backend
309 Returns:
310 The names of the views in the file, on all ranks. A complex pair of views
311 ``real_{name}`` and ``imag_{name}`` is listed as ``name``.
312 """
313 check_file_exists(filename)
315 def read(preexisting_views: list[int]) -> tuple[None, list[str]]:
316 gmsh.merge(str(filename))
317 views = [view for view in gmsh.view.getTags() if view not in preexisting_views]
318 names = [gmsh.view.option.getString(view, "Name") for view in views]
319 complex_names = {n[5:] for n in names if n.startswith("real_") and f"imag_{n[5:]}" in names}
320 functions = [
321 n[5:] if n.startswith(("real_", "imag_")) and n[5:] in complex_names else n
322 for n in names
323 ]
324 return None, list(dict.fromkeys(functions))
326 return _run_on_root(comm, read)[1]
329def write_mesh(
330 filename: Path | str,
331 comm: MPI.Comm,
332 mesh: MeshData,
333 backend_args: dict[str, Any] | None = None,
334 mode: FileMode = FileMode.write,
335 time: float = 0.0,
336):
337 """Write a mesh to a Gmsh file.
339 Nodes and cells get the tag one above their global index. The cells form a single
340 entity in the physical group with tag 1, the default cell marker, as every cell must
341 be in a physical group (see :func:`dolfinx.io.gmsh.model_to_mesh`). Cell tags written
342 with :func:`write_meshtags` replace it.
344 Args:
345 filename: Path to file to write to, with suffix ``.msh``
346 comm: MPI communicator used in storage
347 mesh: Internal data structure for the mesh data to save to file
348 backend_args: Arguments to backend, ``"binary"`` selects binary MSH
349 mode: Only :attr:`FileMode.write`, as an MSH file holds one mesh
350 time: Ignored, as an MSH file holds one mesh
351 """
352 if mode != FileMode.write:
353 raise NotImplementedError("The gmsh backend writes a single mesh per file.")
354 if Path(filename).suffix != ".msh":
355 raise ValueError(f"Expected a filename with suffix .msh, got {filename}.")
356 if mesh.degree > 2 and mesh.lagrange_variant != int(basix.LagrangeVariant.equispaced):
357 raise NotImplementedError("Gmsh only supports equispaced Lagrange geometries.")
358 cell = (mesh.cell_type, mesh.degree)
359 if (gmsh_type := next((t for t, c in _gmsh_to_cells.items() if c == cell), None)) is None:
360 raise NotImplementedError(
361 f"Gmsh has no element for a {mesh.cell_type} of degree {mesh.degree}."
362 )
363 backend_args = get_default_backend_args(backend_args)
365 # Rank order is global index order
366 nodes = comm.gather(mesh.local_geometry, root=0)
367 cells = comm.gather(mesh.local_topology, root=0)
369 def write(preexisting_views: list[int]) -> tuple[None, None]:
370 assert nodes is not None and cells is not None
371 x = np.vstack(nodes)
372 x_3d = np.zeros((x.shape[0], 3), dtype=np.float64)
373 x_3d[:, : x.shape[1]] = x
374 topology = np.vstack(cells)
375 gmsh_topology = np.empty_like(topology)
376 cell_type = dolfinx.mesh.to_type(mesh.cell_type)
377 gmsh_topology[:, cell_perm_array(cell_type, topology.shape[1])] = topology
379 tdim = dolfinx.mesh.cell_dim(cell_type)
380 gmsh.model.addDiscreteEntity(tdim, 1)
381 gmsh.model.mesh.addNodes(tdim, 1, np.arange(1, x.shape[0] + 1), x_3d.reshape(-1))
382 gmsh.model.mesh.addElementsByType(
383 1,
384 gmsh_type,
385 np.arange(1, topology.shape[0] + 1),
386 gmsh_topology.reshape(-1) + 1,
387 )
388 gmsh.model.addPhysicalGroup(tdim, [1], 1)
389 gmsh.option.setNumber("Mesh.Binary", int(bool(backend_args["binary"])))
390 gmsh.option.setNumber("Mesh.SaveAll", 1)
391 gmsh.write(str(filename))
392 return None, None
394 _run_on_root(comm, write)
397def read_mesh_data(
398 filename: Path | str,
399 comm: MPI.Comm,
400 time: str | float | None = None,
401 read_from_partition: bool = False,
402 backend_args: dict[str, Any] | None = None,
403) -> ReadMeshData:
404 """Read mesh data from a Gmsh file.
406 As in :func:`dolfinx.io.gmsh.model_to_mesh`, the cells are the elements of the
407 highest dimension in a physical group, and every cell must be in exactly one physical
408 group. The cells must be of one type, and node tags contiguous. All data is read on
409 rank 0.
411 Args:
412 filename: Path to file to read from
413 comm: MPI communicator to distribute the mesh over
414 time: Ignored, as an MSH file holds one mesh
415 read_from_partition: Read mesh with partition from file (not supported)
416 backend_args: Arguments to backend, see :func:`get_default_backend_args`
418 Returns:
419 Internal data structure for the mesh data read from file
420 """
421 if read_from_partition:
422 raise NotImplementedError("The gmsh backend cannot read partitioning data.")
423 check_file_exists(filename)
424 backend_args = get_default_backend_args(backend_args)
425 gdim = backend_args["gdim"]
427 def read(
428 preexisting_views: list[int],
429 ) -> tuple[tuple[npt.NDArray[np.float64], npt.NDArray[np.int64]], tuple[str, int, int]]:
430 gmsh.merge(str(filename))
431 x = extract_geometry(gmsh.model)[:, :gdim]
432 tagged = _tagged_entities()
433 tdim = max(tagged)
434 if (cell_data := tagged[tdim]) is None:
435 raise NotImplementedError(f"Mixed or unsupported cell types in {filename}.")
436 tags, cells, _, cell_type, degree = cell_data
437 num_cells = sum(len(t) for t in gmsh.model.mesh.getElements(dim=tdim, tag=-1)[1])
438 if len(tags) != num_cells or len(np.unique(tags)) != len(tags):
439 raise ValueError(
440 f"All cells are expected to be in exactly one physical group, as in "
441 f"dolfinx.io.gmsh.model_to_mesh; found {len(tags)} tags for {num_cells} cells."
442 )
443 return (x, cells), (cell_type, degree, cells.shape[1])
445 local, (cell_type, degree, num_nodes_per_cell) = _run_on_root(comm, read)
446 if local is None:
447 x = np.zeros((0, gdim), dtype=np.float64)
448 cells = np.zeros((0, num_nodes_per_cell), dtype=np.int64)
449 else:
450 x, cells = local
451 return ReadMeshData(
452 cells=cells,
453 cell_type=cell_type,
454 x=x.astype(backend_args["dtype"]),
455 lvar=int(basix.LagrangeVariant.equispaced),
456 degree=degree,
457 )
460def write_meshtags(
461 filename: str | Path,
462 comm: MPI.Comm,
463 data: MeshTagsData,
464 backend_args: dict[str, Any] | None = None,
465):
466 """Write mesh tags to a Gmsh file as physical groups.
468 Each tag value becomes a physical group with the value as tag, named
469 ``f"{data.name}_{value}"``. Gmsh physical tags are positive and unique per
470 dimension, so tag values must be positive.
472 Cell tags replace the marker of the cells they cover, moving them to the physical
473 group of their value, while other cells keep their marker, e.g. the default marker 1
474 of :func:`write_mesh`. Markers no cell has any more are removed. Entities of lower
475 dimension are added as new elements, and only one set of mesh tags can be written
476 per lower dimension.
478 Args:
479 filename: Path to file to write to, holding a mesh written with
480 :func:`write_mesh`
481 comm: MPI communicator used in storage
482 data: Internal data structure for the mesh tags to save to file
483 backend_args: Arguments to backend, ``"binary"`` selects binary MSH
484 """
485 check_file_exists(filename)
486 backend_args = get_default_backend_args(backend_args)
487 indices = comm.gather(data.indices, root=0)
488 values = comm.gather(data.values, root=0)
490 def write(preexisting_views: list[int]) -> tuple[None, None]:
491 assert indices is not None and values is not None
492 entities = np.vstack(indices).astype(np.int64)
493 tag_values = np.hstack(values).astype(np.int32)
494 if np.any(tag_values <= 0):
495 raise ValueError(f"Gmsh physical tags must be positive, got {np.unique(tag_values)}.")
496 if (binary := backend_args["binary"]) is None:
497 # The second line of an MSH file is "version file-type data-size"
498 with open(filename, "rb") as file:
499 file.readline()
500 binary = file.readline().split()[1] == b"1"
502 gmsh.merge(str(filename))
503 tagged = _tagged_entities()
504 tdim = max(tagged)
505 if (cell_data := tagged[tdim]) is None:
506 raise NotImplementedError(f"Mixed or unsupported cell types in {filename}.")
507 cell_tags, cells, _, _, degree = cell_data
508 dim = data.dim
509 groups = [tag for d, tag in gmsh.model.getPhysicalGroups() if d == dim]
510 if dim < tdim and len(groups) > 0:
511 raise NotImplementedError(
512 f"{filename} already has physical groups of dimension {dim}, and only one "
513 "set of mesh tags can be written per dimension below the cells."
514 )
516 assert data.cell_type is not None
517 entity_degree = 0 if data.cell_type == "point" else degree
518 cell = (data.cell_type, entity_degree)
519 if (gmsh_type := next((t for t, c in _gmsh_to_cells.items() if c == cell), None)) is None:
520 raise NotImplementedError(
521 f"Gmsh has no element for a {data.cell_type} of degree {entity_degree}."
522 )
523 gmsh_entities = np.empty_like(entities)
524 perm = cell_perm_array(dolfinx.mesh.to_type(data.cell_type), entities.shape[1])
525 gmsh_entities[:, perm] = entities
527 if dim == tdim:
528 # Find the element tag of each tagged cell by matching node sets
529 _, inverse = np.unique(
530 np.sort(np.vstack([cells, entities]), axis=1), axis=0, return_inverse=True
531 )
532 inverse = inverse.reshape(-1)
533 cell_position = np.empty(inverse.max() + 1, dtype=np.int64)
534 cell_position[inverse[: len(cells)]] = np.arange(len(cells))
535 element_tags = cell_tags[cell_position[inverse[len(cells) :]]]
536 # Gmsh rejects removing elements an entity does not hold
537 for d, entity in gmsh.model.getEntities(tdim):
538 entity_tags = gmsh.model.mesh.getElements(d, entity)[1]
539 held = np.intersect1d(np.hstack([*entity_tags, []]), element_tags)
540 if len(held) > 0:
541 gmsh.model.mesh.removeElements(d, entity, held.astype(np.int64))
542 else:
543 first = gmsh.model.mesh.getMaxElementTag() + 1
544 element_tags = np.arange(first, first + len(entities), dtype=np.int64)
546 for value in np.unique(tag_values):
547 selected = tag_values == value
548 if value in groups:
549 # A marker the cells already have, e.g. the default marker 1
550 entity = gmsh.model.getEntitiesForPhysicalGroup(dim, int(value))[0]
551 gmsh.model.setPhysicalName(dim, int(value), f"{data.name}_{value}")
552 else:
553 entity = gmsh.model.addDiscreteEntity(dim)
554 gmsh.model.addPhysicalGroup(dim, [entity], int(value), f"{data.name}_{value}")
555 gmsh.model.mesh.addElementsByType(
556 entity, gmsh_type, element_tags[selected], gmsh_entities[selected].reshape(-1) + 1
557 )
558 if dim == tdim:
559 # Remove the markers no cell has any more
560 for tag in groups:
561 group_entities = gmsh.model.getEntitiesForPhysicalGroup(dim, tag)
562 num_cells = sum(
563 len(t) for e in group_entities for t in gmsh.model.mesh.getElements(dim, e)[1]
564 )
565 if num_cells == 0:
566 gmsh.model.removePhysicalGroups([(dim, tag)])
568 gmsh.option.setNumber("Mesh.Binary", int(binary))
569 gmsh.option.setNumber("Mesh.SaveAll", 1)
570 gmsh.write(str(filename))
571 # Append the views without repeating the mesh
572 gmsh.option.setNumber("PostProcessing.SaveMesh", 0)
573 for view in gmsh.view.getTags():
574 if view not in preexisting_views:
575 gmsh.view.write(view, str(filename), append=True)
576 return None, None
578 _run_on_root(comm, write)
581def read_meshtags_data(
582 filename: str | Path,
583 comm: MPI.Comm,
584 name: str,
585 backend_args: dict[str, Any] | None = None,
586) -> MeshTagsData:
587 """Read mesh tags from the physical groups of a Gmsh file.
589 Args:
590 filename: Path to file to read from
591 comm: MPI communicator used in storage
592 name: ``"cell"``, ``"facet"``, ``"ridge"`` or ``"peak"`` for all physical groups of
593 that codimension, as in :func:`dolfinx.io.gmsh.model_to_mesh`. Otherwise the
594 physical group called ``name``, or the groups called ``f"{name}_{tag}"``, as
595 written by :func:`write_meshtags`.
596 backend_args: Arguments to backend
598 Returns:
599 Internal data structure for the mesh tags read from file, with the physical tag as
600 value. All data is on rank 0.
601 """
602 check_file_exists(filename)
604 def read(
605 preexisting_views: list[int],
606 ) -> tuple[tuple[npt.NDArray[np.int64], npt.NDArray[np.int32]], tuple[int, int]]:
607 gmsh.merge(str(filename))
608 tagged = _tagged_entities()
609 tdim = max(tagged)
610 groups = gmsh.model.getPhysicalGroups()
611 if name in codim_names:
612 dim = tdim - codim_names[name]
613 group_tags = [tag for d, tag in groups if d == dim]
614 else:
615 matches = [
616 (d, tag)
617 for d, tag in groups
618 if gmsh.model.getPhysicalName(d, tag) in (name, f"{name}_{tag}")
619 ]
620 if len(matches) == 0:
621 available = sorted(gmsh.model.getPhysicalName(d, tag) for d, tag in groups)
622 raise KeyError(
623 f"Could not find physical group '{name}' in {filename}. Available: {available}"
624 )
625 if len(dims := {d for d, _ in matches}) != 1:
626 raise RuntimeError(f"Physical groups '{name}' have several dimensions {dims}.")
627 dim = dims.pop()
628 group_tags = [tag for _, tag in matches]
629 if dim not in tagged:
630 return (np.zeros((0, 0), dtype=np.int64), np.zeros(0, dtype=np.int32)), (dim, 0)
631 if (entity_data := tagged[dim]) is None:
632 raise NotImplementedError(f"Mixed or unsupported element types of dimension {dim}.")
633 _, entities, tag_values, _, _ = entity_data
634 keep = np.isin(tag_values, group_tags)
635 return (entities[keep], tag_values[keep]), (dim, entities.shape[1])
637 local, (dim, num_nodes_per_entity) = _run_on_root(comm, read)
638 if local is None:
639 local = (np.zeros((0, num_nodes_per_entity), dtype=np.int64), np.zeros(0, dtype=np.int32))
640 return MeshTagsData(name=name, values=local[1], indices=local[0], dim=dim)
643def read_dofmap(
644 filename: str | Path,
645 comm: MPI.Comm,
646 name: str,
647 backend_args: dict[str, Any] | None,
648) -> dolfinx.graph.AdjacencyList:
649 """Read the dofmap of a function with a given name.
651 Args:
652 filename: Path to file to read from
653 comm: MPI communicator used in storage
654 name: Name of the function to read the dofmap for
655 backend_args: Arguments to backend
657 Returns:
658 Dofmap as an :class:`dolfinx.graph.AdjacencyList`
659 """
660 raise NotImplementedError("The gmsh backend cannot make checkpoints.")
663def read_dofs(
664 filename: str | Path,
665 comm: MPI.Comm,
666 name: str,
667 time: float,
668 backend_args: dict[str, Any] | None,
669) -> tuple[npt.NDArray[np.float32 | np.float64 | np.complex64 | np.complex128], int]:
670 """Read the dofs (values) of a function with a given name from a given timestep.
672 Args:
673 filename: Path to file to read from
674 comm: MPI communicator used in storage
675 name: Name of the function to read the dofs for
676 time: Time stamp associated with the function to read
677 backend_args: Arguments to backend
679 Returns:
680 Contiguous sequence of degrees of freedom (with respect to input data)
681 and the global starting point on the process.
682 Process 0 has [0, M), process 1 [M, N), process 2 [N, O) etc.
683 """
684 raise NotImplementedError("The gmsh backend cannot make checkpoints.")
687def read_cell_perms(
688 comm: MPI.Comm, filename: Path | str, backend_args: dict[str, Any] | None
689) -> npt.NDArray[np.uint32]:
690 """
691 Read cell permutation from file with given communicator,
692 Split in continuous chunks based on number of cells in the input data.
694 Args:
695 comm: MPI communicator used in storage
696 filename: Path to file to read from
697 backend_args: Arguments to backend
699 Returns:
700 Contiguous sequence of permutations (with respect to input data)
701 Process 0 has [0, M), process 1 [M, N), process 2 [N, O) etc.
702 """
703 raise NotImplementedError("The gmsh backend cannot make checkpoints.")
706def write_function(
707 filename: Path | str,
708 comm: MPI.Comm,
709 u: FunctionData,
710 time: float,
711 mode: FileMode,
712 backend_args: dict[str, Any] | None,
713):
714 """Write a function to file.
716 Args:
717 filename: Path to file to write to
718 comm: MPI communicator used in storage
719 u: Internal data structure for the function data to save to file
720 time: Time stamp associated with function
721 mode: File-mode to store the function
722 backend_args: Arguments to backend
723 """
724 raise NotImplementedError("The gmsh backend cannot make checkpoints.")
727def read_legacy_mesh(
728 filename: Path | str, comm: MPI.Comm, group: str
729) -> tuple[npt.NDArray[np.int64], npt.NDArray[np.floating], str | None]:
730 """Read in the mesh topology, geometry and (optionally) cell type from a
731 legacy DOLFIN HDF5-file.
733 Args:
734 filename: Path to file to read from
735 comm: MPI communicator used in storage
736 group: Group in HDF5 file where mesh is stored
738 Returns:
739 Tuple containing:
740 - Topology as a (num_cells, num_vertices_per_cell) array of global vertex indices
741 - Geometry as a (num_vertices, geometric_dimension) array of vertex coordinates
742 - Cell type as a string (e.g. "tetrahedron") or None if not found
743 """
744 raise NotImplementedError("The gmsh backend cannot read legacy DOLFIN meshes.")
747def snapshot_checkpoint(
748 filename: Path | str,
749 mode: FileMode,
750 u: dolfinx.fem.Function,
751 backend_args: dict[str, Any] | None,
752):
753 """Create a snapshot checkpoint of a dolfinx function.
755 Args:
756 filename: Path to file to read from
757 mode: File-mode to store the function
758 u: dolfinx function to create a snapshot checkpoint for
759 backend_args: Arguments to backend
760 """
761 raise NotImplementedError("The gmsh backend cannot make checkpoints.")
764def read_hdf5_array(
765 comm: MPI.Comm,
766 filename: Path | str,
767 group: str,
768 backend_args: dict[str, Any] | None,
769) -> tuple[np.ndarray, int]:
770 """Read an array from an HDF5 file.
772 Args:
773 comm: MPI communicator used in storage
774 filename: Path to file to read from
775 group: Group in HDF5 file where array is stored
776 backend_args: Arguments to backend
778 Returns:
779 Tuple containing:
780 - Numpy array read from file
781 - Global starting point on the process.
782 Process 0 has [0, M), process 1 [M, N), process 2 [N, O) etc.
783 """
784 raise NotImplementedError("The gmsh backend cannot read HDF5 arrays.")
787def read_point_data(
788 filename: Path | str,
789 name: str,
790 comm: MPI.Comm,
791 time: str | float | None,
792 backend_args: dict[str, Any] | None,
793) -> tuple[np.ndarray, int]:
794 """Read point data (``$NodeData``) from a Gmsh file.
796 Args:
797 filename: Path to file
798 name: Name of point data
799 comm: Communicator to launch IO on.
800 time: Time of the step to read, ``None`` for the first step
801 backend_args: Arguments to backend, ``"dtype"`` sets the type of the data, or of
802 its real and imaginary parts for complex data
804 Returns:
805 Tuple containing:
806 - The data of shape ``(num_nodes, num_components)``, ordered by node tag, on
807 rank 0. Nodes without data get zeros. Complex data, stored as the views
808 ``real_{name}`` and ``imag_{name}``, is combined.
809 - The global start of the data on this process, ``0``
810 """
811 check_file_exists(filename)
812 backend_args = get_default_backend_args(backend_args)
814 def read(
815 preexisting_views: list[int],
816 ) -> tuple[npt.NDArray[np.float64 | np.complex128], tuple[int, bool]]:
817 gmsh.merge(str(filename))
818 if (found := _find_view(name, preexisting_views)) is None:
819 raise KeyError(f"Could not find point data '{name}' in {filename}.")
820 views, times = found
821 if time is None:
822 step = 0
823 else:
824 steps = np.flatnonzero(np.isclose(times, float(time)))
825 if len(steps) != 1:
826 raise KeyError(
827 f"Expected one step of '{name}' at time {time} in {filename}, found "
828 f"{len(steps)}. Available times: {times.tolist()}"
829 )
830 step = int(steps[0])
831 parts = []
832 for view in views:
833 data_type, tags, data, _, num_components = gmsh.view.getModelData(view, step)
834 if data_type != "NodeData":
835 raise ValueError(f"'{name}' in {filename} is {data_type}, not NodeData.")
836 part = np.zeros((gmsh.model.mesh.getMaxNodeTag(), num_components), dtype=np.float64)
837 part[np.asarray(tags, dtype=np.int64) - 1] = np.asarray(data).reshape(
838 len(tags), num_components
839 )
840 parts.append(part)
841 values = parts[0] if len(parts) == 1 else parts[0] + 1j * parts[1]
842 return values, (num_components, len(parts) == 2)
844 values, (num_components, is_complex) = _run_on_root(comm, read)
845 dtype = backend_args["dtype"]
846 if is_complex:
847 # complex64 from float32 parts, complex128 from float64
848 dtype = np.result_type(dtype, np.complex64)
849 if values is None:
850 values = np.zeros((0, num_components), dtype=dtype)
851 return values.astype(dtype), 0
854def read_cell_data(
855 filename: Path | str,
856 name: str,
857 comm: MPI.Comm,
858 time: str | float | None,
859 backend_args: dict[str, Any] | None,
860) -> tuple[npt.NDArray[np.int64], np.ndarray]:
861 """Read cell data (``$ElementData``) from a Gmsh file.
863 Args:
864 filename: Path to file
865 name: Name of cell data
866 comm: Communicator to launch IO on.
867 time: Time of the step to read, ``None`` for the first step
868 backend_args: Arguments to backend, ``"dtype"`` sets the type of the data, or of
869 its real and imaginary parts for complex data
871 Returns:
872 Tuple containing:
873 - The nodes of the cells with data, as zero-based node indices in DOLFINx
874 ordering, on rank 0
875 - The data of these cells, of shape ``(num_cells, num_components)``, on rank 0.
876 Complex data, stored as the views ``real_{name}`` and ``imag_{name}``, is
877 combined.
878 """
879 check_file_exists(filename)
880 backend_args = get_default_backend_args(backend_args)
882 def read(
883 preexisting_views: list[int],
884 ) -> tuple[
885 tuple[npt.NDArray[np.int64], npt.NDArray[np.float64 | np.complex128]],
886 tuple[int, int, bool],
887 ]:
888 gmsh.merge(str(filename))
889 if (found := _find_view(name, preexisting_views)) is None:
890 raise KeyError(f"Could not find cell data '{name}' in {filename}.")
891 views, times = found
892 if time is None:
893 step = 0
894 else:
895 steps = np.flatnonzero(np.isclose(times, float(time)))
896 if len(steps) != 1:
897 raise KeyError(
898 f"Expected one step of '{name}' at time {time} in {filename}, found "
899 f"{len(steps)}. Available times: {times.tolist()}"
900 )
901 step = int(steps[0])
902 tagged = _tagged_entities()
903 if (cell_data := tagged[max(tagged)]) is None:
904 raise NotImplementedError(f"Mixed or unsupported cell types in {filename}.")
905 element_tags, cells = cell_data[:2]
906 parts, part_tags = [], []
907 for view in views:
908 data_type, tags, data, _, num_components = gmsh.view.getModelData(view, step)
909 if data_type != "ElementData":
910 raise ValueError(f"'{name}' in {filename} is {data_type}, not ElementData.")
911 part_tags.append(np.asarray(tags, dtype=np.int64))
912 parts.append(np.asarray(data, dtype=np.float64).reshape(len(tags), num_components))
913 if len(parts) == 2 and not np.array_equal(part_tags[0], part_tags[1]):
914 raise ValueError(f"The real and imaginary parts of '{name}' are on different cells.")
915 rows = np.searchsorted(element_tags, part_tags[0])
916 if np.any(rows >= len(element_tags)) or np.any(element_tags[rows] != part_tags[0]):
917 raise ValueError(f"'{name}' in {filename} holds data on elements that are not cells.")
918 values = parts[0] if len(parts) == 1 else parts[0] + 1j * parts[1]
919 return (cells[rows], values), (num_components, cells.shape[1], len(parts) == 2)
921 local, (num_components, num_nodes_per_cell, is_complex) = _run_on_root(comm, read)
922 dtype = backend_args["dtype"]
923 if is_complex:
924 # complex64 from float32 parts, complex128 from float64
925 dtype = np.result_type(dtype, np.complex64)
926 if local is None:
927 local = (
928 np.zeros((0, num_nodes_per_cell), dtype=np.int64),
929 np.zeros((0, num_components), dtype=dtype),
930 )
931 return local[0], local[1].astype(dtype)
934def write_data(
935 filename: Path | str,
936 array_data: ArrayData,
937 comm: MPI.Comm,
938 time: str | float | None,
939 mode: FileMode,
940 backend_args: dict[str, Any] | None,
941):
942 """Write point or cell data to a Gmsh file, as a step of the view named after it.
944 Point data is written as ``$NodeData`` and cell data as ``$ElementData``, keyed by
945 global index + 1. Complex data is split into the views ``real_{name}`` and
946 ``imag_{name}``, as :class:`dolfinx.io.XDMFFile` does. The mesh must have been
947 written to the file with :func:`write_mesh`, from the same mesh.
949 Args:
950 filename: Path to file
951 array_data: Data to write to file
952 comm: MPI communicator to open the file with
953 time: Time of the step, ``None`` for ``0``
954 mode: :attr:`FileMode.append` adds the step, or replaces the step at the same
955 time; :attr:`FileMode.write` replaces all steps of the data
956 backend_args: Arguments to backend, ``"binary"`` selects binary MSH
957 """
958 if (data_type := {"Point": "NodeData", "Cell": "ElementData"}.get(array_data.type)) is None:
959 raise ValueError(f"Unknown data type {array_data.type}, expected 'Point' or 'Cell'.")
960 check_file_exists(filename)
961 backend_args = get_default_backend_args(backend_args)
962 step_time = 0.0 if time is None else float(time)
963 values = comm.gather(array_data.values, root=0)
965 def write(preexisting_views: list[int]) -> tuple[None, None]:
966 assert values is not None
967 if (binary := backend_args["binary"]) is None:
968 # The second line of an MSH file is "version file-type data-size"
969 with open(filename, "rb") as file:
970 file.readline()
971 binary = file.readline().split()[1] == b"1"
973 gmsh.merge(str(filename))
974 data = np.vstack(values)
975 name = array_data.name
976 if np.iscomplexobj(data):
977 parts = [(f"real_{name}", data.real), (f"imag_{name}", data.imag)]
978 else:
979 parts = [(name, data)]
980 for part_name, part in parts:
981 found = _find_view(part_name, preexisting_views)
982 if found is not None and len(found[0]) > 1:
983 raise ValueError(f"'{part_name}' in {filename} is complex data.")
984 if found is not None and mode == FileMode.write:
985 gmsh.view.remove(found[0][0])
986 found = None
987 if found is None:
988 view = gmsh.view.add(part_name)
989 step = 0
990 else:
991 view, times = found[0][0], found[1]
992 if (existing := gmsh.view.getModelData(view, 0)[0]) != data_type:
993 raise ValueError(f"'{part_name}' in {filename} is {existing}.")
994 same_time = np.flatnonzero(np.isclose(times, step_time))
995 step = int(same_time[0]) if len(same_time) > 0 else len(times)
996 gmsh.view.addModelData(
997 view,
998 step,
999 gmsh.model.getCurrent(),
1000 data_type,
1001 np.arange(1, part.shape[0] + 1),
1002 part.astype(np.float64),
1003 time=step_time,
1004 numComponents=part.shape[1],
1005 )
1007 gmsh.option.setNumber("Mesh.Binary", int(binary))
1008 gmsh.option.setNumber("Mesh.SaveAll", 1)
1009 gmsh.write(str(filename))
1010 # Append the views without repeating the mesh
1011 gmsh.option.setNumber("PostProcessing.SaveMesh", 0)
1012 for v in gmsh.view.getTags():
1013 if v not in preexisting_views:
1014 gmsh.view.write(v, str(filename), append=True)
1015 return None, None
1017 _run_on_root(comm, write)