Load a Gmsh .msh and use it in analysis#

This tutorial shows the minimal flow to load a .msh created by Gmsh and use it in FluxFEM analysis.

Key points#

  • ff.load_gmsh_mesh supports both hex and tet meshes.

  • Boundary conditions can be selected by the Physical Surface tag ID in Gmsh.

  • If the .msh lacks boundary faces (quad/tri), they are inferred from geometry.

Prerequisites#

  • meshio is required (already listed in pyproject.toml dependencies).

  • Tag the faces used in analysis as Physical Surface in Gmsh.

Prepare the .msh (Gmsh side)#

  1. Tag the faces used in analysis as Physical Surface. Example: set dirichlet_tag=1 for fixed boundary, traction_tag=2 for load boundary.

  2. Export the .msh (either v2 or v4 is OK).

These tag IDs are later used to select boundary conditions via facet_tags.

Load in FluxFEM#

import numpy as np
import jax.numpy as jnp
import fluxfem as ff
import fluxfem.helpers_wf as h_wf

mesh, facets, facet_tags = ff.load_gmsh_mesh("mesh.msh")

# Switch space by tet/hex
if isinstance(mesh, ff.TetMesh):
    space = ff.make_tet_space(mesh, dim=3, intorder=2)
else:
    space = ff.make_hex_space(mesh, dim=3, intorder=2)

# 3D linear elastic material matrix
E = 210_000.0
nu = 0.3
D = ff.isotropic_3d_D(E, nu)

bilinear_form = ff.BilinearForm.volume(
    lambda u, v, D_mat: h_wf.ddot(v.sym_grad, h_wf.matmul_std(D_mat, u.sym_grad))
    * h_wf.dOmega()
)
K = space.assemble(bilinear_form, params=D)

Boundary conditions (Physical Surface tags)#

Extract boundaries from facet_tags and set Dirichlet and Neumann conditions.

if facets is None or facet_tags is None:
    raise RuntimeError("Boundary faces or tags are not included in the .msh.")

dirichlet_tag = 1
traction_tag = 2

facets_np = np.asarray(facets)
tags_np = np.asarray(facet_tags)

# Dirichlet: fix nodes on tagged faces
dirichlet_nodes = np.unique(facets_np[tags_np == dirichlet_tag].reshape(-1))
dir_dofs = mesh.node_dofs(dirichlet_nodes, components="xyz")

# Neumann: surface traction on tagged faces
neumann_facets = facets_np[tags_np == traction_tag]
surf = ff.make_surface_from_facets(np.asarray(mesh.coords), neumann_facets)
surface_form = ff.LinearForm.surface(lambda v, p: (v | p) * h_wf.ds())
traction_vec = np.array([1.0, 0.0, 0.0], dtype=float)
F = surf.assemble_linear_form_on_space(
    space, surface_form.get_compiled(), params=traction_vec
)

solver = ff.LinearSolver(method="spsolve")
u, _ = solver.solve(K, F, dirichlet=(dir_dofs, None), dirichlet_mode="condense")

Fallback when tags are missing#

If the .msh has no Physical Surface tags, you can select boundaries by coordinate. For example, fix x = xmin and apply traction on x = xmax:

coords_np = np.asarray(mesh.coords)
xmin = float(coords_np[:, 0].min())
xmax = float(coords_np[:, 0].max())

dir_dofs = mesh.boundary_dofs_where(
    lambda pts: np.isclose(pts[:, 0], xmin, atol=1e-8),
    components="xyz",
)

neumann_facets = np.asarray(
    mesh.boundary_facets_where(lambda pts: np.allclose(pts[:, 0], xmax, atol=1e-8))
)
surf = ff.make_surface_from_facets(coords_np, neumann_facets)