Evolving FEniCS: The Extension Ecosystem at Simula Scientific Computing

FEniCS 2026 at University of Chicago in Paris
Jørgen S. Dokken
Henrik N.T. Finsberg

It all started over 21 years ago

Packages developed or maintained

• UFL/FFC(x)/DOLFIN(x)
• DOLFINx_MPC
• scifem
• io4dolfinx (adios4dolfinx)
• fenicsx_ii
• DOLFIN(x)-adjoint
• networks_fenicsx

Scifem - FEM prototyping playground

Not all ideas are good ideas in the beginning

  r_el = basix.ufl.real_element(mesh.basix_cell(), shape=(2, 3))
R = dolfinx.fem.functionspace(mesh, r_el)

  dolfinx.fem.petsc.NonlinearProblem

  dolfinx.mesh.transfer_meshtags_to_submesh

What is next?

mesh = dolfinx.mesh.create_unit_square(comm, 10, 10)
tdim = mesh.topology.dim
tol = 1e-14

# Divide cell into three regions
tags = (4,5,8)
cell_map = mesh.topology.index_map(tdim)
num_cells_local = cell_map.size_local + cell_map.num_ghosts
markers = np.full(num_cells_local, tags[0],  dtype=np.int32)
markers[dolfinx.mesh.locate_entities(
    mesh, tdim, lambda x: x[0] <= 0.5+tol)] = tags[1]
markers[dolfinx.mesh.locate_entities(
    mesh, tdim, lambda x: x[1] <= 0.5+tol)] = tags[2]
cells = np.arange(num_cells_local, dtype=np.int32)
ct = dolfinx.mesh.meshtags(mesh, tdim, cells, markers)

# Create a piecewise constant (per region) function space
V = create_space_of_simple_functions(mesh, ct, tags)
u = dolfinx.fem.Function(V)
u.x.array[0] = 3.2
u.x.array[1] = 5.5
u.x.array[2] = 4.2
assert len(u.x.array) == 3

What is next?

from scifem import closest_point_projection
points, ref_points = closest_point_projection(
    mesh,
    closest_cells,
    points,
    tol_x=1e-7,
)

Based on simplex projections^2,3,4

IO4DOLFINx - a unified IO?

Visualization and checkpointing (write/read functions) has diverged due to N+1 different file formats and finite elements.

IO4DOLFINx - a unified IO?

ADIOS4DOLFINx⁶ introduced a specific split between readable and visualizable functions.

IO4DOLFINx - a unified IO?

Reality is that most users use iso-parameteric finite elements (often P1).

IO4DOLFINx is a backend agnositic interface to many mesh formats.

• gmsh, PyVista, XDMF, VTKHDF
  • {read/write}_{point/mesh}_data
• adios2, h5py
  • {read/write}_checkpoint

FEniCSx_ii

Example based on⁸

$$\begin{align*} - \nabla \cdot (\alpha_1 \nabla u) + \xi (\Pi_R(u) - p))\delta_\Gamma &= f && \text{in } \Omega, \\ - d_s(A d_s p) + P \xi (p - \Pi(u)) &= A \hat f && \text{in } \Lambda, \\ u&=g &&\text{on } \partial\Omega,\\ A d_s p &=0 && \text{at } s\in\{0, 1\}.\\ \int_\Omega \alpha_1 \nabla u \cdot \nabla v~\mathrm{d}x + \int_\Gamma P\xi (\Pi_R(u) - p)\Pi_R(v)~\mathrm{d}s &= \int_\Omega f\cdot v~\mathrm{d}x\\ \int_\Gamma A d_s p \cdot d_s q~\mathrm{d}s + \int_\Gamma P\xi (p - \Pi_R(u))q~\mathrm{d}s &= \int_\Gamma A\hat f\cdot q~\mathrm{d}s \end{align*}$$

Re-implementation of FEniCS_ii⁹

from fenicsx_ii import Average, Circle, LinearProblem, assemble_scalar
V = dolfinx.fem.functionspace(omega, ("Lagrange", 1))
Q = dolfinx.fem.functionspace(lmbda, ("Lagrange", 1))
W = ufl.MixedFunctionSpace(*[V, Q])

R, q_degree = 0.05, 20
restriction_trial = Circle(lmbda, R, degree=q_degree)
restriction_test = Circle(lmbda, R, degree=q_degree)

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

q_el = basix.ufl.quadrature_element(lmbda.basix_cell(), value_shape=(), degree=q_degree)
Rs = dolfinx.fem.functionspace(lmbda, q_el)
avg_u = Average(u, restriction_trial, Rs)
avg_v = Average(v, restriction_test, Rs)

Uses intermediate non-matching
interpolation matrices as⁹

dx_3D = ufl.Measure("dx", domain=omega)
dx_1D = ufl.Measure("dx", domain=lmbda)

A = ufl.pi * R**2
P = 2 * ufl.pi * R
xi = dolfinx.fem.Constant(omega, 1.0)
x = ufl.SpatialCoordinate(omega)
a = ufl.inner(ufl.grad(u), ufl.grad(v)) * dx_3D
a += P * xi * ufl.inner(avg_u - p, avg_v) * dx_1D
a += A * ufl.inner(ufl.grad(p), ufl.grad(q)) * dx_1D
a += P * xi * ufl.inner(p - avg_u, q) * dx_1D
L = f_vol * v * dx_3D
L += f_line * q * dx_1D

Networks_FEniCSx

MPI compatible FEniCSx+Networkx based on ^10,11

from networks_fenicsx import HydraulicNetworkAssembler, NetworkMesh, Solver
from networks_fenicsx.network_generation import make_arterial_tree
from networks_fenicsx.post_processing import export_functions, extract_global_flux
n = 5
G = make_arterial_tree(N=n, direction=np.array([0.1, 1, 0]))
network_mesh = NetworkMesh(
    G, N=40, color_strategy=nx.coloring.strategy_largest_first)
assembler = HydraulicNetworkAssembler(
    network_mesh, flux_degree=1, pressure_degree=0)
assembler.compute_forms(p_bc_ex=p_bc_expr)
solver = Solver(assembler, kind="nest")
solver.assemble()
sol = solver.solve()
global_flux = extract_global_flux(network_mesh, sol)

DOLFINx-adjoint

u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
F = ufl.inner(kappa * ufl.grad(u), ufl.grad(v)) * ufl.dx - f * v * ufl.dx
a, L = ufl.system(F)
uh = dolfinx_adjoint.Function(V, name="State")
petsc_options = {
    "ksp_type": "preonly",
    "pc_type": "lu",
    "pc_factor_mat_solver_type": "mumps",
    "ksp_error_if_not_converged": True,
}
problem = dolfinx_adjoint.LinearProblem(
    a,
    L,
    u=uh,
    bcs=[bc],
    petsc_options=petsc_options,
    adjoint_petsc_options=petsc_options,
    tlm_petsc_options=petsc_options,  # type: ignore
)
problem.solve()

DOLFINx-adjoint

J_symbolic = 0.5 * ufl.inner(uh - d, uh - d) * ufl.dx
J_symbolic += 0.5 * alpha * ufl.inner(f, f) * ufl.dx
J = dolfinx_adjoint.assemble_scalar(J_symbolic)

control = pyadjoint.Control(f)
Jhat = pyadjoint.ReducedFunctional(J, control)

optimization_problem = pyadjoint.MoolaOptimizationProblem(Jhat)
f_moola = DolfinxPrimalVector(f)
optimization_opts = {
  "jtol": 0, "gtol": 1e-9,
  "Hinit": "default", "maxiter": 100,
  "mem_lim": 10, "rjtol": 0}
solver = moola.BFGS(optimization_problem, f_moola, options=optimization_opts)
solution = solver.solve()

Irksome - time derivatives in UFL^11,12,13

from irksome import Dt, MeshConstant
el_u = basix.ufl.element("Lagrange", ct, 3, shape=(gdim,))
el_p = basix.ufl.element("Lagrange", ct, 2)
W = dolfinx.fem.functionspace(msh, basix.ufl.mixed_element([el_u, el_p]))

MC = MeshConstant(msh, backend="dolfinx")
t, dt = MC.Constant(0.0), MC.Constant(1.0 / N)
z = dolfinx.fem.Function(W)
u, p = split(z)
(v, q) = TestFunctions(W)
F = inner(Dt(u), v) * dx + inner(grad(u), grad(v)) * dx \
  - inner(p, div(v)) * dx - inner(div(u), q) * dx \
  - inner(f, v) * dx

Dirichlet BCs with UFL-expressions

from irksome.backends.dolfinx import dirichletbc
bc = dirichletbc(bc_u_as_ufl_expr, boundary_dofs, W.sub(0))
bc_p = dirichletbc(bc_p_as_ufl_expr, corner_dof, W.sub(1))
bcs = [bc, bc_p]

We can now choose the accuracy of the
time-derivative

from irksome import GaussLegendre
from irksome.stage_derivative import StageDerivativeTimeStepper
from irksome.tools import AI

butcher_tableau = GaussLegendre(num_stages=num_stages)
linear_stepper = StageDerivativeTimeStepper(
        F, butcher_tableau, t, dt, z,
        bcs=bcs,
        Fp=None,
        bc_type="DAE",
        splitting=AI,
        solver_parameters=solver_parameters,
        backend="dolfinx",
    )
linear_stepper.advance()

Thanks to all my collaborators

Simula 25 years - FEniCS workshop

Hybrid workshop September 8th - 9th

Confirmed invited speakers

Antonio B. Svizzero (Undabit)
Cécile Daversin-Catty (SRL)
Chris Richardson (Cantab)
David Ham (IC)
Francesco Ballarin (UNICATT)
Henrik N. T. Finsberg (SRL)
Hyunsun Alicia Kim (UCSD)
Jack S. Hale (UNI.LU.)
Jeremy Bleyer (ENPC)
Joakim Sundnes (SRL)

Johan Hoffman (KTH)
Kent Andre Mardal (SRL/UIO)
Martin Řehoř (Rafinex)
Neeraj Cherukunnath (Rolls Royce)
Padmini Rangamani (UCSD)
Remi Delaporte-Mathurin (MIT)
Robert Kirby (BU)
Simon W. Funke (formely SRL)
Susanne Claus (ONERA)
Thomas M. Surowiec (SRL)

Organizers: Jørgen S. Dokken, Cécile Daversin-Catty, Ada Johanne Ellingsrud, Eirik Valseth