FEniCS Finite Element Method for PDEs

TL;DR — Solve partial differential equations with the Finite Element Method (FEM) using FEniCS/dolfinx. Derive the weak form, generate meshes with gmsh, apply Dirichlet/Neumann boundary conditions, solve Poisson or elasticity problems, and export results to XDMF/VTK for ParaView.

When to Use

Use this Skill when you need to:

• Solve elliptic, parabolic, or hyperbolic PDEs on complex geometries
• Implement custom weak forms for multi-physics problems
• Apply mixed Dirichlet/Neumann/Robin boundary conditions
• Perform convergence studies on successively refined meshes
• Export solutions for publication-quality visualization in ParaView

Do not use this Skill when:

• You need a quick 1D finite-difference solution → use SciPy solve_bvp
• You want a spectral method for periodic domains → use pseudo-spectral libraries
• You need GPU-accelerated large-scale CFD → consider OpenFOAM or Fluidity

Background & Key Concepts

Variational (Weak) Form

The FEM converts a strong-form PDE into an integral equation by multiplying by a test function v and integrating by parts. For Poisson's equation:

Strong form: −∇²u = f in Ω, u = uD on ΓD, ∇u·n = g on ΓN

Weak form: Find u ∈ H¹(Ω) such that for all v ∈ H¹₀(Ω): ∫_Ω ∇u·∇v dx = ∫_Ω f v dx + ∫_ΓN g v ds

Function Spaces

Convergence and Error

For Lagrange P1 elements on a quasi-uniform mesh of size h:

• L² error: O(h²) (one order above approximation degree)
• H¹ error: O(h)

Environment Setup

# Recommended: use conda with conda-forge (dolfinx + gmsh are complex to compile)
conda create -n fenics-env python=3.11 -y
conda activate fenics-env
conda install -c conda-forge fenics-dolfinx mpich petsc4py gmsh pyvista -y

# Verify installation
python -c "import dolfinx; print('dolfinx version:', dolfinx.__version__)"
python -c "import gmsh; print('gmsh version:', gmsh.__version__)"

# For Docker users (simplest approach)
docker pull dolfinx/dolfinx:stable
docker run -it --rm -v $(pwd):/work dolfinx/dolfinx:stable bash

Core Workflow

Step 1 — Poisson Equation on Unit Square

"""
Solve −∇²u = f on the unit square [0,1]×[0,1]
with homogeneous Dirichlet BC u=0 on ∂Ω.

Manufactured solution: u_exact = sin(πx)sin(πy)
Source term:           f = 2π² sin(πx)sin(πy)
"""

from mpi4py import MPI
import numpy as np
from dolfinx import mesh, fem, io
from dolfinx.fem.petsc import LinearProblem
import ufl


def solve_poisson_unit_square(n_cells: int = 32, degree: int = 1) -> dict:
    """
    Solve the Poisson equation −∇²u = f on the unit square.

    Args:
        n_cells: Number of cells in each direction (total cells = 2*n_cells²).
        degree:  Polynomial degree of Lagrange finite elements.

    Returns:
        Dictionary with keys: uh (solution), L2_error, H1_error.
    """
    # --- Mesh ---
    domain = mesh.create_unit_square(
        MPI.COMM_WORLD, n_cells, n_cells, mesh.CellType.triangle
    )

    # --- Function space ---
    V = fem.functionspace(domain, ("Lagrange", degree))

    # --- Exact solution for manufactured source term ---
    x = ufl.SpatialCoordinate(domain)
    u_exact_expr = ufl.sin(ufl.pi * x[0]) * ufl.sin(ufl.pi * x[1])
    f_expr       = 2.0 * ufl.pi**2 * u_exact_expr

    # --- Dirichlet BC: u = 0 on all boundaries ---
    def boundary_all(x):
        return (
            np.isclose(x[0], 0.0) | np.isclose(x[0], 1.0) |
            np.isclose(x[1], 0.0) | np.isclose(x[1], 1.0)
        )

    boundary_dofs = fem.locate_dofs_geometrical(V, boundary_all)
    u0 = fem.Function(V)
    u0.x.array[:] = 0.0
    bc = fem.dirichletbc(u0, boundary_dofs)

    # --- Variational problem ---
    u = ufl.TrialFunction(V)
    v = ufl.TestFunction(V)
    a = ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx
    L = f_expr * v * ufl.dx

    # --- Solve ---
    problem = LinearProblem(a, L, bcs=[bc],
                            petsc_options={"ksp_type": "cg", "pc_type": "hypre"})
    uh = problem.solve()

    # --- Compute errors ---
    diff = uh - u_exact_expr
    L2_error = float(fem.assemble_scalar(fem.form(ufl.inner(diff, diff) * ufl.dx)) ** 0.5)
    H1_error = float(fem.assemble_scalar(
        fem.form(ufl.inner(ufl.grad(diff), ufl.grad(diff)) * ufl.dx)
    ) ** 0.5)

    print(f"n_cells={n_cells}, degree={degree}: "
          f"L2={L2_error:.2e}, H1={H1_error:.2e}")
    return {"uh": uh, "L2_error": L2_error, "H1_error": H1_error}


def convergence_study() -> None:
    """Run convergence study with successive mesh refinement."""
    print("Convergence study for Poisson on unit square:")
    print(f"{'N':>6}  {'L2 error':>12}  {'H1 error':>12}  {'L2 rate':>8}")
    prev_L2 = None
    for n in [4, 8, 16, 32, 64]:
        result = solve_poisson_unit_square(n_cells=n, degree=1)
        L2 = result["L2_error"]
        rate = (np.log(prev_L2 / L2) / np.log(2.0)) if prev_L2 else float("nan")
        print(f"{n:>6}  {L2:>12.4e}  {result['H1_error']:>12.4e}  {rate:>8.2f}")
        prev_L2 = L2

Step 2 — Linear Elasticity

"""
Solve linear elasticity on a 2D beam under body force (gravity).

Strong form: −div(σ(u)) = f in Ω
σ(u) = λ tr(ε(u)) I + 2μ ε(u)   (Hooke's law, Lamé form)
ε(u) = ½(∇u + ∇uᵀ)             (small-strain tensor)
"""

from mpi4py import MPI
import numpy as np
from dolfinx import mesh, fem, io
from dolfinx.fem.petsc import LinearProblem
import ufl


def solve_linear_elasticity(
    nx: int = 40,
    ny: int = 10,
    E: float = 210e9,     # Young's modulus (Pa), steel
    nu: float = 0.3,      # Poisson's ratio
    rho: float = 7850.0,  # Density (kg/m³)
    g: float = 9.81,      # Gravitational acceleration (m/s²)
    output_xdmf: str = "elasticity.xdmf",
) -> None:
    """
    Solve 2D linear elasticity (plane stress) on a rectangular beam.

    The left end is clamped (u=0), a body force f=(0,-ρg) is applied.
    Exports solution to XDMF for ParaView.

    Args:
        nx:          Cells in x direction.
        ny:          Cells in y direction.
        E:           Young's modulus in Pa.
        nu:          Poisson's ratio.
        rho:         Material density in kg/m³.
        g:           Gravity magnitude in m/s².
        output_xdmf: Output file path.
    """
    # Lamé parameters
    lam = E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu))
    mu  = E / (2.0 * (1.0 + nu))

    domain = mesh.create_rectangle(
        MPI.COMM_WORLD,
        [np.array([0.0, 0.0]), np.array([1.0, 0.25])],
        [nx, ny],
        cell_type=mesh.CellType.triangle,
    )

    # Vector function space (displacement u ∈ R²)
    V = fem.functionspace(domain, ("Lagrange", 1, (2,)))

    # Strain and stress tensors
    def epsilon(u):
        return ufl.sym(ufl.nabla_grad(u))

    def sigma(u):
        return lam * ufl.nabla_div(u) * ufl.Identity(2) + 2 * mu * epsilon(u)

    # Body force: gravity
    f = fem.Constant(domain, np.array([0.0, -rho * g]))

    # Clamped BC on left boundary x=0
    def left_boundary(x):
        return np.isclose(x[0], 0.0)

    boundary_dofs = fem.locate_dofs_geometrical(V, left_boundary)
    u_D = fem.Function(V)
    u_D.x.array[:] = 0.0
    bc = fem.dirichletbc(u_D, boundary_dofs)

    # Variational formulation
    u = ufl.TrialFunction(V)
    v = ufl.TestFunction(V)
    a = ufl.inner(sigma(u), epsilon(v)) * ufl.dx
    L = ufl.inner(f, v) * ufl.dx

    # Solve
    problem = LinearProblem(a, L, bcs=[bc],
                            petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
    uh = problem.solve()
    uh.name = "Displacement"

    # Compute von Mises stress for visualization
    s = sigma(uh) - (1.0 / 3) * ufl.tr(sigma(uh)) * ufl.Identity(2)
    von_mises = ufl.sqrt(1.5 * ufl.inner(s, s))
    W = fem.functionspace(domain, ("DG", 0))
    vm_expr = fem.Expression(von_mises, W.element.interpolation_points())
    vm_func = fem.Function(W)
    vm_func.interpolate(vm_expr)
    vm_func.name = "VonMises"

    # Export to XDMF
    with io.XDMFFile(MPI.COMM_WORLD, output_xdmf, "w") as xdmf:
        xdmf.write_mesh(domain)
        xdmf.write_function(uh)
        xdmf.write_function(vm_func)

    max_disp = np.max(np.abs(uh.x.array))
    print(f"Max displacement: {max_disp:.4e} m")
    print(f"Solution written to {output_xdmf}")


if __name__ == "__main__":
    solve_linear_elasticity()

Step 3 — gmsh Mesh Generation and Import into dolfinx

"""
Use gmsh to generate a structured mesh of a disk with a hole,
then import into dolfinx for FEM analysis.
"""

import gmsh
import numpy as np
from mpi4py import MPI
from dolfinx.io.gmshio import model_to_mesh
from dolfinx import fem, io
import ufl
from dolfinx.fem.petsc import LinearProblem


def create_annular_mesh(
    r_inner: float = 0.2,
    r_outer: float = 1.0,
    mesh_size: float = 0.05,
    output_msh: str = "annulus.msh",
) -> None:
    """
    Create a 2D annular mesh (disk with circular hole) using gmsh.

    Args:
        r_inner:   Inner radius (hole).
        r_outer:   Outer radius.
        mesh_size: Target mesh element size.
        output_msh: Output .msh file path.
    """
    gmsh.initialize()
    gmsh.model.add("annulus")

    # Outer disk
    outer = gmsh.model.occ.addDisk(0, 0, 0, r_outer, r_outer)
    # Inner disk (hole)
    inner = gmsh.model.occ.addDisk(0, 0, 0, r_inner, r_inner)

    # Boolean cut: outer minus inner
    gmsh.model.occ.cut([(2, outer)], [(2, inner)])
    gmsh.model.occ.synchronize()

    # Mesh size field
    gmsh.model.mesh.setSize(gmsh.model.getEntities(0), mesh_size)

    # Physical groups (needed for dolfinx boundary conditions)
    surfaces = gmsh.model.getEntities(2)
    for _, tag in surfaces:
        gmsh.model.addPhysicalGroup(2, [tag], tag)
        gmsh.model.setPhysicalName(2, tag, f"domain_{tag}")

    curves = gmsh.model.getBoundary(surfaces, oriented=False)
    for dim, tag in curves:
        gmsh.model.addPhysicalGroup(1, [abs(tag)], abs(tag))
        gmsh.model.setPhysicalName(1, abs(tag), f"boundary_{abs(tag)}")

    gmsh.model.mesh.generate(2)
    gmsh.model.mesh.optimize("Netgen")
    gmsh.write(output_msh)
    gmsh.finalize()
    print(f"Annular mesh written to {output_msh}")


def solve_poisson_annulus(msh_file: str = "annulus.msh") -> None:
    """
    Import gmsh mesh and solve Poisson equation on the annular domain.

    BC: u=1 on inner boundary, u=0 on outer boundary.
    """
    gmsh.initialize()
    gmsh.open(msh_file)

    domain, cell_tags, facet_tags = model_to_mesh(
        gmsh.model, MPI.COMM_WORLD, 0, gdim=2
    )
    gmsh.finalize()

    V = fem.functionspace(domain, ("Lagrange", 2))

    # Locate facets for BCs by geometric criterion
    def inner_boundary(x):
        return np.sqrt(x[0]**2 + x[1]**2) < 0.25

    def outer_boundary(x):
        return np.sqrt(x[0]**2 + x[1]**2) > 0.9

    dofs_inner = fem.locate_dofs_geometrical(V, inner_boundary)
    dofs_outer = fem.locate_dofs_geometrical(V, outer_boundary)

    u_inner = fem.Function(V); u_inner.x.array[:] = 1.0
    u_outer = fem.Function(V); u_outer.x.array[:] = 0.0

    bcs = [
        fem.dirichletbc(u_inner, dofs_inner),
        fem.dirichletbc(u_outer, dofs_outer),
    ]

    u = ufl.TrialFunction(V)
    v = ufl.TestFunction(V)
    a = ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx
    L = fem.Constant(domain, 0.0) * v * ufl.dx  # Laplace: f=0

    problem = LinearProblem(a, L, bcs=bcs)
    uh = problem.solve()
    uh.name = "u"

    with io.XDMFFile(MPI.COMM_WORLD, "annulus_solution.xdmf", "w") as xdmf:
        xdmf.write_mesh(domain)
        xdmf.write_function(uh)
    print("Annulus Laplace solution saved to annulus_solution.xdmf")

Advanced Usage

Steady Stokes Flow (Velocity-Pressure)

"""
Solve the steady Stokes equations (viscous flow at Re→0) on a channel.

−μ ∇²u + ∇p = f
∇·u = 0

Uses Taylor-Hood P2/P1 elements (LBB-stable mixed formulation).
"""

from mpi4py import MPI
import numpy as np
from dolfinx import mesh, fem
from dolfinx.fem.petsc import LinearProblem
import ufl


def solve_stokes_channel(
    nx: int = 64,
    ny: int = 16,
    mu: float = 1.0,      # dynamic viscosity
    U_max: float = 1.0,   # max inlet velocity
) -> None:
    """
    Solve steady Stokes flow in a 2D channel.
    Parabolic inlet profile, no-slip walls, stress-free outlet.
    """
    domain = mesh.create_rectangle(
        MPI.COMM_WORLD,
        [np.array([0.0, 0.0]), np.array([4.0, 1.0])],
        [nx, ny],
        cell_type=mesh.CellType.triangle,
    )

    # Taylor-Hood P2/P1 elements
    P2 = fem.functionspace(domain, ("Lagrange", 2, (2,)))
    P1 = fem.functionspace(domain, ("Lagrange", 1))

    # Mixed space
    V_el = ufl.VectorElement("Lagrange", domain.ufl_cell(), 2)
    Q_el = ufl.FiniteElement("Lagrange", domain.ufl_cell(), 1)
    W = fem.functionspace(domain, ufl.MixedElement([V_el, Q_el]))

    (u, p) = ufl.TrialFunctions(W)
    (v, q) = ufl.TestFunctions(W)

    f = fem.Constant(domain, np.array([0.0, 0.0]))

    a = (mu * ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx
         - ufl.div(v) * p * ufl.dx
         + q * ufl.div(u) * ufl.dx)
    L = ufl.inner(f, v) * ufl.dx

    # No-slip on top and bottom walls
    def walls(x):
        return np.isclose(x[1], 0.0) | np.isclose(x[1], 1.0)

    # Parabolic inlet profile: u_x = 4*U_max*y*(1-y), u_y=0
    def inlet_velocity(x):
        vals = np.zeros((2, x.shape[1]))
        vals[0] = 4.0 * U_max * x[1] * (1.0 - x[1])
        return vals

    def inlet(x):
        return np.isclose(x[0], 0.0)

    W0, _ = W.sub(0).collapse()
    dofs_walls = fem.locate_dofs_geometrical((W.sub(0), W0), walls)
    dofs_inlet = fem.locate_dofs_geometrical((W.sub(0), W0), inlet)

    u_no_slip = fem.Function(W0); u_no_slip.x.array[:] = 0.0
    u_inflow  = fem.Function(W0); u_inflow.interpolate(inlet_velocity)

    bcs = [
        fem.dirichletbc(u_no_slip, dofs_walls,  W.sub(0)),
        fem.dirichletbc(u_inflow,  dofs_inlet,  W.sub(0)),
    ]

    problem = LinearProblem(a, L, bcs=bcs,
                            petsc_options={"ksp_type": "minres", "pc_type": "hypre"})
    wh = problem.solve()
    print("Stokes channel flow solved. Extract wh.sub(0) for velocity, wh.sub(1) for pressure.")

Time-Dependent Heat Equation

"""
Solve the unsteady heat equation with backward Euler time integration.
∂u/∂t − α∇²u = 0  (α = thermal diffusivity)
"""

from mpi4py import MPI
import numpy as np
from dolfinx import mesh, fem
from dolfinx.fem.petsc import LinearProblem
import ufl


def solve_heat_equation(
    nx: int = 40,
    alpha: float = 0.01,
    T_final: float = 1.0,
    dt: float = 0.01,
) -> None:
    """Solve 2D heat equation with Gaussian initial condition."""
    domain = mesh.create_unit_square(MPI.COMM_WORLD, nx, nx)
    V = fem.functionspace(domain, ("Lagrange", 1))

    u_n = fem.Function(V)
    x = fem.Expression(
        ufl.exp(-50.0 * ((ufl.SpatialCoordinate(domain)[0] - 0.5)**2
                         + (ufl.SpatialCoordinate(domain)[1] - 0.5)**2)),
        V.element.interpolation_points()
    )
    u_n.interpolate(x)

    u = ufl.TrialFunction(V)
    v = ufl.TestFunction(V)
    dt_const = fem.Constant(domain, dt)

    a = (u * v + dt_const * alpha * ufl.inner(ufl.grad(u), ufl.grad(v))) * ufl.dx
    L = u_n * v * ufl.dx

    problem = LinearProblem(a, L, bcs=[],
                            petsc_options={"ksp_type": "cg", "pc_type": "hypre"})

    t = 0.0
    n_steps = int(T_final / dt)
    for step in range(n_steps):
        t += dt
        uh = problem.solve()
        u_n.x.array[:] = uh.x.array[:]
        if step % 20 == 0:
            max_u = np.max(np.abs(uh.x.array))
            print(f"t={t:.3f}: max(u)={max_u:.4f}")
    print("Heat equation time integration complete.")

Troubleshooting

External Resources

• FEniCS/dolfinx documentation: https://docs.fenicsproject.org/
• FEniCS tutorial (Langtangen & Logg): https://fenicsproject.org/pub/tutorial/
• gmsh documentation: https://gmsh.info/doc/texinfo/gmsh.html
• ParaView visualization: https://www.paraview.org
• Brenner & Scott, "The Mathematical Theory of Finite Element Methods", 3rd ed.

Examples

Example 1 — Full Poisson Convergence Study

if __name__ == "__main__":
    print("=== Poisson Convergence Study ===")
    convergence_study()
    print()
    print("=== Poisson on Unit Square (n=64) ===")
    result = solve_poisson_unit_square(n_cells=64, degree=2)
    print(f"L2 error with P2 elements: {result['L2_error']:.2e}")

Example 2 — Cantilever Beam Under Gravity

if __name__ == "__main__":
    print("=== Linear Elasticity: Steel Cantilever Beam ===")
    solve_linear_elasticity(
        nx=80, ny=20,
        E=210e9,    # Steel Young's modulus
        nu=0.3,
        rho=7850.0,
        g=9.81,
        output_xdmf="steel_beam.xdmf",
    )
    print("Open steel_beam.xdmf in ParaView to visualize displacement and von Mises stress.")

Example 3 — gmsh Annulus Mesh and Laplace Solve

if __name__ == "__main__":
    print("=== gmsh Annular Mesh + Laplace Equation ===")
    create_annular_mesh(r_inner=0.2, r_outer=1.0, mesh_size=0.05)
    solve_poisson_annulus("annulus.msh")
    print("Open annulus_solution.xdmf in ParaView.")

Changelog