Mesh with GMSH

Mesh with GMSH#

GMSH is a powerful mesh generation tool that can be used to create complex geometries for FESTIM simulations. It supports a wide range of shapes, physical labels, and CAD import/export, making it ideal for defining detailed 2D or 3D domains.

In this tutorial, we will cover:

Using GMSH directly from a Python script
Converting a GMSH model into a dolfinx mesh that can be used with FESTIM
Generating a mesh from a CAD geometry (e.g. STEP file)

Tip

GMSH can be installed with conda install -c conda-forge python-gmsh

DFG 3D example#

This GMSH example is taken directly from Jørgen Dokken’s GMSH tutorial.

The geometry corresponds to the domain used in the DFG 3D CFD benchmark case. While we do not explore the physics of the benchmark here, the example serves as a practical demonstration of how to script mesh generation in GMSH, convert the resulting geometry into a DOLFINx mesh, and use it in a FESTIM simulation.

import gmsh
import numpy as np
import os

# Initialize the GMSH API
gmsh.initialize()
gmsh.model.add("DFG 3D")

# Define geometry parameters (length L, breadth B, height H, cylinder radius r)
L, B, H, r = 2.5, 0.41, 0.41, 0.05

# Create the main channel as a rectangular box
channel = gmsh.model.occ.addBox(0, 0, 0, L, B, H)

# Create the obstacle cylinder inside the channel
cylinder = gmsh.model.occ.addCylinder(0.5, 0, 0.2, 0, B, 0, r)

# Subtract cylinder from channel to get the fluid region
fluid = gmsh.model.occ.cut([(3, channel)], [(3, cylinder)])
gmsh.model.occ.synchronize()

# Mark the fluid volume for later identification
volumes = gmsh.model.getEntities(dim=3)
fluid_marker = 11
gmsh.model.addPhysicalGroup(volumes[0][0], [volumes[0][1]], fluid_marker)
gmsh.model.setPhysicalName(volumes[0][0], fluid_marker, "Fluid volume")

# Identify and tag boundary surfaces based on their center of mass
surfaces = gmsh.model.occ.getEntities(dim=2)
inlet, outlet = None, None
walls, obstacles = [], []

inlet_marker, outlet_marker = 1, 3
wall_marker, obstacle_marker = 5, 7

for dim, tag in surfaces:
    com = gmsh.model.occ.getCenterOfMass(dim, tag)
    if np.allclose(com, [0, B / 2, H / 2]):
        gmsh.model.addPhysicalGroup(dim, [tag], inlet_marker)
        gmsh.model.setPhysicalName(dim, inlet_marker, "Fluid inlet")
        inlet = tag
    elif np.allclose(com, [L, B / 2, H / 2]):
        gmsh.model.addPhysicalGroup(dim, [tag], outlet_marker)
        gmsh.model.setPhysicalName(dim, outlet_marker, "Fluid outlet")
    elif np.isclose(com[2], 0) or np.isclose(com[1], B) or \
         np.isclose(com[2], H) or np.isclose(com[1], 0):
        walls.append(tag)
    else:
        obstacles.append(tag)

# Tag wall and obstacle surfaces
gmsh.model.addPhysicalGroup(2, walls, wall_marker)
gmsh.model.setPhysicalName(2, wall_marker, "Walls")
gmsh.model.addPhysicalGroup(2, obstacles, obstacle_marker)
gmsh.model.setPhysicalName(2, obstacle_marker, "Obstacle")

# Define mesh size field to refine near the obstacle
distance = gmsh.model.mesh.field.add("Distance")
gmsh.model.mesh.field.setNumbers(distance, "FacesList", obstacles)
resolution = r / 10
threshold = gmsh.model.mesh.field.add("Threshold")
gmsh.model.mesh.field.setNumber(threshold, "IField", distance)
gmsh.model.mesh.field.setNumber(threshold, "LcMin", resolution)
gmsh.model.mesh.field.setNumber(threshold, "LcMax", 20 * resolution)
gmsh.model.mesh.field.setNumber(threshold, "DistMin", 0.5 * r)
gmsh.model.mesh.field.setNumber(threshold, "DistMax", r)

# Optionally refine mesh near inlet
inlet_dist = gmsh.model.mesh.field.add("Distance")
gmsh.model.mesh.field.setNumbers(inlet_dist, "FacesList", [inlet])
inlet_thre = gmsh.model.mesh.field.add("Threshold")
gmsh.model.mesh.field.setNumber(inlet_thre, "IField", inlet_dist)
gmsh.model.mesh.field.setNumber(inlet_thre, "LcMin", 5 * resolution)
gmsh.model.mesh.field.setNumber(inlet_thre, "LcMax", 10 * resolution)
gmsh.model.mesh.field.setNumber(inlet_thre, "DistMin", 0.1)
gmsh.model.mesh.field.setNumber(inlet_thre, "DistMax", 0.5)

# Apply the minimal field combining both refinement regions
minimum = gmsh.model.mesh.field.add("Min")
gmsh.model.mesh.field.setNumbers(minimum, "FieldsList", [threshold, inlet_thre])
gmsh.model.mesh.field.setAsBackgroundMesh(minimum)

# Synchronize and generate 3D mesh
gmsh.model.occ.synchronize()
gmsh.model.mesh.generate(3)


# Ensure the output folder exists
os.makedirs("gmsh", exist_ok=True)

# Save the mesh in GMSH format for downstream use
gmsh.write("gmsh/mesh3D.msh")

Show code cell output

Hide code cell output

Info    : [  0%] Difference                                                                                  
Info    : [ 10%] Difference                                                                                  
Info    : [ 20%] Difference                                                                                  
Info    : [ 30%] Difference                                                                                  
Info    : [ 50%] Difference                                                                                  
Info    : [ 70%] Difference                                                                                  
Info    : [ 80%] Difference - Adding holes                                                                                
Info    : [ 90%] Difference - Classify solids                                                                                
                                                                                
Info    : Meshing 1D...
Info    : [  0%] Meshing curve 13 (Circle)
Info    : [ 10%] Meshing curve 14 (Line)
Info    : [ 20%] Meshing curve 15 (Circle)
Info    : [ 30%] Meshing curve 16 (Line)
Info    : [ 30%] Meshing curve 17 (Line)
Info    : [ 40%] Meshing curve 18 (Line)
Info    : [ 50%] Meshing curve 19 (Line)
Info    : [ 50%] Meshing curve 20 (Line)
Info    : [ 60%] Meshing curve 21 (Line)
Info    : [ 70%] Meshing curve 22 (Line)
Info    : [ 70%] Meshing curve 23 (Line)
Info    : [ 80%] Meshing curve 24 (Line)
Info    : [ 90%] Meshing curve 25 (Line)
Info    : [ 90%] Meshing curve 26 (Line)
Info    : [100%] Meshing curve 27 (Line)
Info    : Done meshing 1D (Wall 0.0176903s, CPU 0.017505s)
Info    : Meshing 2D...
Info    : [  0%] Meshing surface 7 (Cylinder, Frontal-Delaunay)

Info    : [ 20%] Meshing surface 8 (Plane, Frontal-Delaunay)
Info    : [ 30%] Meshing surface 9 (Plane, Frontal-Delaunay)
Info    : [ 50%] Meshing surface 10 (Plane, Frontal-Delaunay)
Info    : [ 60%] Meshing surface 11 (Plane, Frontal-Delaunay)
Info    : [ 80%] Meshing surface 12 (Plane, Frontal-Delaunay)
Info    : [ 90%] Meshing surface 13 (Plane, Frontal-Delaunay)
Info    : Done meshing 2D (Wall 0.38194s, CPU 0.375832s)
Info    : Meshing 3D...
Info    : 3D Meshing 1 volume with 1 connected component
Info    : Tetrahedrizing 10220 nodes...

Info    : Done tetrahedrizing 10228 nodes (Wall 0.166304s, CPU 0.161614s)
Info    : Reconstructing mesh...
Info    :  - Creating surface mesh
Info    :  - Identifying boundary edges
Info    :  - Recovering boundary

Info    : Done reconstructing mesh (Wall 0.419184s, CPU 0.400162s)
Info    : Found volume 1

Info    : It. 0 - 0 nodes created - worst tet radius 19.0147 (nodes removed 0 0)
Info    : It. 500 - 500 nodes created - worst tet radius 4.00155 (nodes removed 0 0)

Info    : It. 1000 - 1000 nodes created - worst tet radius 2.83444 (nodes removed 0 0)
Info    : It. 1500 - 1500 nodes created - worst tet radius 2.47059 (nodes removed 0 0)
Info    : It. 2000 - 2000 nodes created - worst tet radius 2.25416 (nodes removed 0 0)
Info    : It. 2500 - 2500 nodes created - worst tet radius 2.10153 (nodes removed 0 0)
Info    : It. 3000 - 3000 nodes created - worst tet radius 1.97652 (nodes removed 0 0)

Info    : It. 3500 - 3500 nodes created - worst tet radius 1.88378 (nodes removed 0 0)
Info    : It. 4000 - 4000 nodes created - worst tet radius 1.80401 (nodes removed 0 0)
Info    : It. 4500 - 4500 nodes created - worst tet radius 1.73927 (nodes removed 0 0)
Info    : It. 5000 - 5000 nodes created - worst tet radius 1.68427 (nodes removed 0 0)
Info    : It. 5500 - 5500 nodes created - worst tet radius 1.63021 (nodes removed 0 0)

Info    : It. 6000 - 6000 nodes created - worst tet radius 1.68971 (nodes removed 0 0)
Info    : It. 6500 - 6500 nodes created - worst tet radius 1.55032 (nodes removed 0 0)
Info    : It. 7000 - 7000 nodes created - worst tet radius 1.51584 (nodes removed 0 0)
Info    : It. 7500 - 7500 nodes created - worst tet radius 1.48152 (nodes removed 0 0)
Info    : It. 8000 - 8000 nodes created - worst tet radius 1.44974 (nodes removed 0 0)
Info    : It. 8500 - 8500 nodes created - worst tet radius 1.42403 (nodes removed 0 0)

Info    : It. 9000 - 9000 nodes created - worst tet radius 1.39914 (nodes removed 0 0)
Info    : It. 9500 - 9500 nodes created - worst tet radius 1.37613 (nodes removed 0 0)
Info    : It. 10000 - 10000 nodes created - worst tet radius 1.3548 (nodes removed 0 0)
Info    : It. 10500 - 10500 nodes created - worst tet radius 1.33339 (nodes removed 0 0)
Info    : It. 11000 - 11000 nodes created - worst tet radius 1.31304 (nodes removed 0 0)
Info    : It. 11500 - 11500 nodes created - worst tet radius 1.2962 (nodes removed 0 0)

Info    : It. 12000 - 12000 nodes created - worst tet radius 1.27909 (nodes removed 0 0)
Info    : It. 12500 - 12500 nodes created - worst tet radius 1.26284 (nodes removed 0 0)
Info    : It. 13000 - 13000 nodes created - worst tet radius 1.24715 (nodes removed 0 0)
Info    : It. 13500 - 13500 nodes created - worst tet radius 1.23131 (nodes removed 0 0)
Info    : It. 14000 - 14000 nodes created - worst tet radius 1.21714 (nodes removed 0 0)

Info    : It. 14500 - 14500 nodes created - worst tet radius 1.20466 (nodes removed 0 0)
Info    : It. 15000 - 15000 nodes created - worst tet radius 1.19238 (nodes removed 0 0)
Info    : It. 15500 - 15500 nodes created - worst tet radius 1.18012 (nodes removed 0 0)
Info    : It. 16000 - 16000 nodes created - worst tet radius 1.17075 (nodes removed 0 0)
Info    : It. 16500 - 16500 nodes created - worst tet radius 1.15701 (nodes removed 0 0)
Info    : It. 17000 - 17000 nodes created - worst tet radius 1.14583 (nodes removed 0 0)

Info    : It. 17500 - 17500 nodes created - worst tet radius 1.1918 (nodes removed 0 0)
Info    : It. 18000 - 18000 nodes created - worst tet radius 1.12567 (nodes removed 0 0)
Info    : It. 18500 - 18500 nodes created - worst tet radius 1.11558 (nodes removed 0 0)
Info    : It. 19000 - 19000 nodes created - worst tet radius 1.10585 (nodes removed 0 0)

Info    : It. 19500 - 19500 nodes created - worst tet radius 1.09623 (nodes removed 0 0)
Info    : It. 20000 - 20000 nodes created - worst tet radius 1.08762 (nodes removed 0 0)
Info    : It. 20500 - 20500 nodes created - worst tet radius 1.07828 (nodes removed 0 0)
Info    : It. 21000 - 21000 nodes created - worst tet radius 1.07112 (nodes removed 0 0)
Info    : It. 21500 - 21500 nodes created - worst tet radius 1.06352 (nodes removed 0 0)
Info    : It. 22000 - 22000 nodes created - worst tet radius 1.05604 (nodes removed 0 0)

Info    : It. 22500 - 22500 nodes created - worst tet radius 1.04834 (nodes removed 0 0)
Info    : It. 23000 - 23000 nodes created - worst tet radius 1.04087 (nodes removed 0 0)
Info    : It. 23500 - 23500 nodes created - worst tet radius 1.03381 (nodes removed 0 0)
Info    : It. 24000 - 24000 nodes created - worst tet radius 1.02679 (nodes removed 0 0)
Info    : It. 24500 - 24500 nodes created - worst tet radius 1.01997 (nodes removed 0 0)
Info    : It. 25000 - 25000 nodes created - worst tet radius 1.0419 (nodes removed 0 0)

Info    : It. 25500 - 25500 nodes created - worst tet radius 1.00758 (nodes removed 0 0)
Info    : It. 26000 - 26000 nodes created - worst tet radius 1.00108 (nodes removed 0 0)
Info    : 3D refinement terminated (36305 nodes total):
Info    :  - 0 Delaunay cavities modified for star shapeness
Info    :  - 0 nodes could not be inserted
Info    :  - 197061 tetrahedra created in 2.3427 sec. (84116 tets/s)
Info    : 267 node relocations

Info    : Done meshing 3D (Wall 3.47088s, CPU 3.46117s)
Info    : Optimizing mesh...
Info    : Optimizing volume 1

Info    : Optimization starts (volume = 0.417034) with worst = 0.00560819 / average = 0.764714:
Info    : 0.00 < quality < 0.10 :       429 elements
Info    : 0.10 < quality < 0.20 :      1338 elements
Info    : 0.20 < quality < 0.30 :      2543 elements
Info    : 0.30 < quality < 0.40 :      4130 elements
Info    : 0.40 < quality < 0.50 :      7055 elements
Info    : 0.50 < quality < 0.60 :     11771 elements
Info    : 0.60 < quality < 0.70 :     22477 elements
Info    : 0.70 < quality < 0.80 :     45585 elements
Info    : 0.80 < quality < 0.90 :     67291 elements
Info    : 0.90 < quality < 1.00 :     34442 elements
Info    : 4170 edge swaps, 206 node relocations (volume = 0.417034): worst = 0.158901 / average = 0.777081 (Wall 0.133574s, CPU 0.133915s)
Info    : 4229 edge swaps, 222 node relocations (volume = 0.417034): worst = 0.163219 / average = 0.777158 (Wall 0.176119s, CPU 0.176555s)
Info    : 4231 edge swaps, 226 node relocations (volume = 0.417034): worst = 0.163219 / average = 0.777166 (Wall 0.216531s, CPU 0.217126s)
Info    : No ill-shaped tets in the mesh :-)
Info    : 0.00 < quality < 0.10 :         0 elements
Info    : 0.10 < quality < 0.20 :         1 elements
Info    : 0.20 < quality < 0.30 :        13 elements
Info    : 0.30 < quality < 0.40 :      4174 elements
Info    : 0.40 < quality < 0.50 :      6962 elements
Info    : 0.50 < quality < 0.60 :     11613 elements
Info    : 0.60 < quality < 0.70 :     22610 elements
Info    : 0.70 < quality < 0.80 :     46146 elements
Info    : 0.80 < quality < 0.90 :     67766 elements
Info    : 0.90 < quality < 1.00 :     34142 elements

Info    : Done optimizing mesh (Wall 0.578277s, CPU 0.554622s)
Info    : 36305 nodes 214413 elements
Info    : Writing 'gmsh/mesh3D.msh'...

Info    : Done writing 'gmsh/mesh3D.msh'

Reading GMSH models#

DOLFINx provides convenient tools to convert GMSH models directly into DOLFINx meshes and associated mesh tags, which can then be used within FESTIM.

The function gmshio.model_to_mesh() takes a GMSH model object and converts it into a DOLFINx mesh along with cell and facet markers. This is useful when working directly with GMSH from a Python script, without writing intermediate files.

from dolfinx.io import gmsh as gmshio
from mpi4py import MPI

model_rank = 0
mesh_data = gmshio.model_to_mesh(
    gmsh.model, MPI.COMM_WORLD, model_rank
)

Alternatively, if you have saved your mesh to a .msh file, you can load it later using gmshio.read_from_msh(), specifying the mesh dimension (gdim).

mesh_data = gmshio.read_from_msh(
    "gmsh/mesh3D.msh", MPI.COMM_WORLD, 0, gdim=3
)
mesh = mesh_data.mesh
assert mesh_data.facet_tags is not None
facet_tags = mesh_data.facet_tags
facet_tags.name = "Facet markers"

assert mesh_data.cell_tags is not None
cell_tags = mesh_data.cell_tags
cell_tags.name = "Cell markers"

After loading the mesh and mesh tags, you can inspect the unique identifiers assigned to cells and facets by printing their values. This helps verify that physical groups have been correctly imported.

print(f"Cell tags: {np.unique(cell_tags.values)}")
print(f"Facet tags: {np.unique(facet_tags.values)}")

Cell tags: [11]
Facet tags: [1 3 5 7]

You can visualise the mesh along with the cell and facet tags using PyVista. This provides an intuitive way to inspect the mesh structure and verify that subdomains and boundaries are correctly marked.

2026-07-15 17:33:37.086 (   0.877s) [    73067F8E4440]vtkXOpenGLRenderWindow.:1458  WARN| bad X server connection. DISPLAY=

FESTIM Model Setup#

We now present a complete FESTIM simulation using the mesh generated from GMSH.

The steady-state diffusion equation to solve is:

\[ \nabla \cdot (D \nabla c) = 0 \]

where the diffusion coefficient, \(D=1\).

The Dirichlet boundary conditions are applied as follows:

\[ c = 1 \quad \text{on} \ \Gamma_{\mathrm{top}} \]

\[ c = 2 \quad \text{on} \ \Gamma_{\mathrm{bottom}} \]

\[ c = 0 \quad \text{on} \ \Gamma_{\mathrm{obstacle}} \]

Here, \(\Gamma_{\mathrm{top}}\), \(\Gamma_{\mathrm{bottom}}\), and \(\Gamma_{\mathrm{obstacle}}\) correspond to the physical boundaries marked in the mesh.

import festim as F

material = F.Material(D_0=1, E_D=0)

top_volume = F.VolumeSubdomain(id=11, material=material)

tube_surf = F.SurfaceSubdomain(id=7)
walls = F.SurfaceSubdomain(id=5)
top_surface = F.SurfaceSubdomain(id=1)
bottom_surface = F.SurfaceSubdomain(id=3)

my_model = F.HydrogenTransportProblem()

my_model.mesh = F.Mesh(mesh)

# we need to pass the meshtags to the model directly
my_model.facet_meshtags = facet_tags
my_model.volume_meshtags = cell_tags

my_model.subdomains = [top_surface, bottom_surface, tube_surf, walls, top_volume]

H = F.Species("H")
my_model.species = [H]

my_model.temperature = 400

my_model.boundary_conditions = [
    F.FixedConcentrationBC(subdomain=tube_surf, value=0, species=H),
    F.FixedConcentrationBC(subdomain=top_surface, value=1, species=H),
    F.FixedConcentrationBC(subdomain=bottom_surface, value=2, species=H),
]

my_model.settings = F.Settings(atol=1e-10, rtol=1e-10, transient=False)

my_model.initialise()
my_model.run()

Visualisation#

Import CAD in GMSH#

For complex geometries, GMSH allows importing CAD files such as STEP or IGES formats.

In this example, we use a CAD model from the GMSH tutorial, generate a mesh from it, and then import the mesh into a FESTIM simulation.

import gmsh
import os

gmsh.initialize()

# download cad from https://gitlab.onelab.info/gmsh/gmsh/-/raw/gmsh_4_8_4/tutorial/t20_data.step?inline=false
import requests

if not os.path.exists(os.path.join(os.pardir, "gmsh/t20_data.step")):
    url = "https://gitlab.onelab.info/gmsh/gmsh/-/raw/gmsh_4_8_4/tutorial/t20_data.step?inline=false"
    response = requests.get(url)
    with open("gmsh/t20_data.step", "wb") as f:
        f.write(response.content)

gmsh.model.add("t20")
v = gmsh.model.occ.importShapes("gmsh/t20_data.step")

gmsh.model.occ.synchronize()
volumes = gmsh.model.getEntities(dim=3)
vol_marker = 1
gmsh.model.addPhysicalGroup(volumes[0][0], [volumes[0][1]], vol_marker)
gmsh.model.setPhysicalName(volumes[0][0], vol_marker, "Volume")

surfaces = gmsh.model.occ.getEntities(dim=2)
gmsh.model.addPhysicalGroup(2, [surfaces[0][1]], 1)
gmsh.model.setPhysicalName(2, 1, "Surf1")

gmsh.model.addPhysicalGroup(2, [surfaces[3][1]], 2)
gmsh.model.setPhysicalName(2, 2, "Surf2")

# Finally, let's specify a global mesh size and mesh the partitioned model:
gmsh.option.setNumber("Mesh.MeshSizeMin", 3)
gmsh.option.setNumber("Mesh.MeshSizeMax", 3)
gmsh.model.mesh.generate(3)
gmsh.write("gmsh/t20.msh")
gmsh.finalize()

Show code cell output

Hide code cell output

Info    :  - Label 'Shapes/Rhino Product' (3D)
Info    : Meshing 1D...
Info    : [  0%] Meshing curve 1 (BSpline)
Info    : [ 10%] Meshing curve 2 (BSpline)
Info    : [ 10%] Meshing curve 3 (BSpline)
Info    : [ 10%] Meshing curve 4 (BSpline)
Info    : [ 10%] Meshing curve 5 (BSpline)
Info    : [ 20%] Meshing curve 6 (BSpline)
Info    : [ 20%] Meshing curve 7 (BSpline)
Info    : [ 20%] Meshing curve 8 (BSpline)
Info    : [ 20%] Meshing curve 9 (BSpline)
Info    : [ 20%] Meshing curve 10 (BSpline)
Info    : [ 30%] Meshing curve 11 (BSpline)
Info    : [ 30%] Meshing curve 12 (BSpline)
Info    : [ 30%] Meshing curve 13 (BSpline)
Info    : [ 30%] Meshing curve 14 (BSpline)
Info    : [ 30%] Meshing curve 15 (BSpline)
Info    : [ 40%] Meshing curve 16 (BSpline)
Info    : [ 40%] Meshing curve 17 (BSpline)
Info    : [ 40%] Meshing curve 18 (BSpline)
Info    : [ 40%] Meshing curve 19 (BSpline)
Info    : [ 40%] Meshing curve 20 (BSpline)
Info    : [ 50%] Meshing curve 21 (BSpline)
Info    : [ 50%] Meshing curve 22 (BSpline)
Info    : [ 50%] Meshing curve 23 (BSpline)
Info    : [ 50%] Meshing curve 24 (BSpline)
Info    : [ 60%] Meshing curve 25 (BSpline)
Info    : [ 60%] Meshing curve 26 (BSpline)
Info    : [ 60%] Meshing curve 27 (BSpline)
Info    : [ 60%] Meshing curve 28 (BSpline)
Info    : [ 60%] Meshing curve 29 (BSpline)
Info    : [ 70%] Meshing curve 30 (BSpline)
Info    : [ 70%] Meshing curve 31 (BSpline)
Info    : [ 70%] Meshing curve 32 (BSpline)
Info    : [ 70%] Meshing curve 33 (BSpline)
Info    : [ 70%] Meshing curve 34 (BSpline)
Info    : [ 80%] Meshing curve 35 (BSpline)
Info    : [ 80%] Meshing curve 36 (BSpline)
Info    : [ 80%] Meshing curve 37 (BSpline)
Info    : [ 80%] Meshing curve 38 (BSpline)
Info    : [ 80%] Meshing curve 39 (BSpline)
Info    : [ 90%] Meshing curve 40 (BSpline)
Info    : [ 90%] Meshing curve 41 (BSpline)
Info    : [ 90%] Meshing curve 42 (BSpline)
Info    : [ 90%] Meshing curve 43 (BSpline)
Info    : [ 90%] Meshing curve 44 (BSpline)
Info    : [100%] Meshing curve 45 (BSpline)
Info    : [100%] Meshing curve 46 (BSpline)
Info    : [100%] Meshing curve 47 (BSpline)
Info    : [100%] Meshing curve 48 (BSpline)
Info    : Done meshing 1D (Wall 0.0923492s, CPU 0.093768s)
Info    : Meshing 2D...
Info    : [  0%] Meshing surface 1 (BSpline surface, Frontal-Delaunay)
Info    : [ 10%] Meshing surface 2 (BSpline surface, Frontal-Delaunay)
Info    : [ 10%] Meshing surface 3 (BSpline surface, Frontal-Delaunay)
Info    : [ 20%] Meshing surface 4 (BSpline surface, Frontal-Delaunay)
Info    : [ 20%] Meshing surface 5 (BSpline surface, Frontal-Delaunay)
Info    : [ 30%] Meshing surface 6 (BSpline surface, Frontal-Delaunay)
Info    : [ 30%] Meshing surface 7 (BSpline surface, Frontal-Delaunay)
Info    : [ 40%] Meshing surface 8 (BSpline surface, Frontal-Delaunay)
Info    : [ 40%] Meshing surface 9 (BSpline surface, Frontal-Delaunay)
Info    : [ 50%] Meshing surface 10 (BSpline surface, Frontal-Delaunay)
Info    : [ 50%] Meshing surface 11 (BSpline surface, Frontal-Delaunay)
Info    : [ 60%] Meshing surface 12 (BSpline surface, Frontal-Delaunay)
Info    : [ 60%] Meshing surface 13 (BSpline surface, Frontal-Delaunay)
Info    : [ 70%] Meshing surface 14 (BSpline surface, Frontal-Delaunay)
Info    : [ 70%] Meshing surface 15 (BSpline surface, Frontal-Delaunay)
Info    : [ 80%] Meshing surface 16 (BSpline surface, Frontal-Delaunay)
Info    : [ 80%] Meshing surface 17 (BSpline surface, Frontal-Delaunay)
Info    : [ 90%] Meshing surface 18 (BSpline surface, Frontal-Delaunay)
Info    : [ 90%] Meshing surface 19 (BSpline surface, Frontal-Delaunay)
Info    : [100%] Meshing surface 20 (BSpline surface, Frontal-Delaunay)
Info    : [100%] Meshing surface 21 (BSpline surface, Frontal-Delaunay)
Info    : Done meshing 2D (Wall 0.055716s, CPU 0.056961s)
Info    : Meshing 3D...
Info    : 3D Meshing 1 volume with 1 connected component
Info    : Tetrahedrizing 1081 nodes...
Info    : Done tetrahedrizing 1089 nodes (Wall 0.00973786s, CPU 0.009796s)
Info    : Reconstructing mesh...
Info    :  - Creating surface mesh
Info    :  - Identifying boundary edges
Info    :  - Recovering boundary
Info    : Done reconstructing mesh (Wall 0.0209287s, CPU 0.021252s)
Info    : Found volume 1
Info    : It. 0 - 0 nodes created - worst tet radius 1.55238 (nodes removed 0 0)
Info    : 3D refinement terminated (1312 nodes total):
Info    :  - 0 Delaunay cavities modified for star shapeness
Info    :  - 12 nodes could not be inserted
Info    :  - 4669 tetrahedra created in 0.0098332 sec. (474819 tets/s)
Info    : 0 node relocations
Info    : Done meshing 3D (Wall 0.0463324s, CPU 0.047119s)
Info    : Optimizing mesh...
Info    : Optimizing volume 1

Info    : Optimization starts (volume = 18459.8) with worst = 0.0192265 / average = 0.731441:
Info    : 0.00 < quality < 0.10 :        19 elements
Info    : 0.10 < quality < 0.20 :        59 elements
Info    : 0.20 < quality < 0.30 :        91 elements
Info    : 0.30 < quality < 0.40 :        91 elements
Info    : 0.40 < quality < 0.50 :       176 elements
Info    : 0.50 < quality < 0.60 :       320 elements
Info    : 0.60 < quality < 0.70 :       752 elements
Info    : 0.70 < quality < 0.80 :      1338 elements
Info    : 0.80 < quality < 0.90 :      1311 elements
Info    : 0.90 < quality < 1.00 :       512 elements
Info    : 136 edge swaps, 1 node relocations (volume = 18459.8): worst = 0.0873398 / average = 0.746397 (Wall 0.00183683s, CPU 0.002113s)
Info    : 137 edge swaps, 2 node relocations (volume = 18459.8): worst = 0.0873398 / average = 0.746527 (Wall 0.00234748s, CPU 0.002675s)
Info    : No ill-shaped tets in the mesh :-)
Info    : 0.00 < quality < 0.10 :        12 elements
Info    : 0.10 < quality < 0.20 :        24 elements
Info    : 0.20 < quality < 0.30 :         6 elements
Info    : 0.30 < quality < 0.40 :       100 elements
Info    : 0.40 < quality < 0.50 :       171 elements
Info    : 0.50 < quality < 0.60 :       311 elements
Info    : 0.60 < quality < 0.70 :       739 elements
Info    : 0.70 < quality < 0.80 :      1359 elements
Info    : 0.80 < quality < 0.90 :      1316 elements
Info    : 0.90 < quality < 1.00 :       518 elements
Info    : Done optimizing mesh (Wall 0.00566312s, CPU 0.005975s)
Info    : 1312 nodes 7056 elements
Info    : Writing 'gmsh/t20.msh'...
Info    : Done writing 'gmsh/t20.msh'

FESTIM model#

model_rank = 0
mesh_data = gmshio.read_from_msh(
    "gmsh/t20.msh", MPI.COMM_WORLD, model_rank, gdim=3
)

mesh = mesh_data.mesh
assert mesh_data.facet_tags is not None
facet_tags = mesh_data.facet_tags
facet_tags.name = "Facet markers"
assert mesh_data.cell_tags is not None
cell_tags = mesh_data.cell_tags
cell_tags.name = "Cell markers"

print(f"Cell tags: {np.unique(cell_tags.values)}")
print(f"Facet tags: {np.unique(facet_tags.values)}")

my_model = F.HydrogenTransportProblem()

my_model.mesh = F.Mesh(mesh)

material = F.Material(D_0=1, E_D=0)

vol = F.VolumeSubdomain(id=1, material=material)

surf1 = F.SurfaceSubdomain(id=1)
surf2 = F.SurfaceSubdomain(id=2)

# we need to pass the meshtags to the model directly
my_model.facet_meshtags = facet_tags
my_model.volume_meshtags = cell_tags

my_model.subdomains = [surf1, surf2, vol]

H = F.Species("H")
my_model.species = [H]

my_model.temperature = 400

my_model.boundary_conditions = [
    F.FixedConcentrationBC(subdomain=surf1, value=1, species=H),
    F.FixedConcentrationBC(subdomain=surf2, value=0, species=H),
]

my_model.settings = F.Settings(atol=1e-10, rtol=1e-10, transient=False)

my_model.initialise()
my_model.run()

Info    : Reading 'gmsh/t20.msh'...
Info    : 98 entities
Info    : 1312 nodes
Info    : 4806 elements
Info    : Done reading 'gmsh/t20.msh'
Cell tags: [1]
Facet tags: [1 2]

Visualisation#

Meshing a multi-volume geometry#

This section discusses how to mesh multiple volumes in GMSH, which users may need to do for a multi-material hydrogen transport problem in FESTIM. We use CadQuery, a powerful Python-based library that can build 3D parametric models. See its documentation to learn more.

First, we use CadQuery to create a U-shaped pipe and then export to .brep:

import cadquery as cq

tube_r = 5      # inner radius
wall = 2         # wall thickness
leg_length = 40  # vertical leg length
bend_r = 15      # centerline bend radius

path = (
    cq.Workplane("XZ")
    .vLine(leg_length)
    .radiusArc((2*bend_r, leg_length), bend_r)  # 180° bend
    .vLine(-leg_length)
    .consolidateWires()
)

# Inner (fluid) volume
inner = (
    cq.Workplane("XY")
    .circle(tube_r)
    .sweep(path)
)

# Outer solid
outer = (
    cq.Workplane("XY")
    .circle(tube_r + wall)
    .sweep(path)
)

# Tube wall
shell = outer.cut(inner)

assembly = cq.Assembly()
assembly.add(shell, name="wall", color=cq.Color("gray"))
assembly.add(inner, name="fluid", color=cq.Color("blue"))

assembly.toCompound().exportBrep("tube.brep")

This should create an output file named tube.brep.

Note

In this tutorial we do not run the code section above, since cadquery and dolfinx.0.10 has a dependency conflict. We recommend having a separate environment to run CadQuery, as shown in this repo.

Now, we mesh the CAD in GMSH. First, we add the CAD and get the volumes/surfaces using getEntities (3 for volumes, 2 for surfaces):

import gmsh

gmsh.initialize()
gmsh.option.setString("Geometry.OCCTargetUnit", "M")
gmsh.model.add("tube")

gmsh.merge("tube.brep")
gmsh.model.occ.synchronize()

volumes = gmsh.model.getEntities(dim=3)
surfaces = gmsh.model.getEntities(dim=2)

print(f"Found {len(volumes)} volumes")
print(f"Found {len(surfaces)} surfaces")

Info    : Reading 'tube.brep'...
Info    : Done reading 'tube.brep'
Found 2 volumes
Found 10 surfaces

Tip

For complex geometries, it may be helpful to run gmsh.model.occ.fragment to break up volumes and splitting overlapping entities. We do not use it here since our geometry is fairly simple.

Here we separate individual volumes and surfaces for physical grouping:

fluid_vol = volumes[1][1]
wall_vol  = volumes[0][1]

outer_surfaces = [1, 2, 3, 5, 7]
interfaces     = [4, 6, 8]
inlet          = [9]
outlet         = [10]

# Setting markers

fluid_vol_marker = 2
wall_vol_marker = 1
outer_surfaces_marker = 10
interfaces_marker = 11
inlet_marker = 12
outlet_marker = 13

Tip

You can find these surface IDs by running the GMSH script without adding physical groups. Once you create your mesh, open it in GMSH and go to Tools->Visibility to get the correct IDs for the corresponding surfaces and volumes. Then, use those tags in your GMSH script to make physical groups.

Let’s assign physical groups and markers:

# Volumes
gmsh.model.addPhysicalGroup(3, [fluid_vol], tag=fluid_vol_marker, name="fluid")
gmsh.model.addPhysicalGroup(3, [wall_vol], tag=wall_vol_marker, name="wall")

# Surfaces
gmsh.model.addPhysicalGroup(2, outer_surfaces, tag=outer_surfaces_marker, name="outer_surfaces")
gmsh.model.addPhysicalGroup(2, interfaces, tag=interfaces_marker, name="fluid_wall_interface")
gmsh.model.addPhysicalGroup(2, inlet, tag=inlet_marker, name="inlet")
gmsh.model.addPhysicalGroup(2, outlet, tag=outlet_marker, name="outlet")

For curved meshes, users can use MeshSizeFromCurvature to specify element sizes (higher values lead to more refinement):

gmsh.option.setNumber("Mesh.MeshSizeFromCurvature", 20)

Now we generate the mesh:

gmsh.model.mesh.generate(3)
gmsh.write("tube.msh")
gmsh.finalize()

Show code cell output

Hide code cell output

Info    : Meshing 1D...
Info    : [  0%] Meshing curve 1 (Line)
Info    : [ 10%] Meshing curve 2 (Circle)
Info    : [ 20%] Meshing curve 3 (Circle)
Info    : [ 30%] Meshing curve 4 (Circle)
Info    : [ 30%] Meshing curve 5 (Circle)
Info    : [ 40%] Meshing curve 6 (Circle)
Info    : [ 50%] Meshing curve 7 (Line)
Info    : [ 60%] Meshing curve 8 (Circle)
Info    : [ 60%] Meshing curve 9 (Line)
Info    : [ 70%] Meshing curve 10 (Circle)
Info    : [ 80%] Meshing curve 11 (Circle)
Info    : [ 80%] Meshing curve 12 (Circle)
Info    : [ 90%] Meshing curve 13 (Circle)
Info    : [100%] Meshing curve 14 (Line)
Info    : Done meshing 1D (Wall 0.00421445s, CPU 0.004753s)
Info    : Meshing 2D...
Info    : [  0%] Meshing surface 1 (Unknown, Frontal-Delaunay)
Info    : [ 20%] Meshing surface 2 (Plane, Frontal-Delaunay)
Info    : [ 30%] Meshing surface 3 (Unknown, Frontal-Delaunay)
Info    : [ 40%] Meshing surface 4 (Unknown, Frontal-Delaunay)
Info    : [ 50%] Meshing surface 5 (Unknown, Frontal-Delaunay)
Info    : [ 60%] Meshing surface 6 (Unknown, Frontal-Delaunay)
Info    : [ 70%] Meshing surface 7 (Plane, Frontal-Delaunay)
Info    : [ 80%] Meshing surface 8 (Unknown, Frontal-Delaunay)
Info    : [ 90%] Meshing surface 9 (Plane, Frontal-Delaunay)
Info    : [100%] Meshing surface 10 (Plane, Frontal-Delaunay)

Info    : Done meshing 2D (Wall 0.172785s, CPU 0.173713s)
Info    : Meshing 3D...
Info    : 3D Meshing 2 volumes with 1 connected component
Info    : Tetrahedrizing 3593 nodes...

Info    : Done tetrahedrizing 3601 nodes (Wall 0.0416954s, CPU 0.042537s)
Info    : Reconstructing mesh...
Info    :  - Creating surface mesh
Info    :  - Identifying boundary edges
Info    :  - Recovering boundary
Info    : Done reconstructing mesh (Wall 0.103711s, CPU 0.103004s)
Info    : Found volume 2
Info    : Found volume 1

Info    : It. 0 - 0 nodes created - worst tet radius 3.46544 (nodes removed 0 0)

Info    : It. 500 - 500 nodes created - worst tet radius 1.22407 (nodes removed 0 0)
Info    : It. 1000 - 1000 nodes created - worst tet radius 1.04394 (nodes removed 0 0)
Info    : 3D refinement terminated (4760 nodes total):
Info    :  - 0 Delaunay cavities modified for star shapeness
Info    :  - 0 nodes could not be inserted
Info    :  - 24061 tetrahedra created in 0.077466 sec. (310600 tets/s)
Info    : 0 node relocations
Info    : Done meshing 3D (Wall 0.307145s, CPU 0.308644s)
Info    : Optimizing mesh...
Info    : Optimizing volume 1
Info    : Optimization starts (volume = 9477.7) with worst = 0.0137843 / average = 0.796014:
Info    : 0.00 < quality < 0.10 :        34 elements
Info    : 0.10 < quality < 0.20 :        51 elements
Info    : 0.20 < quality < 0.30 :        55 elements
Info    : 0.30 < quality < 0.40 :        63 elements
Info    : 0.40 < quality < 0.50 :       193 elements
Info    : 0.50 < quality < 0.60 :       279 elements
Info    : 0.60 < quality < 0.70 :       817 elements
Info    : 0.70 < quality < 0.80 :      4096 elements
Info    : 0.80 < quality < 0.90 :      2513 elements
Info    : 0.90 < quality < 1.00 :      2604 elements
Info    : 139 edge swaps, 0 node relocations (volume = 9477.7): worst = 0.163673 / average = 0.804864 (Wall 0.00230863s, CPU 0.002384s)
Info    : 141 edge swaps, 0 node relocations (volume = 9477.7): worst = 0.163673 / average = 0.804832 (Wall 0.00307633s, CPU 0.003213s)
Info    : No ill-shaped tets in the mesh :-)
Info    : 0.00 < quality < 0.10 :         0 elements
Info    : 0.10 < quality < 0.20 :         1 elements
Info    : 0.20 < quality < 0.30 :         2 elements
Info    : 0.30 < quality < 0.40 :        67 elements
Info    : 0.40 < quality < 0.50 :       192 elements
Info    : 0.50 < quality < 0.60 :       257 elements
Info    : 0.60 < quality < 0.70 :       786 elements
Info    : 0.70 < quality < 0.80 :      4128 elements
Info    : 0.80 < quality < 0.90 :      2552 elements
Info    : 0.90 < quality < 1.00 :      2589 elements
Info    : Optimizing volume 2
Info    : Optimization starts (volume = 9866.29) with worst = 0.019376 / average = 0.764427:
Info    : 0.00 < quality < 0.10 :        39 elements
Info    : 0.10 < quality < 0.20 :        93 elements
Info    : 0.20 < quality < 0.30 :       184 elements
Info    : 0.30 < quality < 0.40 :       261 elements
Info    : 0.40 < quality < 0.50 :       385 elements
Info    : 0.50 < quality < 0.60 :       619 elements
Info    : 0.60 < quality < 0.70 :      1578 elements
Info    : 0.70 < quality < 0.80 :      3466 elements
Info    : 0.80 < quality < 0.90 :      4748 elements
Info    : 0.90 < quality < 1.00 :      1983 elements
Info    : 316 edge swaps, 7 node relocations (volume = 9866.29): worst = 0.188728 / average = 0.77935 (Wall 0.00531054s, CPU 0.004502s)
Info    : 319 edge swaps, 7 node relocations (volume = 9866.29): worst = 0.300579 / average = 0.779449 (Wall 0.00639394s, CPU 0.005664s)
Info    : No ill-shaped tets in the mesh :-)
Info    : 0.00 < quality < 0.10 :         0 elements
Info    : 0.10 < quality < 0.20 :         0 elements
Info    : 0.20 < quality < 0.30 :         0 elements
Info    : 0.30 < quality < 0.40 :       259 elements
Info    : 0.40 < quality < 0.50 :       355 elements
Info    : 0.50 < quality < 0.60 :       586 elements
Info    : 0.60 < quality < 0.70 :      1611 elements
Info    : 0.70 < quality < 0.80 :      3478 elements
Info    : 0.80 < quality < 0.90 :      4780 elements
Info    : 0.90 < quality < 1.00 :      1997 elements
Info    : Done optimizing mesh (Wall 0.0281524s, CPU 0.027857s)
Info    : 4760 nodes 31154 elements
Info    : Writing 'tube.msh'...
Info    : Done writing 'tube.msh'

We should have a new file named tube.msh. Similar to the example above, we can read the mesh information for DOLFINx:

from mpi4py import MPI
import festim as F

mesh_data = gmshio.read_from_msh(
    "tube.msh", MPI.COMM_WORLD, 0, gdim=3
)
mesh = mesh_data.mesh
facet_tags = mesh_data.facet_tags
facet_tags.name = "Facet markers"

cell_tags = mesh_data.cell_tags
cell_tags.name = "Cell markers"

Info    : Reading 'tube.msh'...
Info    : 34 entities
Info    : 4760 nodes
Info    : 30864 elements
Info    : Done reading 'tube.msh'

We can now use this mesh in FESTIM for a discontinuous problem:

my_model = F.HydrogenTransportProblemDiscontinuous()
my_model.mesh = F.Mesh(mesh)

fluid = F.Material(D_0=1, E_D=0, K_S_0=1, E_K_S=0)
wall = F.Material(D_0=0.01, E_D=0, K_S_0=2, E_K_S=0)

fluid_vol = F.VolumeSubdomain(id=fluid_vol_marker, material=fluid)
wall_vol = F.VolumeSubdomain(id=wall_vol_marker, material=wall)

inlet = F.SurfaceSubdomain(id=inlet_marker)
outlet = F.SurfaceSubdomain(id=outlet_marker)
outer_surfaces = F.SurfaceSubdomain(id=outer_surfaces_marker)

my_model.facet_meshtags = facet_tags
my_model.volume_meshtags = cell_tags

my_model.subdomains = [inlet, outlet, outer_surfaces, fluid_vol, wall_vol]
my_model.surface_to_volume = {
    outer_surfaces: wall_vol,
    inlet: fluid_vol,
    outlet: fluid_vol
}
my_model.interfaces =[
    F.Interface(id=interfaces_marker, subdomains=[wall_vol, fluid_vol], penalty_term=1e5)
]
H = F.Species("H", subdomains=[fluid_vol, wall_vol])
my_model.species = [H]

my_model.temperature = 400

my_model.boundary_conditions = [
    F.FixedConcentrationBC(subdomain=inlet, value=1, species=H),
    F.FixedConcentrationBC(subdomain=outlet, value=0, species=H),
]

my_model.settings = F.Settings(atol=1e-10, rtol=1e-10, transient=False)
my_model.initialise()
my_model.run()

Mesh with GMSH

Contents

Mesh with GMSH#

DFG 3D example#

Reading GMSH models#

FESTIM Model Setup#

Visualisation#

Import CAD in GMSH#

FESTIM model#

Visualisation#

Meshing a multi-volume geometry#