Transfer meshtags to submeshes

Transfer meshtags to submeshes#

Author: Jørgen S. Dokken

SPDX-License-Identifier: MIT

In DOLFINx, one can create sub-meshes of entities of any co-dimension \((0,\dots,\mathrm{tdim})\). For complex meshes, it is not easy to locate subdomains or boundaries of interest. In this example, we will import a mesh from GMSH, where we have marked different parts of the domain with different markers, and transfer these markers to a submesh.

import os
import sys
import pyvista
from scifem import transfer_meshtags_to_submesh, compat
from mpi4py import MPI
import gmsh
import dolfinx
try:
    from dolfinx.io import gmsh as gmshio
except ImportError:
    import dolfinx.io.gmshio as gmshio

We start by embedding an ellipsoid within another ellipsoid using GMSH. For more details about creating this mesh in GMSH, see: FEniCS Workshop: External meshes

gmsh.initialize()
center = (0, 0, 0)
aspect_ratio = 0.5
R_i = 0.3
R_e = 0.8
inner_disk = gmsh.model.occ.addDisk(*center, R_i, aspect_ratio * R_i)
outer_disk = gmsh.model.occ.addDisk(*center, R_e, R_e)
whole_domain, map_to_input = gmsh.model.occ.fragment(
    [(2, outer_disk)], [(2, inner_disk)])
gmsh.model.occ.synchronize()
circle_inner = [idx for (dim, idx) in map_to_input[1] if dim == 2]
circle_outer = [idx for (dim, idx) in map_to_input[0]
                if dim == 2 and idx not in circle_inner]
gmsh.model.addPhysicalGroup(2, circle_inner, tag=3)
gmsh.model.addPhysicalGroup(2, circle_outer, tag=7)
inner_boundary = gmsh.model.getBoundary(
    [(2, e) for e in circle_inner], recursive=False, oriented=False)
outer_boundary = gmsh.model.getBoundary(
    [(2, e) for e in circle_outer], recursive=False, oriented=False)
interface = [idx for (dim, idx) in inner_boundary if dim == 1]
ext_boundary = [
    idx for (dim, idx) in outer_boundary if idx not in interface and dim == 1]
gmsh.model.addPhysicalGroup(1, interface, tag=12)
gmsh.model.addPhysicalGroup(1, ext_boundary, tag=15)
15
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", 0.2)
gmsh.model.mesh.generate(2)
gmsh.model.mesh.setOrder(3)
gmsh.model.mesh.optimize("Netgen")
Info    : Meshing 1D...
Info    : [  0%] Meshing curve 1 (Ellipse)
Info    : [ 60%] Meshing curve 2 (Ellipse)
Info    : Done meshing 1D (Wall 0.00256396s, CPU 0.002952s)
Info    : Meshing 2D...
Info    : [  0%] Meshing surface 1 (Plane, Frontal-Delaunay)
Info    : [ 60%] Meshing surface 2 (Plane, Frontal-Delaunay)
Info    : Done meshing 2D (Wall 0.00183043s, CPU 0.001924s)
Info    : 90 nodes 188 elements
Info    : Meshing order 3 (curvilinear on)...
Info    : [  0%] Meshing curve 1 order 3
Info    : [ 30%] Meshing curve 2 order 3
Info    : [ 60%] Meshing surface 1 order 3
Info    : [ 80%] Meshing surface 2 order 3
Info    : Surface mesh: worst distortion = 0.426319 (0 elements in ]0, 0.2]); worst gamma = 0.847892
Info    : Done meshing order 3 (Wall 0.00154099s, CPU 0.001623s)
Info    : Optimizing mesh (Netgen)...
Info    : Done optimizing mesh (Wall 9.76e-07s, CPU 1e-06s)

Next, we read this mesh and corresponding markers into DOLFINx

mesh_data = gmshio.model_to_mesh(
gmsh.model, MPI.COMM_WORLD, 0, gdim=2)
if hasattr(mesh_data, "mesh"):
    circular_mesh = mesh_data.mesh
    cell_marker = mesh_data.cell_tags
    facet_marker = mesh_data.facet_tags
else:
    circular_mesh, cell_marker, facet_marker = mesh_data

We visualize the mesh and its markers

Hide code cell source

def plot_mesh(mesh: dolfinx.mesh.Mesh, values=None):
    """
    Given a DOLFINx mesh, create a `pyvista.UnstructuredGrid`,
    and plot it and the mesh nodes
    """
    plotter = pyvista.Plotter()
    V_linear = dolfinx.fem.functionspace(mesh, ("Lagrange", 1))
    linear_grid = pyvista.UnstructuredGrid(*dolfinx.plot.vtk_mesh(V_linear))
    if compat.cmap(mesh).degree > 1:
        ugrid = pyvista.UnstructuredGrid(*dolfinx.plot.vtk_mesh(mesh))
        if values is not None:
            ugrid.cell_data["Marker"] = values
        plotter.add_mesh(ugrid, style="points", color="b", point_size=10)
        ugrid = ugrid.tessellate()
        plotter.add_mesh(ugrid, show_edges=False)
        plotter.add_mesh(linear_grid, style="wireframe", color="black")

    else:
        if values is not None:
            linear_grid.cell_data["Marker"] = values
        plotter.add_mesh(linear_grid, show_edges=True)
    plotter.show_axes()
    plotter.view_xy()
    plotter.show()


plot_mesh(circular_mesh, cell_marker.values)
2026-08-05 06:37:39.610 (   1.118s) [    7EFFC63C0140]vtkXOpenGLRenderWindow.:1460  WARN| bad X server connection. DISPLAY=:99.0
---------------------------------------------------------------------------
FileNotFoundError                         Traceback (most recent call last)
File /dolfinx-env/lib/python3.12/site-packages/aiohttp/web_urldispatcher.py:546, in StaticResource.__init__(self, prefix, directory, name, expect_handler, chunk_size, show_index, follow_symlinks, append_version)
    545 try:
--> 546     directory = Path(directory).expanduser().resolve(strict=True)
    547 except FileNotFoundError as error:

File /usr/lib/python3.12/pathlib.py:1242, in Path.resolve(self, strict)
   1241 try:
-> 1242     s = self._flavour.realpath(self, strict=strict)
   1243 except OSError as e:

File <frozen posixpath>:432, in realpath(filename, strict)

File <frozen posixpath>:477, in _joinrealpath(path, rest, strict, seen)

FileNotFoundError: [Errno 2] No such file or directory: '/dolfinx-env/lib/python3.12/site-packages/trame_client/module/vue3-www'

The above exception was the direct cause of the following exception:

ValueError                                Traceback (most recent call last)
Cell In[8], line 27
     23     plotter.view_xy()
     24     plotter.show()
     25 
     26 
---> 27 plot_mesh(circular_mesh, cell_marker.values)

Cell In[8], line 24, in plot_mesh(mesh, values)
     20             linear_grid.cell_data["Marker"] = values
     21         plotter.add_mesh(linear_grid, show_edges=True)
     22     plotter.show_axes()
     23     plotter.view_xy()
---> 24     plotter.show()

File /dolfinx-env/lib/python3.12/site-packages/pyvista/_deprecate_positional_args.py:243, in _deprecate_positional_args.<locals>._inner_deprecate_positional_args.<locals>.inner_f(*args, **kwargs)
    239             warn_external(msg, PyVistaDeprecationWarning)
    241         warn_positional_args()
--> 243 return f(*args, **kwargs)

File /dolfinx-env/lib/python3.12/site-packages/pyvista/plotting/plotter.py:8231, in Plotter.show(self, title, window_size, interactive, auto_close, interactive_update, full_screen, screenshot, return_img, cpos, jupyter_backend, return_viewer, return_cpos, before_close_callback, store_image_depth, **kwargs)
   8228         jupyter_backend = self._theme.jupyter_backend
   8230     if jupyter_backend is None or jupyter_backend.lower() != 'none':
-> 8231         jupyter_disp = handle_plotter(self, backend=jupyter_backend, **jupyter_kwargs)
   8233 self.render()
   8235 # initial double render needed for certain passes when offscreen

File /dolfinx-env/lib/python3.12/site-packages/pyvista/jupyter/notebook.py:74, in handle_plotter(plotter, backend, screenshot, **kwargs)
     71     if backend in ['server', 'client', 'trame', 'html']:
     72         from pyvista.trame.jupyter import show_trame  # noqa: PLC0415
---> 74         return show_trame(plotter, mode=backend, **kwargs)
     76 except ImportError as e:
     77     # Trame was explicitly requested but not available
     78     _ensure_entry_points()

File /dolfinx-env/lib/python3.12/site-packages/pyvista/trame/jupyter.py:394, in show_trame(plotter, mode, name, server_proxy_enabled, server_proxy_prefix, jupyter_extension_enabled, collapse_menu, add_menu, add_menu_items, default_server_rendering, handler, **kwargs)
    391 kwargs.setdefault('height', dh)
    393 if mode == 'html':
--> 394     return EmbeddableWidget(plotter, **kwargs)
    396 if jupyter_extension_enabled is None:
    397     jupyter_extension_enabled = pv.global_theme.trame.jupyter_extension_enabled

File /dolfinx-env/lib/python3.12/site-packages/pyvista/trame/jupyter.py:137, in EmbeddableWidget.__init__(self, plotter, width, height, **kwargs)
    135     msg = 'Please install `ipywidgets`.'
    136     raise ImportError(msg)
--> 137 scene = plotter.export_html(filename=None)
    138 src = scene.getvalue().replace('"', '&quot;')
    139 # eventually we could maybe expose this, but for now make sure we're at least
    140 # consistent with matplotlib's color (light gray)

File /dolfinx-env/lib/python3.12/site-packages/pyvista/plotting/plotter.py:784, in BasePlotter.export_html(self, filename)
    781     msg = 'Please install trame dependencies: pip install "pyvista[jupyter]"'
    782     raise ImportError(msg)
--> 784 data = self.export_vtksz(filename=None)
    785 buffer = io.StringIO()
    786 write_html(data, buffer)

File /dolfinx-env/lib/python3.12/site-packages/pyvista/plotting/plotter.py:837, in BasePlotter.export_vtksz(self, filename, format)
    835 server = get_server(pv.global_theme.trame.jupyter_server_name)
    836 if not server.running:
--> 837     elegantly_launch(pv.global_theme.trame.jupyter_server_name)
    839 view = PyVistaLocalView(self, trame_server=server)
    841 content = view.export(format=format)

File /dolfinx-env/lib/python3.12/site-packages/pyvista/trame/jupyter.py:487, in elegantly_launch(*args, **kwargs)
    484 # Basically monkey patches asyncio to support this
    485 nest_asyncio2.apply()
--> 487 return asyncio.run(launch_it())

File /dolfinx-env/lib/python3.12/site-packages/nest_asyncio2.py:115, in run(main, debug, loop_factory)
    113 task = asyncio.ensure_future(main, loop=loop)
    114 try:
--> 115     return loop.run_until_complete(task)
    116 finally:
    117     if not task.done():

File /dolfinx-env/lib/python3.12/site-packages/nest_asyncio2.py:230, in _patch_loop.<locals>.run_until_complete(self, future)
    227 if self._ready:
    228     self._run_once()
--> 230 return f.result()

File /usr/lib/python3.12/asyncio/futures.py:203, in Future.result(self)
    201 self.__log_traceback = False
    202 if self._exception is not None:
--> 203     raise self._exception.with_traceback(self._exception_tb)
    204 return self._result

File /usr/lib/python3.12/asyncio/tasks.py:314, in Task.__step_run_and_handle_result(***failed resolving arguments***)
    310 try:
    311     if exc is None:
    312         # We use the `send` method directly, because coroutines
    313         # don't have `__iter__` and `__next__` methods.
--> 314         result = coro.send(None)
    315     else:
    316         result = coro.throw(exc)

File /dolfinx-env/lib/python3.12/site-packages/pyvista/trame/jupyter.py:482, in elegantly_launch.<locals>.launch_it()
    481 async def launch_it():
--> 482     await launch_server(*args, **kwargs).ready

File /dolfinx-env/lib/python3.12/site-packages/pyvista/trame/jupyter.py:215, in launch_server(server, port, host, wslink_backend, **kwargs)
    213 if server._running_stage == 0:
    214     server.controller.on_server_ready.add(on_ready)
--> 215     server.start(
    216         exec_mode='task',
    217         host=host,
    218         port=port,
    219         open_browser=False,
    220         show_connection_info=False,
    221         disable_logging=True,
    222         timeout=0,
    223         backend=wslink_backend,
    224     )
    225 # else, server is already running or launching
    226 return server

File /dolfinx-env/lib/python3.12/site-packages/trame_server/core.py:725, in Server.start(self, port, thread, open_browser, show_connection_info, disable_logging, backend, follow_symlinks, exec_mode, timeout, host, **kwargs)
    722 CoreServer.configure(options)
    724 self._running_stage = 1
--> 725 task = CoreServer.server_start(
    726     options,
    727     **{  # Do a proper merging/override
    728         **kwargs,
    729         "disableLogging": disable_logging,
    730         "backend": backend,
    731         "exec_mode": exec_mode,
    732     },
    733 )
    735 # Manage exit life cycle unless coroutine
    736 if exec_mode == "main":

File /dolfinx-env/lib/python3.12/site-packages/trame_server/protocol.py:52, in CoreServer.server_start(options, disableLogging, backend, exec_mode, **kwargs)
     43 @staticmethod
     44 def server_start(
     45     options,
   (...)     50 ):
     51     # NOTE: **kwargs to wslink's start_webserver are currently unused
---> 52     return server.start_webserver(
     53         options=options,
     54         protocol=CoreServer,
     55         disableLogging=disableLogging,
     56         backend=backend,
     57         exec_mode=exec_mode,
     58         **kwargs,
     59     )

File /dolfinx-env/lib/python3.12/site-packages/wslink/server.py:257, in start_webserver(options, protocol, disableLogging, backend, exec_mode, **_)
    254     server_config["handle_signals"] = not options.nosignalhandlers
    256 # Create the webserver and start it
--> 257 ws_server = create_webserver(server_config, backend=backend)
    259 # Once we have python 3.7 minimum, we can start the server with asyncio.run()
    260 # asyncio.run(ws_server.start())
    261 
    262 # Until then, we can start the server this way
    263 try:

File /dolfinx-env/lib/python3.12/site-packages/wslink/server.py:168, in create_webserver(server_config, backend)
    167 def create_webserver(server_config, backend="aiohttp"):
--> 168     return backends.create_webserver(server_config, backend=backend)

File /dolfinx-env/lib/python3.12/site-packages/wslink/backends/__init__.py:5, in create_webserver(server_config, backend)
      2 if backend == "aiohttp":
      3     from .aiohttp import create_webserver  # noqa: PLC0415
----> 5     return create_webserver(server_config)
      7 if backend == "generic":
      8     from .generic import create_webserver  # noqa: PLC0415

File /dolfinx-env/lib/python3.12/site-packages/wslink/backends/aiohttp/__init__.py:218, in create_webserver(server_config)
    215     return ReverseWebAppServer(server_config)
    217 # Normal web server
--> 218 return WebAppServer(server_config)

File /dolfinx-env/lib/python3.12/site-packages/wslink/backends/aiohttp/__init__.py:118, in WebAppServer.__init__(self, server_config)
    116     # Resolve / => index.html
    117     self.app.router.add_route("GET", "/", _root_handler)
--> 118     self.app.add_routes(routes)
    120 self.app[STATE_KEY] = {}

File /dolfinx-env/lib/python3.12/site-packages/aiohttp/web_app.py:379, in Application.add_routes(self, routes)
    378 def add_routes(self, routes: Iterable[AbstractRouteDef]) -> list[AbstractRoute]:
--> 379     return self.router.add_routes(routes)

File /dolfinx-env/lib/python3.12/site-packages/aiohttp/web_urldispatcher.py:1259, in UrlDispatcher.add_routes(self, routes)
   1257 registered_routes = []
   1258 for route_def in routes:
-> 1259     registered_routes.extend(route_def.register(self))
   1260 return registered_routes

File /dolfinx-env/lib/python3.12/site-packages/aiohttp/web_routedef.py:87, in StaticDef.register(self, router)
     86 def register(self, router: UrlDispatcher) -> list[AbstractRoute]:
---> 87     resource = router.add_static(self.prefix, self.path, **self.kwargs)
     88     routes = resource.get_info().get("routes", {})
     89     return list(routes.values())

File /dolfinx-env/lib/python3.12/site-packages/aiohttp/web_urldispatcher.py:1183, in UrlDispatcher.add_static(self, prefix, path, name, expect_handler, chunk_size, show_index, follow_symlinks, append_version)
   1181 if prefix.endswith("/"):
   1182     prefix = prefix[:-1]
-> 1183 resource = StaticResource(
   1184     prefix,
   1185     path,
   1186     name=name,
   1187     expect_handler=expect_handler,
   1188     chunk_size=chunk_size,
   1189     show_index=show_index,
   1190     follow_symlinks=follow_symlinks,
   1191     append_version=append_version,
   1192 )
   1193 self.register_resource(resource)
   1194 return resource

File /dolfinx-env/lib/python3.12/site-packages/aiohttp/web_urldispatcher.py:548, in StaticResource.__init__(self, prefix, directory, name, expect_handler, chunk_size, show_index, follow_symlinks, append_version)
    546     directory = Path(directory).expanduser().resolve(strict=True)
    547 except FileNotFoundError as error:
--> 548     raise ValueError(f"'{directory}' does not exist") from error
    549 if not directory.is_dir():
    550     raise ValueError(f"'{directory}' is not a directory")

ValueError: '/dolfinx-env/lib/python3.12/site-packages/trame_client/module/vue3-www' does not exist
# Next, we create a submesh, only extracting the upper half of the mesh
tdim = circular_mesh.topology.dim
submesh, cell_map, vertex_map, node_map = dolfinx.mesh.create_submesh(
    circular_mesh, tdim, dolfinx.mesh.locate_entities(circular_mesh, tdim, lambda x: x[1] > 0))

We transfer the cell markers to the submesh

sub_cell_marker, sub_cell_map = transfer_meshtags_to_submesh(
    cell_marker, submesh, vertex_map, cell_map)

and visualize it

Hide code cell source

plot_mesh(submesh, sub_cell_marker.values)