Finite Element Solver

A finite element method (FEM) solver written from scratch in Python. It meshes a domain, assembles the discrete system, and solves the Poisson, heat, wave, and linear and nonlinear elasticity equations in 2D and 3D. Boundary conditions can be Dirichlet, Neumann, or Robin, and are described geometrically so one specification survives a remesh. It also includes its own meshing (Delaunay, Ruppert's algorithm, red-green refinement), adaptive refinement driven by a posteriori error estimators, quadratic and curved elements, buckling and modal analysis, adjoint sensitivities, and SIMP topology optimization.

▶ See it running: the demo gallery

The gallery renders every demo beside the code that produced it and is rebuilt on each push to main. The figures below are a subset of it, with the full captions in the gallery.

---

Meshing a domain

The library includes its own meshing. An Outline of lines, arcs, circles, and Bezier curves, drawn by hand or traced from an SVG, is simplified with Douglas-Peucker and triangulated with Ruppert's algorithm (outline.mesh(min_angle=..., max_area_fraction=...)) to a minimum-angle and maximum-area bound, and can then be adaptively refined where the solution needs it (shown in later demos). Each boundary facet of the result is tagged with the outline it came from, so a hole can take a condition by on_tag(1) rather than by coordinates. Below, the same Poisson problem ($-\nabla^2 u = 1$ with $u = 0$ on the boundary) is solved on four outlines. Each asks something different of the mesher: disconnected islands, Bezier curves, a hole, and sharp notches.

Solving PDEs

Poisson's equation

Poisson's equation, $-\nabla^2 u = f$, models heat transfer, electrostatics, and other steady diffusion. The meshing section above solves it with a constant source; with no source it is Laplace's equation, which here gives a flow. An ideal (incompressible, irrotational) flow has a velocity potential $\phi$ with $\mathbf{v} = \nabla\phi$, and $\phi$ solves Laplace's equation. The obstacle is a NACA 2412 airfoil at a 12-degree angle of attack. A potential difference drives the flow left to right. The wing carries no boundary condition at all, which in the weak form is the natural zero-flux condition, so it becomes a streamline the flow parts around. The equipotentials crowd over the upper surface, where the flow speeds up.

Heat equation

The heat equation, $\partial u / \partial t = \alpha \nabla^2 u$, describes how temperature spreads over time. It is integrated with the theta-method, defaulting to $\theta = \tfrac{1}{2}$ (Crank-Nicolson); $\theta = 1$ is backward Euler. A finned heatsink is held hot underneath its base, and every other surface sheds heat to ambient through a convective (Robin) film, $\partial u / \partial n + \kappa (u - u_\infty) = 0$.

The finned sink is compared against a solid block of the same size. Driven by the same chip power, the block runs about 108°C above ambient while the finned sink runs 58°C, roughly halving the thermal resistance. Held at the same base temperature, the finned sink sheds about 1.8x the heat with two-thirds the metal. The fin efficiency (right), the heat a fin sheds relative to what it would shed at the base temperature throughout, follows the textbook $\tanh(mL)/(mL)$ law and falls as fins lengthen, since a long fin runs cold toward the tip.

Wave equation

The wave equation, $\partial^2 u / \partial t^2 = c^2 \nabla^2 u$, is second order in time and is integrated with Newmark's average-acceleration method. Below, a wave front meets a harbor breakwater: it reflects off the wall and passes the gap, where it spreads into the sheltered water as a circular wave centred on the opening.

Solids & structures

Linear elasticity, in 2D and 3D

The linear elastic solver computes the displacement and the full stress tensor from applied loads and boundary conditions. A cantilever is clamped on the left and pulled down over the middle of the right edge. The bending stress is largest at the clamp, with tension above the neutral axis and compression below. The 3D panel is a tetrahedral cantilever under the same clamp and load, drawn as its boundary surface. It is solved with AMG-preconditioned conjugate gradients, since in 3D a direct factorization's fill-in starts to hurt.

One solve gives one stress tensor, and the four 2D stress panels are rotation-invariant reductions of it: von Mises, mean normal stress, the Tresca measure, and the largest tensile principal value.

Stress at a re-entrant corner, and why fillets exist

An L-bracket clamped at the top and pulled down at the tip concentrates stress at its inner corner. A sharp re-entrant corner is a stress singularity, where the exact elastic stress is infinite, so no mesh resolves it and adaptive refinement into the corner keeps the computed peak climbing. A fillet removes the singularity and the peak settles on a finite value. The plot on the right tracks the corner peak against mesh size for both. The sharp corner climbs without bound, so the stress it reports is a property of the mesh rather than the part, while the fillet converges. This is why real parts round their inner corners.

Linear and nonlinear materials, one stretch

The same clamped block is stretched to the right in four ways, each panel coloured by how hard the material is working inside. The first two are ordinary small-strain elasticity, the familiar model where stress grows in proportion to strain, solved two different ways: directly, and by minimising the block's stored elastic energy. They settle into the same shape to machine precision, but their colour differs because they report stress differently: the first uses engineering stress, force over the block's original cross-section, the second the true stress on the stretched, now-thinner cross-section, which reads higher. The last two switch to nonlinear elasticity, which large deformations actually need. St-Venant-Kirchhoff is the simplest such model but stiffens far too aggressively in tension, while Neo-Hookean, a rubber-like law, stays realistic out to large stretch.

All four panels share one colour scale. The stresses span a wide range, so the scale is logarithmic: brighter means more stress, and a whole panel shifting brighter (as St-Venant-Kirchhoff's does) means that material is carrying more everywhere.

From an outline to a stress concentration

This demo runs the whole pipeline. A plate with a hole is meshed from its outline, given roller and traction conditions (the rim is left traction-free), solved on curved quadratic elements, and adaptively refined toward the stress at the rim. The stress crowds into the material either side of the hole and relaxes to the applied value within about a diameter. At the rim it peaks at 3.03x the applied stress. The classic Kirsch factor is 3 for a hole in an infinite plate, and Howland's value for a hole a tenth of this plate's width is 3.02.

Buckling analysis

A slender column under compression does not fail by crushing. It snaps sideways once the load crosses a critical value. BucklingAnalysis finds that value and the buckled shapes by linearised (eigenvalue) buckling. A reference load sets up a prestress, the geometric stiffness $K_g$ is assembled from it, and the generalized eigenproblem $K \phi = -\lambda K_g \phi$ gives the critical load factors and mode shapes. The column is meshed with P2 elements, which do not lock in bending the way a constant-strain triangle does.

Modal (free-vibration) analysis

Free vibration solves $K \phi = \omega^2 M \phi$ with the consistent mass matrix, using shift-invert about zero to find the lowest frequencies. No load is applied, so the modes are a property of the structure alone. A steel tuning fork is meshed from its outline and held at the stem base. Its low modes come in pairs, and the one whose tines swing in opposite directions leaves the stem still and rings; that is the fork's voice.

Topology optimization

Topology optimization distributes material to minimize compliance (deformation under load). Here a simply supported beam carries a central load, and the SIMP (Solid Isotropic Material with Penalization) method finds the stiffest structure using half the material, penalizing intermediate densities so the design resolves toward solid or void. It finds the classic arch, a compression arch over a tension tie braced by a diagonal web. Compliance is the work the load does, so it measures deflection directly, and the optimized truss is only about 1.7x as compliant as the fully solid block on half the material. What it removed was near the neutral axis, where the material was barely resisting the bending.

The gallery's topology page plays the SIMP iterations frame by frame, from an even grey to the black-and-white truss.

Accuracy & performance

Convergence against manufactured solutions

Against exactly known (manufactured) solutions, P1 elements are second order in space (halve $h$, quarter the error) for both a scalar unknown and a coupled vector one, and P2 elements third. In time, backward Euler is first order and Crank-Nicolson second. How the load is built matters too: sampling the source at the quadrature points instead of reading it at the vertices keeps the rate and improves the constant about 3x. Every rate here also runs as an assertion in the test suite.

Higher-order elements

P2 (quadratic) triangles carry edge-midpoint DOFs that let the solution curve within an element. On the same meshes they are third order in $L^2$ where P1 is second (top left above), and reach a given accuracy with fewer degrees of freedom (top right).

Curved (isoparametric) boundary elements go a step further. On a boundary that carries an analytic curve (a Circle or Arc), an IsoparametricTriangleElement places its edge-midpoint node on the true curve, so the element's boundary edge follows the curve instead of cutting a chord. Meshing carries the curve through, so Ruppert's split points and red-green refinement project onto it and a circular hole stays round under refinement. The meshed area then converges at the element's own order rather than the polygonal $O(h^2)$; the hole-in-plate and bracket demos above run on these elements.

Adaptive refinement

Adaptive refinement re-solves and splits wherever an a posteriori error estimator finds the most error, keeping triangle quality with red-green refinement. There are three estimators. The residual estimator measures the interior residual, the flux jump across edges, and the boundary residual (2D only). The Zienkiewicz-Zhu recovery estimator measures the gap between the discrete flux and a recovered continuous one. The goal-oriented estimator refines toward a chosen quantity of interest through an adjoint solve. Below, the residual estimator on a peaked Poisson source concentrates the mesh where the solution is hardest to approximate, and reaches a given error with about a third of the unknowns uniform refinement needs (the gallery has the chart).

Representation error

Before any PDE is solved, the mesh already limits what its P1 space can represent. The target $\sin(40 r^2)$ has rings that tighten with radius. Projected onto a coarse P1 mesh, the slow inner rings come through but the fast outer ones break up into the triangulation. This representation error is the floor every solver on this mesh starts from, and refining the mesh is what lowers it.

---

Quick Start

from fem import box_mesh, Conditions, Dirichlet, Poisson, Plotter, Source
from fem.regions import everywhere

mesh = box_mesh(corners=[[0, 0], [1, 1]], resolution=(40, 40))

# Conditions are geometric, so the same ones are valid on any mesh of this domain:
# what is applied to the domain (supports, loads, a source) all in one object.
conditions = Conditions(Dirichlet(everywhere(), 0), Source(1))

solution = Poisson().problem(mesh, conditions).solve()

plotter = Plotter(title="Poisson")
plotter.plot(mesh, solution, mode="surface")
plotter.show()

The equation builds a problem from the parts, which compose directly. The same solve, by hand:

from fem import FunctionSpace, LinearProblem, DiffusionForm

space = FunctionSpace(mesh, n_components=1)
problem = LinearProblem(space, DiffusionForm(), conditions)
solution = problem.solve()

problem.solve() picks LinearSolve for a constant tangent and Newton otherwise; strategy= chooses how it iterates, and the problem's backend (problem.with_backend(IterativeBackend())) how each linear system is solved, independently. The result is typed by the physics: Poisson().problem(mesh, conditions).solve() is a DiffusionSolution, an elastic one an ElasticSolution, with no narrowing at the call.

What you choose at each step

Every step of a solve is one choice among a few named objects, all importable from fem. Steps run top to bottom; a step with several parts lists each, and a part marked optional can be skipped. The last column is where each lives. ARCHITECTURE.md explains how they fit; this is the menu.

Step	Part	Options	Where
Outline (optional)		an Outline is closed loops of Line, Arc, Circle, and CubicBezier pieces. It is usually drawn from point lists with Outline.from_polygons(...), or traced from a drawing with Outline.from_svg(...) and thinned with .simplified(tolerance).	fem/mesh/outline.py, fem/mesh/curves.py, fem/mesh/svg.py
Mesh	build	a mesh is generated by triangulating an outline with outline.mesh(min_angle=, max_area_fraction=), or with a structured helper like box_mesh(corners, resolution); arrays from elsewhere are wrapped as Mesh(vertices, elements)	fem/mesh/mesh.py, fem/mesh/structured.py, fem/mesh/pslg.py, fem/mesh/ruppert.py
	refine (optional)	a mesh is refined by hand with RedGreenRefiner, or where the error is by AdaptiveRefinement (the Outer loop below)	fem/mesh/refinement.py
Equation		the PDE and its constants: Projection, Poisson, Heat, Wave, LinearElastic, or FiniteStrainElastic, which also takes a law= of StVenantKirchhoff (the default) or NeohookeanEnergyDensity; LinearElastic also takes a thermal= of ThermalStrain for thermoelasticity, its temperature a constant, a callable, or a Poisson / Heat solution on the same mesh, and a reduction= of plane strain (the default) or plane stress, what a 2D solve means in the third direction	fem/physics/equations.py, fem/physics/forms.py, fem/physics/materials.py, fem/physics/energies.py
Conditions	where	a region of the domain: everywhere, on_plane, in_box, or on_tag for an outline loop, combined with union and intersect; at_indices names nodes directly and does not survive a remesh	fem/regions.py
	what	boundary conditions Dirichlet, Neumann, Robin on a region, a volume Source, PointLoads, and the Initial(u0, v0=) state a solve starts from, collected into one Conditions(...); with no Initial the solve starts at rest	fem/conditions.py, fem/boundary.py, fem/loads.py
	value	each takes a constant, a callable of position, or a TimeDependent callable of position and time	fem/regions.py
Element		the element type sets the shape functions (linear or quadratic), whether the geometry is straight-sided or curved, and the quadrature degree: LinearTriangleElement, LinearTetrahedralElement, QuadraticTriangleElement, QuadraticTetrahedralElement, or IsoparametricTriangleElement; the default is the linear element of the mesh's dimension	fem/elements.py, fem/quadrature.py
Problem		equation.problem(mesh, conditions, element_type=) discretizes the equation on the mesh and resolves the conditions against it	fem/problem.py, fem/space.py
Solve	which solve	a steady solve is problem.solve(); a solve in time is ThetaMethod(dt=, steps=) or NewmarkMethod(dt=, steps=) (with RayleighDamping) .solve(problem); an eigen-solve is BucklingAnalysis or ModalAnalysis .solve(problem); any takes initial= to start from an Initial other than the conditions' own	fem/algebra/integrators.py, fem/analysis/buckling.py, fem/analysis/modal.py
	strategy	how the problem is iterated, via strategy=: LinearSolve, or NewtonSolve with a BacktrackingLineSearch and TangentRegularization; by default, LinearSolve for a constant tangent and NewtonSolve otherwise	fem/algebra/solve.py
	backend	how each linear system is solved, held by the problem (problem.with_backend(IterativeBackend())): DirectBackend (the default), IterativeBackend (the default above a size threshold when the operator is SPD), or MinresBackend; a linear problem factors once and holds the result for every solve	fem/algebra/backends.py, fem/algebra/system.py
Outer loop (optional)		a driver that re-solves: AdaptiveRefinement refines where a ResidualEstimator, RecoveryEstimator, or GoalOrientedEstimator finds error; DesignOptimizer moves the density of a SIMPModel; both use SensitivityAnalysis for gradients	fem/analysis/adaptivity.py, fem/analysis/estimators.py, fem/analysis/design.py, fem/analysis/sensitivity.py
Result		a steady solve returns a DiffusionSolution, ElasticSolution, or FieldSolution, each a NodalField; a solve in time a TransientSolution, an eigen-solve a BucklingSolution or ModalSolution, and the optimizer a DesignHistory	fem/post/solution.py, fem/post/invariants.py, fem/post/recovery.py, fem/post/io.py
Plot		Plotter.plot(target, values, mode=) draws a solution, field, or mesh in one of the modes mesh, boundary, colored, surface, arrows, solid, bc, refinement	fem/plot/

A solution is a typed dataclass. An elastic solve returns an ElasticSolution, which carries the stress and strain as full tensors and derives the scalar measures on demand:

solution = LinearElastic(E=200, nu=0.3).problem(mesh, conditions).solve()

solution.dofs              # (n_nodes * n_components,) the DOF vector
solution.nodal_values      # (n_nodes, n_components) the same by node
solution.evaluate(points)  # (n_points, n_components) the field at any points
solution.deformed_mesh()   # the mesh displaced by the field
solution.stress            # (n_elements, 3, 3) Cauchy stress tensors
solution.von_mises         # (n_elements,) equivalent stress, the usual plot
solution.principal_stress  # (n_elements, 3) principal values, ascending
solution.compliance        # (n_elements,) strain energy per element

fem/post/invariants.py holds those reductions; each is rotation-invariant.

Installation

The project uses uv for environment and dependency management. uv sync creates a project-local .venv, installs the fem package in editable mode, and pins exact versions in uv.lock.

uv sync   # core solver, SVG-outline and 3D tetrahedral meshing, and dev tools (pytest)

Prefer plain pip? It is a standard pyproject.toml package:

pip install -e .

Running demos and tests

Runnable demos live in examples/demos/ (run from the repo root) behind a small CLI. Each demo is a package of two files: physics.py poses and solves the problem and is what the gallery shows as the demo's source, figures.py draws it.

uv run python examples/cli.py list          # see every available demo
uv run python examples/cli.py run poisson   # run one by name
uv run python examples/cli.py gallery       # render them all as a browsable site

The test suite:

uv run pytest

uv run ruff check and uv run pyright gate CI alongside the tests.

The figures in this README are committed images (PNGs and one GIF), refreshed by hand. After changing a demo, regenerate them and commit the result:

uv run python examples/make_readme_figures.py   # rewrites the figures in images/

Methods

• Galerkin finite element method: P1 (linear) and P2 (quadratic) bases on triangles and tetrahedra, plus curved isoparametric triangles, over a Gaussian quadrature layer
• Boundary conditions: Dirichlet, Neumann, Robin, and per-component (roller) constraints
• PDEs: L2 projection, Poisson, variable-coefficient diffusion, heat, wave, and two elastic models: LinearElastic (infinitesimal strain, Navier-Cauchy) and FiniteStrainElastic (geometrically exact Green-Lagrange strain under a hyperelastic law, St Venant-Kirchhoff by default); 2D elasticity is plane strain by default or plane stress by choice; one-way thermoelasticity, a temperature field driving a thermal strain in the small-strain model
• Time integration: theta-method (backward Euler, Crank-Nicolson) for first-order systems; Newmark average-acceleration for second-order
• Linear algebra: sparse throughout; direct (splu) by default, or AMG-preconditioned CG for large SPD systems, with a rigid-body near-kernel for elasticity
• Eigen-analysis: linearised (eigenvalue) buckling and free-vibration (modal) analysis, a geometric-stiffness or mass pencil and a sparse generalized eigensolve
• Derived fields: Cauchy stress and strain tensors, von Mises, principal stresses, compliance
• Error estimation: residual, Zienkiewicz-Zhu recovery, and goal-oriented (adjoint-weighted) estimators, driving closed-loop adaptive refinement
• Sensitivity and optimization: adjoint gradients of a quantity of interest; Newton-Raphson with an optional backtracking line search; optimality criteria (SIMP topology and design optimization)
• Mesh algorithms: Delaunay triangulation, Ruppert's algorithm (line segments to triangle mesh), red-green refinement

Roadmap

BACKLOG.md tracks the detailed open work; this is the direction.

• Broaden the physics: plasticity, fluids (Stokes / Navier-Stokes), electrostatics, advection-diffusion
• Inverse problems and shape optimization on the adjoint core
• Standard formats: STL and OBJ meshes, Gmsh .msh import, VTK/ParaView .vtu export
• Standard benchmark suite: NAFEMS, Cook's membrane, plate-with-hole, L-shaped singularity, Euler columns
• Finish the core: P2-aware plotting in 3D and P2 adaptivity, mixed (u-p) for near-incompressibility, time-varying loads, two-grid preconditioner

References

The Finite Element Method: Theory, Implementation, and Applications by Mats G. Larson and Fredrik Bengzon.

SIMP Method for Topology Optimization by Dassault Systèmes.