From fd3cc8b3b8dfd8c15c667ef76a6d81210f551bb9 Mon Sep 17 00:00:00 2001
From: Themis Skamagkis <themisskamagkis@gmail.com>
Date: Mon, 12 Oct 2026 09:49:17 +0200
Subject: [PATCH 01/55] [all] Ignore bytecode and verification results

---
 .gitignore | 5 +++++
 1 file changed, 5 insertions(+)
 create mode 100644 .gitignore

diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..e86d952
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,5 @@
+# ignore cached bytecode-compiled python scene files
+*__pycache__/
+
+# a run's own output
+verification/results/

From b092d84558518cbad11193e409c78f812acbc82b Mon Sep 17 00:00:00 2001
From: Themis Skamagkis <themisskamagkis@gmail.com>
Date: Mon, 12 Oct 2026 09:49:17 +0200
Subject: [PATCH 02/55] [docs] Introduce the verification suite in the README

---
 README.md | 33 +++++++++++++++++++++++++++++++++
 1 file changed, 33 insertions(+)

diff --git a/README.md b/README.md
index 25d0262..843a7da 100644
--- a/README.md
+++ b/README.md
@@ -1,2 +1,35 @@
 # Verification-Validation
 Initially funded by TIRREX, this repository proposes a framework for verification and validation. It relies on both manufactured solutions and comparison with alternative simulation technologies. This work therefore proves the reliability and compliance of the SOFA models &amp; algorithms.
+
+For how the repository is organised, how the suite works and how to extend it, see
+[docs/architecture.md](docs/architecture.md).
+
+## Verification
+
+**Verification** uses the method of manufactured solutions: an exact displacement is chosen, the
+loads that produce it are derived, and SOFA solves the problem on a sequence of refined meshes.
+The error must shrink at the rate the discretisation predicts. Covered force fields:
+- `LinearSmallStrainFEMForceField`
+
+## Layout
+
+```
+common/                                  shared code: geometry and SOFA helpers
+verification/                            the verification suite
+  cases/<Component>/**.json              one test case per file
+```
+
+## Requirements
+
+- SOFA with SofaPython3
+- numpy, sympy, matplotlib
+
+## Running
+
+```bash
+python verification/run.py <path-to-json-file>  # one test case
+python verification/run.py --all                # every test case
+python verification/run.py --clean              # delete previous results
+```
+
+Results are written to `verification/results/`.

From 46811280ebd7a664a2389cc3d8f876bae19c6c6f Mon Sep 17 00:00:00 2001
From: Themis Skamagkis <themisskamagkis@gmail.com>
Date: Mon, 12 Oct 2026 09:49:17 +0200
Subject: [PATCH 03/55] [common] Add the common package

---
 common/__init__.py | 0
 1 file changed, 0 insertions(+), 0 deletions(-)
 create mode 100644 common/__init__.py

diff --git a/common/__init__.py b/common/__init__.py
new file mode 100644
index 0000000..e69de29

From 36d1afb5defdf9b761d69601094072031fac1f63 Mon Sep 17 00:00:00 2001
From: Themis Skamagkis <themisskamagkis@gmail.com>
Date: Mon, 12 Oct 2026 09:49:17 +0200
Subject: [PATCH 04/55] [geometry] Add the tool-agnostic geometry base

---
 common/geometry/geometry.py | 116 ++++++++++++++++++++++++++++++++++++
 1 file changed, 116 insertions(+)
 create mode 100644 common/geometry/geometry.py

diff --git a/common/geometry/geometry.py b/common/geometry/geometry.py
new file mode 100644
index 0000000..46e44ef
--- /dev/null
+++ b/common/geometry/geometry.py
@@ -0,0 +1,116 @@
+"""Tool-agnostic geometry descriptions for the VnV suite."""
+
+from abc import ABC, abstractmethod
+
+
+def _region_indices_from_mesh(mesh_path):
+    """Node indices tagged by each physical group of a Gmsh .msh file, keyed by group name."""
+    import meshio
+
+    mesh = meshio.read(mesh_path)
+    tag_names = {tag: name for name, (tag, _dim) in mesh.field_data.items()}
+    groups = {}
+    for block, tags in zip(mesh.cells, mesh.cell_data.get("gmsh:physical", [])):
+        for element, tag in zip(block.data, tags):
+            name = tag_names.get(tag)
+            if name is not None:
+                groups.setdefault(name, set()).update(element.tolist())
+    return {name: sorted(indices) for name, indices in groups.items()}
+
+
+class Geometry(ABC):
+    """A shape parameterization plus its named boundary regions."""
+
+    dim: int  # topological dimension of the shape (and of the elements meshing it)
+    PARAMETER_ORDER: tuple  # the family's parameters, ordered by the dimension that introduces each
+
+    def __init__(self, *, dim=None, spatialDimensions=None, meshPath=None, **parameters):
+        """Build this shape, taking each name in its `PARAMETER_NAMES` as a keyword.
+
+        e.g. ``Beam3D(length=2.0, width=1.0, height=1.0)``. Pass `dim` only for a family base,
+        which names no dimension of its own: ``Grid(dim=3, length=2.0, width=1.0, height=1.0)``
+        builds that same box.
+        """
+        stated = getattr(type(self), "dim", None)
+        if stated is None and dim is None:
+            raise TypeError(f"{type(self).__name__} states no dim: pass dim= to build one, or "
+                            f"name the geometry that states it")
+        if stated is not None and dim is not None and dim != stated:
+            raise TypeError(f"{type(self).__name__} is {stated}D and cannot be built with dim={dim}")
+        self.dim = stated if dim is None else dim
+
+        if set(parameters) != set(self.PARAMETER_NAMES):
+            raise TypeError(f"{type(self).__name__} takes {list(self.PARAMETER_NAMES)}, "
+                            f"got {sorted(parameters)}")
+
+        for name in self.PARAMETER_NAMES:
+            setattr(self, name, parameters[name])
+        self._set_spatial_dimensions(spatialDimensions)
+
+        if meshPath is not None:
+            self._regions = {name: (lambda nodes, indices=indices: indices)
+                             for name, indices in _region_indices_from_mesh(meshPath).items()}
+        else:
+            self._regions = {name: (lambda nodes, predicate=predicate:
+                                    [i for i, p in enumerate(nodes) if predicate(p)])
+                             for name, predicate in self._analytic_predicates.items()}
+
+    @property
+    def PARAMETER_NAMES(self) -> tuple:
+        """The parameters this shape takes: the first `dim` of its family's.
+
+        Each dimension adds one, so a member takes as many as it has dimensions: of Grid's
+        ("length", "width", "height"), Bar1D takes ("length",) and Beam2D ("length", "width").
+        A family that does not grow this way sets `PARAMETER_NAMES` itself instead.
+        """
+        return self.PARAMETER_ORDER[:self.dim]
+
+    @property
+    def parameters(self) -> list:
+        """This shape's defining parameters, read off the attributes named in `PARAMETER_NAMES`.
+
+        e.g. Beam2D(length=2.0, width=0.5).parameters == [2.0, 0.5]
+        """
+        return [getattr(self, name) for name in self.PARAMETER_NAMES]
+
+    @property
+    def named_parameters(self) -> dict:
+        """Parameters keyed by name: how a manufactured solution asks for one of them.
+
+        e.g. Beam2D(length=2.0, width=0.5).named_parameters == {"length": 2.0, "width": 0.5}
+        """
+        return dict(zip(self.PARAMETER_NAMES, self.parameters))
+
+    def _set_spatial_dimensions(self, spatial_dimensions):
+        """Set how many dimensions of space this shape is embedded in; defaults to its own."""
+        if spatial_dimensions is None:
+            spatial_dimensions = self.dim
+
+        if spatial_dimensions < self.dim:
+            raise ValueError(
+                f"{type(self).__name__} is {self.dim}D and cannot be embedded in a "
+                f"{spatial_dimensions}D space: the embedding space must be at least as "
+                f"large as the shape."
+            )
+        if spatial_dimensions > 3:
+            raise ValueError(
+                f"{type(self).__name__}: spatial_dimensions cannot exceed 3, "
+                f"got {spatial_dimensions}."
+            )
+
+        self.spatial_dimensions = spatial_dimensions
+
+    @property
+    def region_names(self):
+        """This geometry's named boundary regions."""
+        return self._regions.keys()
+
+    def region_indices(self, region, nodes) -> list:
+        """Indices into `nodes` belonging to `region`."""
+        return self._regions[region](nodes)
+
+    @property
+    @abstractmethod
+    def _analytic_predicates(self) -> dict:
+        """Named boundary regions as ``{name: predicate(point) -> bool}``, from this geometry's
+        own shape parameters."""

From 37dc6a6130d0ce92b41ba6426d8651941d86a640 Mon Sep 17 00:00:00 2001
From: Themis Skamagkis <themisskamagkis@gmail.com>
Date: Mon, 12 Oct 2026 09:49:17 +0200
Subject: [PATCH 05/55] [geometry] Add the bar and beam grid geometries

---
 common/geometry/grid.py | 36 ++++++++++++++++++++++++++++++++++++
 1 file changed, 36 insertions(+)
 create mode 100644 common/geometry/grid.py

diff --git a/common/geometry/grid.py b/common/geometry/grid.py
new file mode 100644
index 0000000..d441498
--- /dev/null
+++ b/common/geometry/grid.py
@@ -0,0 +1,36 @@
+"""Rectilinear structured grid geometry in 1 to 3 dimensions: bars and beams."""
+
+import numpy as np
+
+from .geometry import Geometry
+
+
+class Grid(Geometry):
+    """A rectilinear structured grid, axis i spanning 0..parameters[i]. Always starts at the origin."""
+
+    PARAMETER_ORDER = ("length", "width", "height")
+    AXIS_NAMES = (("left", "right"), ("bottom", "top"), ("front", "back"))
+
+    @property
+    def _analytic_predicates(self) -> dict:
+        predicates = {}
+        for axis, (low, high) in enumerate(self.AXIS_NAMES[:self.dim]):
+            extent = self.parameters[axis]
+            predicates[low] = lambda p, axis=axis: np.isclose(p[axis], 0.0)
+            predicates[high] = lambda p, axis=axis, extent=extent: np.isclose(p[axis], extent)
+        return predicates
+
+
+class Bar1D(Grid):
+    """A 1D bar occupying ``[0, length]`` along x."""
+    dim = 1
+
+
+class Beam2D(Grid):
+    """A 2D rectangle occupying ``[0, length] x [0, width]`` in the xy-plane."""
+    dim = 2
+
+
+class Beam3D(Grid):
+    """A 3D box occupying ``[0, length] x [0, width] x [0, height]``."""
+    dim = 3

From 34408da786d0171ea39bd20e8a58ef7f03591764 Mon Sep 17 00:00:00 2001
From: Themis Skamagkis <themisskamagkis@gmail.com>
Date: Mon, 12 Oct 2026 09:49:17 +0200
Subject: [PATCH 06/55] [geometry] Add the tire geometries

---
 common/geometry/tire.py | 32 ++++++++++++++++++++++++++++++++
 1 file changed, 32 insertions(+)
 create mode 100644 common/geometry/tire.py

diff --git a/common/geometry/tire.py b/common/geometry/tire.py
new file mode 100644
index 0000000..2eef1df
--- /dev/null
+++ b/common/geometry/tire.py
@@ -0,0 +1,32 @@
+"""Tire geometry (Chamberland et al., 2010, Fig. 1)."""
+
+import numpy as np
+
+from .geometry import Geometry
+
+
+class Tire(Geometry):
+    """An annulus in the xy-plane, centered on the origin, extruded along z when it is 3D."""
+
+    PARAMETER_ORDER = ("inner_radius", "outer_radius", "extrusion_length")
+
+    @property
+    def _analytic_predicates(self) -> dict:
+        predicates = {
+            "inner": lambda p: np.isclose(np.hypot(p[0], p[1]), self.inner_radius),
+            "outer": lambda p: np.isclose(np.hypot(p[0], p[1]), self.outer_radius),
+        }
+        if self.dim == 3:
+            predicates["front"] = lambda p: np.isclose(p[2], 0.0)
+            predicates["back"] = lambda p: np.isclose(p[2], self.extrusion_length)
+        return predicates
+
+
+class Tire2D(Tire):
+    """An annulus occupying ``inner_radius <= r <= outer_radius`` of the xy-plane."""
+    dim = 2
+
+
+class Tire3D(Tire):
+    """A Tire2D extruded along z from 0 to extrusion_length."""
+    dim = 3

From d01333f331f8d21bd5f46c0ea0b5221f6a4b1288 Mon Sep 17 00:00:00 2001
From: Themis Skamagkis <themisskamagkis@gmail.com>
Date: Mon, 12 Oct 2026 09:49:17 +0200
Subject: [PATCH 07/55] [geometry] Export the geometries

---
 common/geometry/__init__.py | 3 +++
 1 file changed, 3 insertions(+)
 create mode 100644 common/geometry/__init__.py

diff --git a/common/geometry/__init__.py b/common/geometry/__init__.py
new file mode 100644
index 0000000..f82eb80
--- /dev/null
+++ b/common/geometry/__init__.py
@@ -0,0 +1,3 @@
+from .geometry import Geometry
+from .grid import Grid, Bar1D, Beam2D, Beam3D
+from .tire import Tire, Tire2D, Tire3D

From 10b45f0f34f783ffa18e0015897f71e4062e13c7 Mon Sep 17 00:00:00 2001
From: Themis Skamagkis <themisskamagkis@gmail.com>
Date: Mon, 12 Oct 2026 09:49:17 +0200
Subject: [PATCH 08/55] [sofa] Map element kinds to SOFA components

---
 common/sofa/conventions.py | 38 ++++++++++++++++++++++++++++++++++++++
 1 file changed, 38 insertions(+)
 create mode 100644 common/sofa/conventions.py

diff --git a/common/sofa/conventions.py b/common/sofa/conventions.py
new file mode 100644
index 0000000..d14acd8
--- /dev/null
+++ b/common/sofa/conventions.py
@@ -0,0 +1,38 @@
+"""Central map from tool-agnostic vocabulary to SOFA components."""
+
+from dataclasses import dataclass
+
+# Spatial dimension (DOF/embedding space) -> MechanicalObject / force-field vector template.
+VEC_DIM = {1: "Vec1d", 2: "Vec2d", 3: "Vec3d"}
+
+@dataclass(frozen=True)
+class ElementKind:
+    """Everything the scene-building layer needs to know about one type of finite element."""
+
+    dim: int                    # topological dimension of this element
+    container: str              # topology container class
+    data_name: str              # Data field this kind's own container addresses its elements by
+    cpp: str                    # SOFA geometry name, for compound templates ("Vec3d,Hexahedron")
+    facet_kind: str | None      # the kind of this element's own facets, None if it has none
+    grid_mapping: str | None    # topological mapping generating it from the grid, None if grid-native
+    boundary_mappings: tuple    # (element kind, mapping) steps from this element to its boundary
+
+ELEMENTS = {
+    "edge": ElementKind(dim=1, container="EdgeSetTopologyContainer",        data_name="edges",
+                        cpp="Edge",        facet_kind=None,   grid_mapping=None,
+                        boundary_mappings=()),
+    "tri":  ElementKind(dim=2, container="TriangleSetTopologyContainer",    data_name="triangles",
+                        cpp="Triangle",    facet_kind="edge", grid_mapping="Quad2TriangleTopologicalMapping",
+                        boundary_mappings=(("edge", "Triangle2EdgeTopologicalMapping"),)),
+    "quad": ElementKind(dim=2, container="QuadSetTopologyContainer",        data_name="quads",
+                        cpp="Quad",        facet_kind="edge", grid_mapping=None,
+                        # To arrive at a Quad2Edge map -> triangulate quads and then Tri2Edge
+                        boundary_mappings=(("tri", "Quad2TriangleTopologicalMapping"),
+                                           ("edge", "Triangle2EdgeTopologicalMapping"))),
+    "tet":  ElementKind(dim=3, container="TetrahedronSetTopologyContainer", data_name="tetrahedra",
+                        cpp="Tetrahedron", facet_kind="tri",  grid_mapping="Hexa2TetraTopologicalMapping",
+                        boundary_mappings=(("tri", "Tetra2TriangleTopologicalMapping"),)),
+    "hexa": ElementKind(dim=3, container="HexahedronSetTopologyContainer",  data_name="hexahedra",
+                        cpp="Hexahedron",  facet_kind="quad", grid_mapping=None,
+                        boundary_mappings=(("quad", "Hexa2QuadTopologicalMapping"),)),
+}

From b6f010b8dfa6681fd33a0f7aff5499953bdaf70b Mon Sep 17 00:00:00 2001
From: Themis Skamagkis <themisskamagkis@gmail.com>
Date: Mon, 12 Oct 2026 09:49:17 +0200
Subject: [PATCH 09/55] [sofa] Add the SOFA scene base

---
 common/sofa/scene.py | 29 +++++++++++++++++++++++++++++
 1 file changed, 29 insertions(+)
 create mode 100644 common/sofa/scene.py

diff --git a/common/sofa/scene.py b/common/sofa/scene.py
new file mode 100644
index 0000000..87de443
--- /dev/null
+++ b/common/sofa/scene.py
@@ -0,0 +1,29 @@
+"""Generic SOFA scene abstraction. Scenes in VnV build off of it."""
+
+from abc import ABC, abstractmethod
+
+
+class SofaScene(ABC):
+    """Abstract SOFA scene: loads its plugins, then builds the scene graph.
+
+    A subclass appends the plugins it needs to `plugins`, and defines `body` to build the
+    scene graph.
+    """
+
+    plugins = []  # plugins every scene needs; a subclass appends its own on top
+
+    def __init__(self):
+        # The list of plugins is kept as a copy per-instance.
+        self.plugins = list(self.plugins)
+
+    def add_plugins(self, root):
+        """Add every plugin in `plugins` as one RequiredPlugin. Subclass appends to `self.plugins`."""
+        root.addObject('RequiredPlugin', pluginName=self.plugins)
+
+    @abstractmethod
+    def body(self, root):
+        """Build the scene graph."""
+
+    def build(self, root):
+        self.add_plugins(root)
+        return self.body(root)

From 923095132d95c605b6b3ce99b71acb05c6ae261d Mon Sep 17 00:00:00 2001
From: Themis Skamagkis <themisskamagkis@gmail.com>
Date: Mon, 12 Oct 2026 09:49:17 +0200
Subject: [PATCH 10/55] [sofa] Add the region clamp controller

---
 common/sofa/controllers.py | 47 ++++++++++++++++++++++++++++++++++++++
 1 file changed, 47 insertions(+)
 create mode 100644 common/sofa/controllers.py

diff --git a/common/sofa/controllers.py b/common/sofa/controllers.py
new file mode 100644
index 0000000..4bd21e3
--- /dev/null
+++ b/common/sofa/controllers.py
@@ -0,0 +1,47 @@
+"""Post-init SOFA controllers."""
+
+import Sofa.Core
+
+
+class RegionClamp(Sofa.Core.Controller):
+    """Builds one PartialFixedProjectiveConstraint per boundary region, then applies Dirichlet
+    boundary conditions once the simulation has initialized.
+
+    For each region, the nodes belonging to it are found and handed to its constraint as the indices
+    it fixes. If a prescribed displacement is given, those nodes are also moved from rest to rest +
+    u in the constraint's fixed directions before the constraint clamps them there. rest_position
+    itself is never changed.
+    """
+
+    def __init__(self, geometry, dofs, node, prescribe_displacement_on, vec_type,
+                 displacement=None, *args, **kwargs):
+        super().__init__(*args, **kwargs)
+        self.geometry = geometry
+        self.dofs = dofs
+        self.displacement = displacement
+
+        # Fix dofs partially as determined by MMS `prescribe_displacement_on`
+        self.region_constraints = []
+        for region, fixed_directions in prescribe_displacement_on.items():
+            # Pad the array with zeroes for embedded dimensions: e.g. [1, 0] -> [1, 0, 0]
+            padded = tuple(fixed_directions) + (0,) * (geometry.spatial_dimensions - len(fixed_directions))
+            constraint = node.addObject('PartialFixedProjectiveConstraint', name='clamp' + region,
+                                  template=vec_type, fixedDirections=list(padded))
+            self.region_constraints.append((region, padded, constraint))
+
+    def onSimulationInitDoneEvent(self, event):
+        rest_positions = self.dofs.rest_position.array()
+
+        # Apply non-zero Dirichlet BC: Move each fixed direction to rest + u; rest_position is untouched.
+        with self.dofs.position.writeableArray() as pos:
+            for region, fixed_directions, constraint in self.region_constraints:
+                indices = self.geometry.region_indices(region, rest_positions)
+
+                if self.displacement is not None and indices:
+                    displacement = self.displacement(rest_positions[indices])
+
+                    for axis, fixed in enumerate(fixed_directions):
+                        if fixed:
+                            pos[indices, axis] = rest_positions[indices, axis] + displacement[:, axis]
+
+                constraint.indices.value = indices

From 2df708a6e60a494ac2c7327320938690d4ba0a88 Mon Sep 17 00:00:00 2001
From: Themis Skamagkis <themisskamagkis@gmail.com>
Date: Mon, 12 Oct 2026 09:49:17 +0200
Subject: [PATCH 11/55] [sofa] Add the prefab base: force field and solvers

---
 common/sofa/prefabs/base.py | 62 +++++++++++++++++++++++++++++++++++++
 1 file changed, 62 insertions(+)
 create mode 100644 common/sofa/prefabs/base.py

diff --git a/common/sofa/prefabs/base.py b/common/sofa/prefabs/base.py
new file mode 100644
index 0000000..1603d4a
--- /dev/null
+++ b/common/sofa/prefabs/base.py
@@ -0,0 +1,62 @@
+"""Prefabs a study configures: they take the case blocks they need and build themselves."""
+
+import Sofa
+import Sofa.Core
+
+
+def params(config, *exclude):
+    """Converts a set of parameters in a configuration dictionary to SOFA Data/parameterization.
+    Every key is passed on except 'type' and excluded ones marked by the user through *exclude."""
+    skip = ('type',) + exclude
+    return {k: v for k, v in config.items() if k not in skip}
+
+
+class ScenePrefab(Sofa.Prefab):
+    """Base SOFA prefab: holds `configs` and offers adders a subclass calls to build itself from it."""
+
+    def __init__(self, *args, configs=None, geometry=None, **kwargs):
+        # object.__setattr__: the C++ object Sofa's own __setattr__ needs doesn't exist yet, and
+        # Sofa.Prefab.__init__ below calls init(), which reads configs/geometry
+        object.__setattr__(self, 'configs', configs or {})
+        object.__setattr__(self, 'geometry', geometry)
+        Sofa.Prefab.__init__(self, *args, **kwargs)
+
+    def mechanical_node(self):
+        raise NotImplementedError
+
+    def add_force_field(self, node, config, material, template):
+        """The *FEMForceField component under test, its Data merged from `config` and `material`"""
+        node.addObject(config['type'], name='fem', template=template, topology='@topology',
+                       **params(material), **params(config))
+
+    def add_solvers(self, node, config):
+        """Newton (if present) + linear solver (+ preconditioner) + integration scheme."""
+
+        # Newton Solver
+        newton_config = config.get('newton')
+        if newton_config is not None:
+            node.addObject(newton_config['type'], name='newton', **params(newton_config))
+
+        # Linear Solver & Preconditioner (optional)
+        linear_config = config['linearSolver']
+        precond_config = linear_config.get('preconditioner')
+        if precond_config is None:
+            node.addObject(linear_config['type'], name='linearSolver', **params(linear_config))
+        else:
+            precond_system_config = precond_config['system']
+            node.addObject(precond_system_config['type'], name='precondSystem', **params(precond_system_config))
+            node.addObject(precond_config['type'], name='precond', linearSystem='@precondSystem',
+                           **params(precond_config, 'system'))
+
+            linear_system_config = linear_config['system']
+            node.addObject(linear_system_config['type'], name='solverSystem',
+                           preconditionerSystem='@precondSystem', **params(linear_system_config))
+            node.addObject(linear_config['type'], name='linearSolver', linearSystem='@solverSystem',
+                           preconditioner='@precond', **params(linear_config, 'system', 'preconditioner'))
+
+        # ODE Integration Scheme
+        integration_scheme = params(config['integration'])
+        integration_scheme['linearSolver'] = '@linearSolver'
+        if newton_config is not None:
+            integration_scheme['newtonSolver'] = '@newton'
+        node.addObject(config['integration']['type'], name='ode', **integration_scheme)

From 4b120bb1fb022e4265f5a04b8515abf8eb887b08 Mon Sep 17 00:00:00 2001
From: Themis Skamagkis <themisskamagkis@gmail.com>
Date: Mon, 12 Oct 2026 09:49:18 +0200
Subject: [PATCH 12/55] [sofa] Add the elastic grid prefabs

---
 common/sofa/prefabs/grid.py | 89 +++++++++++++++++++++++++++++++++++++
 1 file changed, 89 insertions(+)
 create mode 100644 common/sofa/prefabs/grid.py

diff --git a/common/sofa/prefabs/grid.py b/common/sofa/prefabs/grid.py
new file mode 100644
index 0000000..8af2906
--- /dev/null
+++ b/common/sofa/prefabs/grid.py
@@ -0,0 +1,89 @@
+"""SOFA prefabs: an elastic rectilinear structured grid in 1 to 3 dimensions, as a bar or a beam."""
+
+from .base import ScenePrefab, params
+from ..conventions import VEC_DIM, ELEMENTS
+from ...geometry import Grid
+
+# Element kinds this prefab meshes, by the topological dimension of the geometry.
+ELEMENT_KINDS = {1: ('edge',), 2: ('quad', 'tri'), 3: ('hexa', 'tet')}
+
+
+def validate_parameters(name, config):
+    """Required ElasticGrid parameters, and an element kind its geometry can be meshed with."""
+    required = ['geometry', 'grid_resolution', 'element']
+    missing = [p for p in required if p not in config]
+    if missing:
+        raise ValueError(f"{name}: missing required parameters {missing}")
+    geometry = config['geometry']
+    if not isinstance(geometry, Grid):
+        raise ValueError(f"{name}: geometry must be a Grid, "
+                         f"got {type(geometry).__name__}")
+    kinds = ELEMENT_KINDS[geometry.dim]
+    if config['element'] not in kinds:
+        raise ValueError(f"{name}: element on a {geometry.dim}D {type(geometry).__name__} must be "
+                         f"one of {list(kinds)}, got {config['element']!r}")
+
+
+class ElasticGrid(ScenePrefab):
+    """Prefab for an elastic rectilinear structured grid model in SOFA."""
+
+    # The mechanical node's name in the scene graph; a named grid calls it what it is.
+    BODY = 'Body'
+
+    prefabParameters = [
+        {'name': 'grid_resolution', 'type': 'Vec3d',  'help': 'grid points per axis [nx, ny, nz]'},
+        {'name': 'element',    'type': 'string', 'help': 'element kind (edge/tri/quad/tet/hexa)'},
+    ]
+
+    def __init__(self, *args, **kwargs):
+        validate_parameters(type(self).__name__, kwargs)
+        super().__init__(*args, **kwargs)
+
+    def init(self):
+        dim = self.geometry.dim
+        extents = list(self.geometry.parameters) + [0.0] * (3 - dim)
+        VecType = VEC_DIM[self.geometry.spatial_dimensions]
+        element_kind = ELEMENTS[self.element.value]
+        grid_resolution = list(self.grid_resolution.value)[:dim] + [1] * (3 - dim)
+
+        # Grid Topology Node
+        with self.addChild('Grid') as grid_node:
+            grid_node.addObject('RegularGridTopology', name='grid',
+                                nx=int(grid_resolution[0]), ny=int(grid_resolution[1]),
+                                nz=int(grid_resolution[2]), min=[0.0, 0.0, 0.0], max=extents)
+
+        # Node containing the body's components
+        with self.addChild(self.BODY) as body:
+            self.body = body
+
+            # Topology
+            if element_kind.grid_mapping is None:
+                body.addObject(element_kind.container, name='topology', position='@../Grid/grid.position',
+                               **{element_kind.data_name: f'@../Grid/grid.{element_kind.data_name}'})
+            else:
+                body.addObject(element_kind.container, name='topology', position='@../Grid/grid.position')
+                body.addObject(element_kind.grid_mapping, input='@../Grid/grid', output='@topology',
+                               **params(self.configs.get('gridMapping', {})))
+                body.addObject(element_kind.container.replace('Container', 'Modifier'))
+            # DOFs
+            body.addObject('MechanicalObject', name='dofs', template=VecType)
+            # The force field under test.
+            self.add_force_field(body, self.configs['forceField'], self.configs['material'],
+                                 template=f"{VecType},{element_kind.cpp}")
+            # ODE & Linear Solvers
+            self.add_solvers(body, self.configs['solvers'])
+
+    def mechanical_node(self):
+        return self.body
+
+
+class ElasticBar(ElasticGrid):
+    """Prefab for an elastic bar model in SOFA."""
+
+    BODY = 'Bar'
+
+
+class ElasticBeam(ElasticGrid):
+    """Prefab for an elastic beam model in SOFA."""
+
+    BODY = 'Beam'

From 16092380e64b8dad22f100ff2e66f7341cd9fbb1 Mon Sep 17 00:00:00 2001
From: Themis Skamagkis <themisskamagkis@gmail.com>
Date: Mon, 12 Oct 2026 09:49:18 +0200
Subject: [PATCH 13/55] [sofa] Add the elastic tire prefab

---
 common/sofa/prefabs/tire.py | 186 ++++++++++++++++++++++++++++++++++++
 1 file changed, 186 insertions(+)
 create mode 100644 common/sofa/prefabs/tire.py

diff --git a/common/sofa/prefabs/tire.py b/common/sofa/prefabs/tire.py
new file mode 100644
index 0000000..0afb488
--- /dev/null
+++ b/common/sofa/prefabs/tire.py
@@ -0,0 +1,186 @@
+"""SOFA prefab: an elastic tire -- an annulus, extruded along z when it is 3D."""
+
+import itertools
+
+import numpy as np
+
+from .base import ScenePrefab
+from ..conventions import VEC_DIM, ELEMENTS
+from ...geometry import Tire
+
+# Element kinds this prefab meshes, by the topological dimension of the geometry.
+ELEMENT_KINDS = {2: ('quad', 'tri'), 3: ('hexa', 'tet')}
+
+# The two triangles one quadrilateral splits into, on one diagonal and on the other.
+_TRIANGLE_PATTERNS = (((1, 2, 0), (3, 0, 2)),
+                      ((0, 1, 3), (2, 3, 1)))
+
+# Swapping the two nodes of every hexahedron edge along one axis, as the permutation it applies to
+# the eight corner slots. sofa::geometry::Hexahedron::{x,y,z}Edges.
+_EDGE_SWAPS = ((1, 0, 3, 2, 5, 4, 7, 6),
+               (3, 2, 1, 0, 7, 6, 5, 4),
+               (4, 5, 6, 7, 0, 1, 2, 3))
+
+# The six tetrahedra one hexahedron splits into, before and after an odd number of those swaps.
+_TETRA_PATTERNS = {
+    False: ((0, 5, 1, 6), (0, 1, 3, 6), (1, 3, 6, 2), (6, 3, 0, 7), (6, 7, 0, 5), (7, 5, 4, 0)),
+    True:  ((0, 5, 6, 1), (0, 1, 6, 3), (1, 3, 2, 6), (6, 3, 7, 0), (6, 7, 5, 0), (7, 5, 0, 4)),
+}
+
+
+def validate_parameters(config):
+    """Required ElasticTire parameters, and an element kind its geometry can be meshed with."""
+    required = ['geometry', 'grid_resolution', 'element']
+    missing = [p for p in required if p not in config]
+    if missing:
+        raise ValueError(f"ElasticTire: missing required parameters {missing}")
+    geometry = config['geometry']
+    if not isinstance(geometry, Tire):
+        raise ValueError(f"ElasticTire: geometry must be a Tire2D or a Tire3D, "
+                         f"got {type(geometry).__name__}")
+    kinds = ELEMENT_KINDS[geometry.dim]
+    if config['element'] not in kinds:
+        raise ValueError(f"ElasticTire: element on a {type(geometry).__name__} must be one of "
+                         f"{list(kinds)}, got {config['element']!r}")
+
+
+def annulus(geometry, cells):
+    """The annulus in the xy-plane: its nodes, and the quadrilaterals over them.
+
+    `cells` counts cells radially and circumferentially. The ring closes on itself, so the
+    circumferential direction holds exactly that many nodes, not one more.
+    """
+    radial, circumferential = cells
+
+    radii = np.linspace(geometry.inner_radius, geometry.outer_radius, radial + 1)
+    angles = np.linspace(0.0, 2.0 * np.pi, circumferential, endpoint=False)
+
+    r, theta = np.meshgrid(radii, angles, indexing='ij')
+    nodes = np.stack([r * np.cos(theta), r * np.sin(theta), np.zeros_like(r)],
+                     axis=-1).reshape(-1, 3)
+
+    node = np.arange(nodes.shape[0]).reshape(r.shape)
+    i, j = np.mgrid[0:radial, 0:circumferential]
+    ring = (j + 1) % circumferential   # the last cell of the ring closes back onto j = 0
+
+    # (radial, circumferential) turns counter-clockwise, the winding GridTopology gives a cell.
+    corners = (node[i, j], node[i + 1, j], node[i + 1, ring], node[i, ring])
+    return nodes, np.stack(corners, axis=-1).reshape(-1, 4)
+
+
+def extrude(nodes, quads, length, axial):
+    """The annulus swept along z into `axial` layers: the swept nodes, and the hexahedra over them.
+
+    Each layer repeats the annulus, so a hexahedron is its quadrilateral followed by the same one
+    a layer up -- again the winding GridTopology gives a cell.
+    """
+    layer = nodes.shape[0]
+    swept = np.tile(nodes, (axial + 1, 1))
+    swept[:, 2] = np.repeat(np.linspace(0.0, length, axial + 1), layer)
+
+    base = np.arange(axial)[:, None, None] * layer
+    hexahedra = np.concatenate([quads[None] + base, quads[None] + base + layer], axis=-1)
+    return swept, hexahedra.reshape(-1, 8)
+
+
+def triangulate(quads, cells):
+    """Two triangles per quadrilateral, split the way Quad2TriangleTopologicalMapping splits one.
+
+    Neighbours share a whole edge, so any choice of diagonal conforms; alternating it cell by cell
+    only keeps the mesh from leaning the same way everywhere.
+    """
+    i, j = np.mgrid[0:cells[0], 0:cells[1]]
+    diagonal = ((i ^ j) & 1).reshape(-1)
+
+    triangles = np.empty((quads.shape[0], 2, 3), dtype=np.int64)
+    for choice, pattern in enumerate(_TRIANGLE_PATTERNS):
+        selected = diagonal == choice
+        triangles[selected] = quads[selected][:, pattern]
+
+    return triangles.reshape(-1, 3)
+
+
+def tetrahedralize(hexahedra, cells):
+    """Six tetrahedra per hexahedron, split the way Hexa2TetraTopologicalMapping splits one.
+
+    The node labelling alternates with the parity of the cell's index, so neighbours agree on the
+    diagonal of every face they share. Across the seam that pairs the last circumferential cell
+    with the first, which only agree when there is an even number of them.
+    """
+    radial, circumferential, axial = cells
+    if circumferential % 2:
+        raise ValueError(f"ElasticTire: an odd number of circumferential cells "
+                         f"({circumferential}) leaves the tetrahedra of the closing cell "
+                         f"disagreeing with the first on their shared diagonal")
+
+    # Layer by layer, as `extrude` stacks them.
+    k, i, j = np.mgrid[0:axial, 0:radial, 0:circumferential]
+    parities = np.stack([i % 2 == 0, j % 2 == 1, k % 2 == 1], axis=-1).reshape(-1, 3)
+
+    tetrahedra = np.empty((hexahedra.shape[0], 6, 4), dtype=np.int64)
+    for signature in itertools.product((False, True), repeat=3):
+        slots, swapped = list(range(8)), False
+        for axis, swap in enumerate(signature):
+            if swap:
+                slots = [slots[p] for p in _EDGE_SWAPS[axis]]
+                swapped = not swapped
+        pattern = [[slots[corner] for corner in tetra] for tetra in _TETRA_PATTERNS[swapped]]
+
+        selected = np.all(parities == signature, axis=-1)
+        tetrahedra[selected] = hexahedra[selected][:, pattern]
+
+    return tetrahedra.reshape(-1, 4)
+
+
+def mesh(geometry, cells, element):
+    """The nodes of `geometry` and its elements of the given kind."""
+    nodes, elements = annulus(geometry, cells[:2])
+    if geometry.dim == 3:
+        nodes, elements = extrude(nodes, elements, geometry.extrusion_length, cells[2])
+
+    if element == 'tri':
+        elements = triangulate(elements, cells)
+    elif element == 'tet':
+        elements = tetrahedralize(elements, cells)
+    return nodes, elements
+
+
+class ElasticTire(ScenePrefab):
+    """Prefab for an elastic tire model in SOFA."""
+
+    prefabParameters = [
+        {'name': 'grid_resolution', 'type': 'Vec3d',
+         'help': 'nodes per axis [radial, circumferential] and [, axial] when extruded'},
+        {'name': 'element',    'type': 'string', 'help': 'element kind (quad/tri/hexa/tet)'},
+    ]
+
+    def __init__(self, *args, **kwargs):
+        validate_parameters(kwargs)
+        super().__init__(*args, **kwargs)
+
+    def init(self):
+        VecType = VEC_DIM[self.geometry.spatial_dimensions]
+        element_kind = ELEMENTS[self.element.value]
+        # A study states nodes per axis, the mesh is laid out in cells.
+        cells = [int(count) - 1 for count in
+                 list(self.grid_resolution.value)[:self.geometry.dim]]
+
+        nodes, elements = mesh(self.geometry, cells, self.element.value)
+
+        # Node containing Tire components
+        with self.addChild('Tire') as tire:
+            self.tire = tire
+
+            # Topology: meshed directly, so no grid and no topological mapping to reach it.
+            tire.addObject(element_kind.container, name='topology', position=nodes.tolist(),
+                           **{element_kind.data_name: elements.tolist()})
+            # DOFs
+            tire.addObject('MechanicalObject', name='dofs', template=VecType)
+            # The force field under test.
+            self.add_force_field(tire, self.configs['forceField'], self.configs['material'],
+                                 template=f"{VecType},{element_kind.cpp}")
+            # ODE & Linear Solvers
+            self.add_solvers(tire, self.configs['solvers'])
+
+    def mechanical_node(self):
+        return self.tire

From 2f9dc3a0369366852c4f5e9ce2952a898056c423 Mon Sep 17 00:00:00 2001
From: Themis Skamagkis <themisskamagkis@gmail.com>
Date: Mon, 12 Oct 2026 09:49:18 +0200
Subject: [PATCH 14/55] [sofa] Map geometries to their prefab

---
 common/sofa/prefabs/__init__.py | 8 ++++++++
 1 file changed, 8 insertions(+)
 create mode 100644 common/sofa/prefabs/__init__.py

diff --git a/common/sofa/prefabs/__init__.py b/common/sofa/prefabs/__init__.py
new file mode 100644
index 0000000..8ea7478
--- /dev/null
+++ b/common/sofa/prefabs/__init__.py
@@ -0,0 +1,8 @@
+from .base import ScenePrefab
+from .grid import ElasticGrid, ElasticBar, ElasticBeam
+from .tire import ElasticTire
+from ...geometry import Grid, Bar1D, Beam2D, Beam3D, Tire, Tire2D, Tire3D
+
+PREFAB_BY_GEOMETRY = {Grid: ElasticGrid, Bar1D: ElasticBar,
+                      Beam2D: ElasticBeam, Beam3D: ElasticBeam,
+                      Tire: ElasticTire, Tire2D: ElasticTire, Tire3D: ElasticTire}

From 1838d6a5d95d30b37fe9f16d565e531eaf32843a Mon Sep 17 00:00:00 2001
From: Themis Skamagkis <themisskamagkis@gmail.com>
Date: Mon, 12 Oct 2026 09:49:18 +0200
Subject: [PATCH 15/55] [sofa] Export the SOFA scene base

---
 common/sofa/__init__.py | 1 +
 1 file changed, 1 insertion(+)
 create mode 100644 common/sofa/__init__.py

diff --git a/common/sofa/__init__.py b/common/sofa/__init__.py
new file mode 100644
index 0000000..3175991
--- /dev/null
+++ b/common/sofa/__init__.py
@@ -0,0 +1 @@
+from .scene import SofaScene

From 19e75106c61205e720f5879ec5cd8ed3505de796 Mon Sep 17 00:00:00 2001
From: Themis Skamagkis <themisskamagkis@gmail.com>
Date: Mon, 12 Oct 2026 09:49:18 +0200
Subject: [PATCH 16/55] [verification] Add the verification package

---
 verification/__init__.py | 0
 1 file changed, 0 insertions(+), 0 deletions(-)
 create mode 100644 verification/__init__.py

diff --git a/verification/__init__.py b/verification/__init__.py
new file mode 100644
index 0000000..e69de29

From d0321d3d852d8438f753996da46779955e68974a Mon Sep 17 00:00:00 2001
From: Themis Skamagkis <themisskamagkis@gmail.com>
Date: Mon, 12 Oct 2026 09:49:18 +0200
Subject: [PATCH 17/55] [verification] Add the linear elastic material law

---
 verification/manufactured/materials.py | 131 +++++++++++++++++++++++++
 1 file changed, 131 insertions(+)
 create mode 100644 verification/manufactured/materials.py

diff --git a/verification/manufactured/materials.py b/verification/manufactured/materials.py
new file mode 100644
index 0000000..d804836
--- /dev/null
+++ b/verification/manufactured/materials.py
@@ -0,0 +1,131 @@
+"""Constitutive laws: a material declares psi(gradient); its stress and tangent are differentiated from that."""
+
+from abc import ABC, abstractmethod
+
+import numpy as np
+import sympy as sp
+
+
+def _toLameParameters1D(youngModulus, poissonRatio):
+    """(mu, lambda) at d = 1, where lambda + 2 mu = E collapses Hooke to sigma = E eps."""
+    return 0.5 * youngModulus, 0.0
+
+
+def _toLameParameters2D(youngModulus, poissonRatio):
+    """(mu, lambda) at d = 2, in plane stress: the convention SOFA's force fields assume."""
+    return (youngModulus / (2 * (1 + poissonRatio)),
+            youngModulus * poissonRatio / (1 - poissonRatio ** 2))
+
+
+def _toLameParameters3D(youngModulus, poissonRatio):
+    """(mu, lambda) at d = 3."""
+    return (youngModulus / (2 * (1 + poissonRatio)),
+            youngModulus * poissonRatio / ((1 + poissonRatio) * (1 - 2 * poissonRatio)))
+
+
+# dim -> Young/Poisson -> (mu, lambda)
+_toLame = {1: _toLameParameters1D,
+         2: _toLameParameters2D,
+         3: _toLameParameters3D}
+
+
+def _tensor(dimensions, label):
+    """A d x d matrix of independent symbols, named "<label>_i_j", one per entry.
+
+    Stand-in for a real tensor (a gradient, a direction, ...) so an expression can be differentiated
+    symbolically before any actual numbers exist.
+    """
+    return sp.Matrix(dimensions, dimensions, lambda i, j: sp.Symbol(f"{label}_{i}_{j}"))
+
+
+class Material(ABC):
+    """A material as one strain energy density psi(gradient); stress and tangent come from it.
+
+    psi is the elastic energy stored per unit volume.
+    It is defined for a given displacement gradient (a d x d matrix, d = spatial dimension).
+    Everything else is a derivative of psi:
+      stress(gradient)                        = d(psi)/d(gradient)                     -- a d x d tensor
+      directional_tangent(gradient, direction) = d(stress)/d(gradient), applied to `direction` -- also d x d
+
+    Preferably, for any given material, only psi is declared. stress and tangent are derived by 
+    symbolic differentiation.
+    """
+
+    @classmethod
+    @abstractmethod
+    def get_material_parameters(cls, parameters, spatial_dimensions):
+        """Return the constant parameters defining the material law from a dictionary."""
+
+    @abstractmethod
+    def energy_density(self, gradient):
+        """psi(gradient): the strain energy density, an expression in `gradient`'s symbolic entries."""
+
+    def stress(self, gradient):
+        """d(psi)/d(gradient), evaluated at the given `gradient`."""
+        gradient_symbols = _tensor(gradient.rows, "G")
+        return self._stress(gradient_symbols).subs(dict(zip(gradient_symbols, gradient)))
+
+    def directional_tangent(self, gradient, direction):
+        """d(stress)/d(gradient) applied to `direction`, evaluated at `gradient`."""
+        gradient_symbols = _tensor(gradient.rows, "G")
+        symbolic_stress = self._stress(gradient_symbols)
+        rows, columns = gradient_symbols.rows, gradient_symbols.cols
+        result = sp.Matrix(rows, columns, lambda i, j: sum(
+            sp.diff(symbolic_stress[i, j], gradient_symbols[k, l]) * direction[k, l]
+            for k in range(rows) for l in range(columns)))
+        return result.subs(dict(zip(gradient_symbols, gradient)))
+
+    def _stress(self, gradient_symbols):
+        """The stress formula in terms of `gradient_symbols`"""
+        psi = self.energy_density(gradient_symbols)
+        return sp.Matrix(gradient_symbols.rows, gradient_symbols.cols,
+                         lambda i, j: sp.diff(psi, gradient_symbols[i, j]))
+
+
+def compile_laws(material, dimensions):
+    """Compile and hand back a material's psi and tangent into plain numpy functions over batches of gradients."""
+    gradient, direction = _tensor(dimensions, "G"), _tensor(dimensions, "D")
+    arguments = list(gradient) + list(direction)
+    energy_law = sp.lambdify(list(gradient), material.energy_density(gradient), "numpy")
+    directional_tangent_law = [sp.lambdify(arguments, entry, "numpy")
+                               for entry in material.directional_tangent(gradient, direction)]
+
+    def columns_of(values):
+        values = np.asarray(values, dtype=float)
+        return values, [values[..., i, j]
+                        for i in range(dimensions) for j in range(dimensions)]
+
+    def compiled_energy_density(gradients):
+        values, columns = columns_of(gradients)
+        return np.broadcast_to(np.asarray(energy_law(*columns), dtype=float), values.shape[:-2])
+
+    def compiled_directional_tangent(at, applied_to):
+        base, base_columns = columns_of(at)
+        applied, applied_columns = columns_of(applied_to)
+        batch = np.broadcast_shapes(base.shape[:-2], applied.shape[:-2])
+        result = np.empty(batch + (dimensions, dimensions), dtype=float)
+        for (i, j), component in zip(np.ndindex(dimensions, dimensions), directional_tangent_law):
+            result[..., i, j] = np.broadcast_to(
+                np.asarray(component(*base_columns, *applied_columns), dtype=float), batch)
+        return result
+
+    return compiled_energy_density, compiled_directional_tangent
+
+
+# --- The materials themselves, one class per type of material. ---
+
+
+class LinearElastic(Material):
+    """Hooke: psi = 1/2 lambda tr(eps)^2 + mu eps:eps on the symmetric gradient eps."""
+
+    @classmethod
+    def get_material_parameters(cls, parameters, spatial_dimensions):
+        return _toLame[spatial_dimensions](parameters["youngModulus"], parameters["poissonRatio"])
+
+    def __init__(self, mu, lam):
+        self.mu = mu
+        self.lam = lam
+
+    def energy_density(self, gradient):
+        strain = (gradient + gradient.T) / 2
+        return self.lam * strain.trace() ** 2 / 2 + self.mu * sum(e ** 2 for e in strain)

From cb7e5b9795895a5c66f856ddae128d6f5601b7ff Mon Sep 17 00:00:00 2001
From: Themis Skamagkis <themisskamagkis@gmail.com>
Date: Mon, 12 Oct 2026 09:49:18 +0200
Subject: [PATCH 18/55] [verification] Add the manufactured solution base

---
 verification/manufactured/base.py | 62 +++++++++++++++++++++++++++++++
 1 file changed, 62 insertions(+)
 create mode 100644 verification/manufactured/base.py

diff --git a/verification/manufactured/base.py b/verification/manufactured/base.py
new file mode 100644
index 0000000..4ec0f83
--- /dev/null
+++ b/verification/manufactured/base.py
@@ -0,0 +1,62 @@
+"""Manufactured solutions: exact displacements declared symbolically, and the regions where
+their BCs apply.
+"""
+
+from abc import ABC, abstractmethod
+
+import sympy as sp
+
+
+class ManufacturedSolution(ABC):
+    """An exact displacement, declared symbolically, and the regions where its BCs apply.
+
+    Stated without reference to a material; `ManufacturedProblem` pairs the two.
+    """
+
+    geometry: type                    # the geometry family this solution is stated on
+    dim: int                          # the PDE's own dimension, not the space it is solved in
+
+    def __init__(self, config, geometry, traction_on=None):
+        if not isinstance(geometry, self.geometry):
+            raise ValueError(f"{type(self).__name__} is stated on a {self.geometry.__name__}, "
+                             f"got a {type(geometry).__name__}")
+        if geometry.dim != self.dim:
+            raise ValueError(f"{type(self).__name__} is {self.dim}D, "
+                             f"got a {geometry.dim}D {type(geometry).__name__}")
+
+        self.config = config                        # the case's "solution" block; a subclass may read more from it
+        self.amplitude = config["amplitude"]        # too large relative to the geometry inverts elements
+        self.parameters = geometry.named_parameters
+        self.regions = list(geometry.region_names)
+        self.traction_on = traction_on or {}        # {region: directions} the case hands to traction
+
+    @property
+    def prescribe_displacement_on(self):
+        """{region: fixed_directions}: every direction of every region, less the ones the case
+        handed to the traction. A region left with nothing fixed is dropped.
+        """
+        prescribed = {}
+        for region in self.regions:
+            free = self.traction_on.get(region, [0] * self.dim)
+            fixed_directions = [1 - direction for direction in free]
+            if any(fixed_directions):
+                prescribed[region] = fixed_directions
+        return prescribed
+
+    @abstractmethod
+    def displacement(self, coordinates):
+        """Exact displacement as sympy expressions, one per component of this solution's own dimension."""
+
+
+class TrigonometricSolution(ManufacturedSolution):
+    """A solution built from sines and cosines fitted to the geometry's extent."""
+
+    def __init__(self, config, geometry, traction_on=None):
+        super().__init__(config, geometry, traction_on)
+        self.periods = config.get("periods", 1)     # more periods need a finer mesh to resolve
+
+    def wavenumber(self, parameter):
+        """The wavenumber fitting `self.periods` full periods across one of the geometry's
+        characteristic parameters, e.g. 2 pi / length for a single period.
+        """
+        return 2 * sp.pi * sp.Rational(str(self.periods)) / sp.Rational(str(self.parameters[parameter]))

From 20db28087326504bf84dbc2b19ee9ce5ea271826 Mon Sep 17 00:00:00 2001
From: Themis Skamagkis <themisskamagkis@gmail.com>
Date: Mon, 12 Oct 2026 09:49:18 +0200
Subject: [PATCH 19/55] [verification] Add the grid manufactured solutions

---
 verification/manufactured/grid.py | 97 +++++++++++++++++++++++++++++++
 1 file changed, 97 insertions(+)
 create mode 100644 verification/manufactured/grid.py

diff --git a/verification/manufactured/grid.py b/verification/manufactured/grid.py
new file mode 100644
index 0000000..b43d346
--- /dev/null
+++ b/verification/manufactured/grid.py
@@ -0,0 +1,97 @@
+"""Manufactured solutions stated on the grid geometry."""
+
+import sympy as sp
+
+from common.geometry import Grid
+
+from .base import ManufacturedSolution, TrigonometricSolution
+
+
+class Quadratic1D(ManufacturedSolution):
+    """u(x) = [A x^2]: body force constant, so the nodal source is the source; clamped at 'left'."""
+
+    geometry = Grid
+    dim = 1
+
+    def displacement(self, coordinates):
+        return [self.amplitude * coordinates[0] ** 2]
+
+
+class Quadratic2D(ManufacturedSolution):
+    """u(x, y) = [A (x^2 + x y), A (y^2 + x y)]
+
+    Body force constant and traction linear, so both are nodally exact.
+    """
+
+    geometry = Grid
+    dim = 2
+
+    def displacement(self, coordinates):
+        x, y = coordinates[0], coordinates[1]
+        return [self.amplitude * (x**2 + x * y), self.amplitude * (y**2 + x * y)]
+
+
+class Quadratic3D(ManufacturedSolution):
+    """u(x, y, z) = [A (x^2 + x y + x z), A (y^2 + y z + x y), A (z^2 + z x + y z)]
+
+    Body force constant and traction linear, so both are nodally exact.
+    """
+
+    geometry = Grid
+    dim = 3
+
+    def displacement(self, coordinates):
+        x, y, z = coordinates
+        return [self.amplitude * (x**2 + x * y + x * z),
+                self.amplitude * (y**2 + y * z + x * y),
+                self.amplitude * (z**2 + z * x + y * z)]
+
+
+class Trigonometric1D(TrigonometricSolution):
+    """u(x) = [A sin(k x)], k = 2 pi / L: prescribed (clamped) at 'left', traction at 'right'."""
+
+    geometry = Grid
+    dim = 1
+
+    def displacement(self, coordinates):
+        return [self.amplitude * sp.sin(self.wavenumber("length") * coordinates[0])]
+
+
+class Trigonometric2D(TrigonometricSolution):
+    """u(x, y) = [A sin(kx x) cos(ky y),
+                 A cos(kx x) sin(ky y)]
+    kx = 2 pi / L, ky = 2 pi / W
+
+    u_x fixed on the left face, u_y on the bottom face, traction elsewhere.
+    """
+
+    geometry = Grid
+    dim = 2
+
+    def displacement(self, coordinates):
+        kx, ky = self.wavenumber("length"), self.wavenumber("width")
+        x, y = coordinates[0], coordinates[1]
+        return [self.amplitude * sp.sin(kx * x) * sp.cos(ky * y),
+                self.amplitude * sp.cos(kx * x) * sp.sin(ky * y)]
+
+
+class Trigonometric3D(TrigonometricSolution):
+    """u(x, y, z) = [A sin(kx x) cos(ky y) cos(kz z),
+                    A cos(kx x) sin(ky y) cos(kz z),
+                    A cos(kx x) cos(ky y) sin(kz z)]
+    kx = 2 pi / L, ky = 2 pi / W, kz = 2 pi / H
+
+    Each component fixed on the low face of its own axis, traction elsewhere.
+    """
+
+    geometry = Grid
+    dim = 3
+
+    def displacement(self, coordinates):
+        kx, ky, kz = (self.wavenumber("length"), self.wavenumber("width"),
+                      self.wavenumber("height"))
+        x, y, z = coordinates
+        return [self.amplitude * sp.sin(kx * x) * sp.cos(ky * y) * sp.cos(kz * z),
+                self.amplitude * sp.cos(kx * x) * sp.sin(ky * y) * sp.cos(kz * z),
+                self.amplitude * sp.cos(kx * x) * sp.cos(ky * y) * sp.sin(kz * z)]
+

From da047fd0516a4fe9ecc1ac9f23a7ab9a17a47dba Mon Sep 17 00:00:00 2001
From: Themis Skamagkis <themisskamagkis@gmail.com>
Date: Mon, 12 Oct 2026 09:49:18 +0200
Subject: [PATCH 20/55] [verification] Add the tire manufactured solutions

---
 verification/manufactured/tire.py | 45 +++++++++++++++++++++++++++++++
 1 file changed, 45 insertions(+)
 create mode 100644 verification/manufactured/tire.py

diff --git a/verification/manufactured/tire.py b/verification/manufactured/tire.py
new file mode 100644
index 0000000..a0a99fc
--- /dev/null
+++ b/verification/manufactured/tire.py
@@ -0,0 +1,45 @@
+"""Manufactured solutions stated on the tire geometry (Chamberland et al., 2010, Sec. 5.2)."""
+
+import sympy as sp
+
+from common.geometry import Tire
+
+from .base import ManufacturedSolution
+
+
+class Quartic2D(ManufacturedSolution):
+    """u(x, y) = [A (x^4 + 2 x y / 5),
+                  A (y^4 / 10 - 2 x y / 5)]
+
+    u prescribed on the inner edge, traction on the outer.
+    """
+
+    geometry = Tire
+    dim = 2
+
+    def displacement(self, coordinates):
+        x, y = coordinates[0], coordinates[1]
+        return [self.amplitude * (x**4 + sp.Rational(2, 5) * x * y),
+                self.amplitude * (sp.Rational(1, 10) * y**4 - sp.Rational(2, 5) * x * y)]
+
+
+class Quartic3D(ManufacturedSolution):
+    """u(X) = [X1^4 + 2 X2 X3 / 5,
+               X2^4 + 2 X1 X3 / 5,
+               X3^4 / 10 - 2 X1 X2 X3 / 5]
+
+    u prescribed on the inner and outer curved surfaces, traction on the two flat ones.
+    """
+
+    geometry = Tire
+    dim = 3
+
+    def displacement(self, coordinates):
+        x, y, z = coordinates
+        X1, X2, X3 = x, z, -y
+
+        u1 = X1**4 + sp.Rational(2, 5) * X2 * X3
+        u2 = X2**4 + sp.Rational(2, 5) * X1 * X3
+        u3 = sp.Rational(1, 10) * X3**4 - sp.Rational(2, 5) * X1 * X2 * X3
+
+        return [self.amplitude * u1, -self.amplitude * u3, self.amplitude * u2]

From 990af596342588d7fdbcde5f396332259a92c0c6 Mon Sep 17 00:00:00 2001
From: Themis Skamagkis <themisskamagkis@gmail.com>
Date: Mon, 12 Oct 2026 09:49:18 +0200
Subject: [PATCH 21/55] [verification] Add the manufactured problem

---
 verification/manufactured/problem.py | 86 ++++++++++++++++++++++++++++
 1 file changed, 86 insertions(+)
 create mode 100644 verification/manufactured/problem.py

diff --git a/verification/manufactured/problem.py b/verification/manufactured/problem.py
new file mode 100644
index 0000000..cf44eb1
--- /dev/null
+++ b/verification/manufactured/problem.py
@@ -0,0 +1,86 @@
+"""Manufactured problems: a manufactured solution paired with a material, and everything derived
+from the pair.
+"""
+
+import numpy as np
+import sympy as sp
+
+from .materials import compile_laws
+
+_COORDINATES = sp.symbols("x y z")
+
+_COMPONENT_NAMES = ("u_x", "u_y", "u_z")
+
+
+class ManufacturedProblem:
+    """A manufactured solution paired with a material. This pairing determines displacement
+    gradient, stress, body-force source, and energy density.
+    """
+
+    def __init__(self, solution, material, spatial_dimensions):
+        self.solution = solution
+        self.material = material
+        dimensions = spatial_dimensions
+
+        self.coordinates = sp.Matrix(_COORDINATES[:dimensions])
+
+        # The solution states its own dim components; embedding into higher spatial dimensions pads
+        # the rest with zero.
+        components = list(self.solution.displacement(_COORDINATES[:self.solution.dim]))
+        if len(components) != self.solution.dim:
+            raise ValueError(f"{type(self.solution).__name__}: displacement must state "
+                             f"{self.solution.dim} components, one per dimension of the solution, "
+                             f"got {len(components)}")
+        components += [0] * (dimensions - self.solution.dim)
+
+        displacement = sp.Matrix(components)
+        stress = material.stress(displacement.jacobian(self.coordinates))
+
+        # Keep a symbolic representation of the displacement: `equation` reads off this directly.
+        self.displacement_expression = displacement
+        self.stress_expression = stress
+
+        gradient = displacement.jacobian(self.coordinates)
+        source = -sp.Matrix([sum(sp.diff(stress[i, j], self.coordinates[j])
+                                 for j in range(dimensions))
+                             for i in range(dimensions)])
+        self.source_expression = source
+
+        self.u =self._compile_field(displacement, (dimensions,))
+        self.grad_u = self._compile_field(gradient, (dimensions, dimensions))
+        self.stress = self._compile_field(stress, (dimensions, dimensions))
+        self.source = self._compile_field(source, (dimensions,))
+        self.energy_density, self.tangent = compile_laws(material, dimensions)
+
+    # The regions to apply BCs are further stated by its solution object.
+    @property
+    def prescribe_displacement_on(self):
+        return self.solution.prescribe_displacement_on
+
+    @property
+    def equation(self):
+        """The manufactured displacement in LaTeX form, one line per component."""
+        names = ("u",) if len(self.displacement_expression) == 1 else _COMPONENT_NAMES
+        return "\n".join(f"${name} = {sp.latex(component)}$"
+                         for name, component in zip(names, self.displacement_expression))
+
+    def _compile_field(self, expression, shape):
+        """Compile a symbolic expression of the coordinates into a numpy field function.
+
+        The result, `evaluate(points)`, takes one point or a whole batch of points (shape
+        (..., spatial_dimensions)) and returns the expression's value at each, shaped (*batch, *shape).
+        """
+        components = [sp.lambdify(list(self.coordinates), entry, "numpy") for entry in expression]
+
+        def compiled_field(points):
+            points = np.asarray(points, dtype=float)
+            batch = points.shape[:-1]
+            columns = [points[..., d] for d in range(len(self.coordinates))]
+            values = np.empty(batch + shape, dtype=float)
+            for index, component in zip(np.ndindex(*shape), components):
+                values[(...,) + index] = np.broadcast_to(
+                    np.asarray(component(*columns), dtype=float), batch)
+            return values
+
+        return compiled_field
+

From 24a33ff47cdc98d72a397edce526ad63c95a6441 Mon Sep 17 00:00:00 2001
From: Themis Skamagkis <themisskamagkis@gmail.com>
Date: Mon, 12 Oct 2026 09:49:18 +0200
Subject: [PATCH 22/55] [verification] Export the manufactured solutions and
 problem

---
 verification/manufactured/__init__.py | 5 +++++
 1 file changed, 5 insertions(+)
 create mode 100644 verification/manufactured/__init__.py

diff --git a/verification/manufactured/__init__.py b/verification/manufactured/__init__.py
new file mode 100644
index 0000000..44adc87
--- /dev/null
+++ b/verification/manufactured/__init__.py
@@ -0,0 +1,5 @@
+from .base import ManufacturedSolution, TrigonometricSolution
+from .problem import ManufacturedProblem
+from .grid import (Quadratic1D, Quadratic2D, Quadratic3D,
+                   Trigonometric1D, Trigonometric2D, Trigonometric3D)
+from .tire import Quartic2D, Quartic3D

From 22cbec5f88dd5ae639f7ac96b0aa6bcba8b09be8 Mon Sep 17 00:00:00 2001
From: Themis Skamagkis <themisskamagkis@gmail.com>
Date: Mon, 12 Oct 2026 09:49:18 +0200
Subject: [PATCH 23/55] [verification] Add the study package

---
 verification/study/__init__.py | 0
 1 file changed, 0 insertions(+), 0 deletions(-)
 create mode 100644 verification/study/__init__.py

diff --git a/verification/study/__init__.py b/verification/study/__init__.py
new file mode 100644
index 0000000..e69de29

From cbf383d5fe8f03d98748e48470b59a4c4d65b9b1 Mon Sep 17 00:00:00 2001
From: Themis Skamagkis <themisskamagkis@gmail.com>
Date: Mon, 12 Oct 2026 09:49:18 +0200
Subject: [PATCH 24/55] [verification] Add the mesh quadrature

---
 verification/study/quadrature.py | 58 ++++++++++++++++++++++++++++++++
 1 file changed, 58 insertions(+)
 create mode 100644 verification/study/quadrature.py

diff --git a/verification/study/quadrature.py b/verification/study/quadrature.py
new file mode 100644
index 0000000..ace160a
--- /dev/null
+++ b/verification/study/quadrature.py
@@ -0,0 +1,58 @@
+"""Dim-generic mesh quadrature: the fields sampled on it, and the norms integrated over them."""
+
+import numpy as np
+
+import Sofa.SofaFEM
+
+
+class MeshQuadrature:
+    """The reference->physical mapping of a whole mesh, evaluated once and shared by every norm.
+
+    Array layout throughout: `e` element, `q` quadrature point, `a` node within an element,
+    `c` field component, `d` spatial direction.
+    """
+
+    def __init__(self, nodes, node_indices, element, degree):
+        self.nodes = np.asarray(nodes)
+        # node_indices[e] = the mesh-node indices that form element e.
+        self.node_indices = np.asarray(node_indices)
+        dim = self.nodes.shape[1]
+
+        # Reference-space data is identical for every element of this type, so fetch it once.
+        reference_weights, self.shape_functions, reference_gradients = \
+            Sofa.SofaFEM.quadrature_data(element, dim, degree)
+        self.shape_function_gradients, measures = Sofa.SofaFEM.element_mapping_batch(
+            element, self.nodes, self.node_indices, reference_gradients)
+
+        self.points = np.einsum('qa,ead->eqd', self.shape_functions, self.nodes[self.node_indices])
+        # The dOmega of every integral: the reference weight carried into physical space.
+        self.weights = reference_weights * measures
+
+    def integrate(self, values):
+        """Integrate a scalar given at every (element, quadrature point)."""
+        return float(np.sum(values * self.weights))
+
+    def values(self, field):
+        """A nodal field at every quadrature point: (e, q, c)."""
+        return np.einsum('qa,eac->eqc', self.shape_functions, np.asarray(field)[self.node_indices])
+
+    def grads(self, field):
+        """Gradient of a nodal field: (e, q, c, d) = d(u_c)/dx_d."""
+        return np.einsum('eac,eqad->eqcd',
+                         np.asarray(field)[self.node_indices], self.shape_function_gradients)
+
+    def sample(self, field):
+        """A manufactured field at every quadrature point: the one seam where it meets this mesh."""
+        return field(self.points)
+
+    def L2Norm(self, values):
+        """L2 norm of a vector field given at every quadrature point: (e, q, c) -> scalar."""
+        return float(np.sqrt(self.integrate(np.sum(values * values, axis=-1))))
+
+    def FrobeniusNorm(self, tensors):
+        """L2 norm of a tensor field, Frobenius over its last two axes: (e, q, c, d) -> scalar."""
+        return float(np.sqrt(self.integrate(np.sum(tensors * tensors, axis=(-2, -1)))))
+
+    def energy(self, gradients, energy_density):
+        """Elastic energy of a displacement gradient field: integral over the mesh of psi."""
+        return self.integrate(energy_density(gradients))

From c651f9eb1a02ebaa3dd07341f03b99c0b3cff2ab Mon Sep 17 00:00:00 2001
From: Themis Skamagkis <themisskamagkis@gmail.com>
Date: Mon, 12 Oct 2026 09:49:18 +0200
Subject: [PATCH 25/55] [verification] Add the error metrics

---
 verification/study/metrics.py | 124 ++++++++++++++++++++++++++++++++++
 1 file changed, 124 insertions(+)
 create mode 100644 verification/study/metrics.py

diff --git a/verification/study/metrics.py b/verification/study/metrics.py
new file mode 100644
index 0000000..dfc24f4
--- /dev/null
+++ b/verification/study/metrics.py
@@ -0,0 +1,124 @@
+"""The quantities a level measures: one class per metric, so a metric is defined in one place."""
+
+from abc import ABC, abstractmethod
+
+import numpy as np
+
+
+class Measurement:
+    """One solved mesh refinement level, sampled at its quadrature points."""
+
+    def __init__(self, displacement, energy_h, manufactured_problem, quadrature):
+        # ManufacturedProblem Class
+        self.manufactured_problem = manufactured_problem
+
+        # MeshQuadrature Class
+        self.quadrature = quadrature
+
+        # The FEM displacement and its gradient, interpolated from the field SOFA solved for to the
+        # quadrature points.
+        self.u_h = quadrature.values(displacement)
+        self.grad_h = quadrature.grads(displacement)
+        # The strain energy SOFA reports for that same solution, from its own element stiffness.
+        self.energy_h = energy_h
+
+        # The manufactured solution's analytical displacement and gradient, evaluated at those same
+        # points.
+        self.u_exact = quadrature.sample(manufactured_problem.u)
+        self.grad_exact = quadrature.sample(manufactured_problem.grad_u)
+
+        # u_error/grad_error are what a metric measures.
+        # u_exact/grad_exact (not the error) set its noise floor, computed through the same norm
+        self.u_error = self.u_h - self.u_exact
+        self.grad_error = self.grad_h - self.grad_exact
+
+        self.energy_exact = quadrature.energy(self.grad_exact, manufactured_problem.energy_density)
+
+
+class Metric(ABC):
+    """A quantity measured at every refinement level: how to measure it, and what it is small relative to."""
+
+    name = ""
+
+    @abstractmethod
+    def measure(self, measurement):
+        """The metric's value on one level."""
+
+    @abstractmethod
+    def scale(self, measurement):
+        """Reference magnitude the noise floor (relativeNoiseFloor * this) is taken against:
+        the norms scale with the field, the energy metrics with the energy."""
+
+
+class L2(Metric):
+    """Root of the integrated squared error:
+
+        ‖e‖_L2 = √(∫ |e|² dΩ),   e = u_h - u
+    """
+
+    name = "L2"
+
+    def measure(self, measurement):
+        return measurement.quadrature.L2Norm(measurement.u_error)
+
+    def scale(self, measurement):
+        return measurement.quadrature.L2Norm(measurement.u_exact)
+
+
+class H1(Metric):
+    """H1 semi-norm of the error, the Frobenius norm of its gradient:
+
+        ‖e‖_H1 = √(∫ ‖∇e‖_F² dΩ),   e = u_h - u
+    """
+
+    name = "H1"
+
+    def measure(self, measurement):
+        return measurement.quadrature.FrobeniusNorm(measurement.grad_error)
+
+    def scale(self, measurement):
+        return measurement.quadrature.FrobeniusNorm(measurement.grad_exact)
+
+
+class ErrorEnergy(Metric):
+    """The strain energy the error field carries:
+
+        √(2 ∫ ψ(∇e) dΩ),   e = u_h - u
+
+    SOFA supplies the displacement; ψ and the integration are python's. Left out of METRICS: it is
+    H1 under a material weighting, so it restates that column instead of testing anything further.
+    """
+
+    name = "ErrorEnergy"
+
+    def measure(self, measurement):
+        return float(np.sqrt(2.0 * measurement.quadrature.energy(
+            measurement.grad_error, measurement.manufactured_problem.energy_density)))
+
+    def scale(self, measurement):
+        return np.sqrt(2.0 * measurement.energy_exact)
+
+
+class EnergyDifference(Metric):
+    """How far the strain energy SOFA reports sits from the one the exact solution carries:
+
+        √(2 |E(u) - E(u_h)|),   E(u) integrated from ψ here, E(u_h) SOFA's getPotentialEnergy
+    """
+
+    name = "EnergyDifference"
+
+    def measure(self, measurement):
+        return float(np.sqrt(2.0 * abs(measurement.energy_exact - measurement.energy_h)))
+
+    def scale(self, measurement):
+        return np.sqrt(2.0 * measurement.energy_exact)
+
+
+# METRICS is the fixed list of metric instances every level is measured against. Order matters
+# because nothing else fixes it: a study iterates over this list to set the column order of every
+# printed table and the key order of every JSON record, so reordering it here reorders both.
+#
+# Expected orders are the case's `expectedOrder`, not a constant here: they follow from the element,
+# not the metric, so a P2 case states different numbers for the same metrics (P1: L2 -> 2,
+# H1 -> 1, ErrorEnergy -> 1, EnergyDifference -> 1).
+METRICS = [L2(), H1(), EnergyDifference()]

From 7de190da0b3330de3d4e431487ce0aa32c18d99e Mon Sep 17 00:00:00 2001
From: Themis Skamagkis <themisskamagkis@gmail.com>
Date: Mon, 12 Oct 2026 09:49:18 +0200
Subject: [PATCH 26/55] [verification] Add the scene package

---
 verification/scene/__init__.py | 0
 1 file changed, 0 insertions(+), 0 deletions(-)
 create mode 100644 verification/scene/__init__.py

diff --git a/verification/scene/__init__.py b/verification/scene/__init__.py
new file mode 100644
index 0000000..e69de29

From ae01dd8ce86b31287ad7b2710a55165f0f29ded8 Mon Sep 17 00:00:00 2001
From: Themis Skamagkis <themisskamagkis@gmail.com>
Date: Mon, 12 Oct 2026 09:49:18 +0200
Subject: [PATCH 27/55] [verification] Add the manufactured source and traction
 controllers

---
 verification/scene/controllers.py | 109 ++++++++++++++++++++++++++++++
 1 file changed, 109 insertions(+)
 create mode 100644 verification/scene/controllers.py

diff --git a/verification/scene/controllers.py b/verification/scene/controllers.py
new file mode 100644
index 0000000..8893f5b
--- /dev/null
+++ b/verification/scene/controllers.py
@@ -0,0 +1,109 @@
+"""Post-init SOFA controllers loading a mesh with a manufactured problem."""
+
+import numpy as np
+
+import Sofa
+import Sofa.Core
+
+from common.sofa.conventions import ELEMENTS
+
+
+class ApplyManufacturedSourceTerm(Sofa.Core.Controller):
+    """Loads a mesh with the body force a manufactured solution puts on the right-hand side."""
+
+    def __init__(self, node, dofs, source, vec_type, element, source_quadrature_degree,
+                 *args, **kwargs):
+        super().__init__(*args, **kwargs)
+        self.node = node
+        self.dofs = dofs
+        self.source = source
+        self.vec_type = vec_type
+        self.element = element
+        self.source_quadrature_degree = source_quadrature_degree
+
+    def init(self):
+        template = f'{self.vec_type},{ELEMENTS[self.element].cpp}'
+
+        self.density = self.node.addObject('NodalSourceDensity', name='sourceDensity',
+                                           template=self.vec_type)
+        # One call for the whole mesh: the manufactured fields evaluate over an array of points.
+        self.density.property.value = self.source(self.dofs.rest_position.array())
+
+        self.node.addObject('VectorSourceTerm', name='bodyForce', template=template,
+                            sourceDensity='@sourceDensity')
+        self.node.addObject('FEMSourceTermIntegrator', name='bodySource', template=template,
+                            topology='@topology', quadratureDegree=self.source_quadrature_degree,
+                            constantSources='@bodyForce')
+
+
+class ApplyManufacturedTraction(Sofa.Core.Controller):
+    """Loads the whole boundary of a mesh with the traction a manufactured solution exerts on it."""
+
+    def __init__(self, geometry, node, dofs, stress, vec_type, element, spatial_dimensions,
+                 source_quadrature_degree, *args, **kwargs):
+        super().__init__(*args, **kwargs)
+        self.geometry = geometry
+        self.node = node
+        self.dofs = dofs
+        self.stress = stress
+        self.vec_type = vec_type
+        self.element = element
+        self.spatial_dimensions = spatial_dimensions
+        self.source_quadrature_degree = source_quadrature_degree
+        self.nodal_stress = None
+        # Point load is used for the particular case of Traction BC in 1D
+        self.point_load = None
+
+    def init(self):
+        rest_positions = self.dofs.rest_position.array()
+        mappings = ELEMENTS[self.element].boundary_mappings
+
+        # Only Edges have no Volume2Boundary mapping. Simply prescribe a nodal force.
+        if not mappings:
+            self.add_point_load(rest_positions)
+            return
+
+        boundary_kind = ELEMENTS[mappings[-1][0]]
+        if boundary_kind.dim + 1 != self.spatial_dimensions:
+            # TODO StressSourceTerm does not support elements of codimension 2. So we cannot apply
+            # traction BCs on Quad/Triangular meshes in 3D. The Edge normals are arbitrary.
+            return
+
+        # Goes from the volume to the boundary elements, applying one topological mapping at a time.
+        boundary = self.node
+        for kind, mapping in mappings:
+            element_kind = ELEMENTS[kind]
+            boundary = boundary.addChild(kind)
+            boundary.addObject(element_kind.container, name='topology')
+            boundary.addObject(element_kind.container.replace('Container', 'Modifier'))
+            boundary.addObject(mapping, input='@../topology', output='@topology')
+
+        # Apply the nodal stress load on the boundary.
+        template = f'{self.vec_type},{boundary_kind.cpp}'
+        # The manufactured stress at every node, converted to the layout NodalStress takes.
+        stress_values = self.stress(rest_positions)
+        rows, columns = np.tril_indices(stress_values.shape[-1])
+        self.nodal_stress = boundary.addObject('NodalStress', name='stress', template=self.vec_type)
+        self.nodal_stress.property.value = np.ascontiguousarray(stress_values[:, rows, columns])
+
+        boundary.addObject('StressSourceTerm', name='traction', template=template, stress='@stress')
+        boundary.addObject('FEMSourceTermIntegrator', name='tractionSource', template=template,
+                           topology='@topology', quadratureDegree=self.source_quadrature_degree,
+                           constantSources='@traction')
+
+    def add_point_load(self, rest_positions):
+        """The nodal traction on every boundary region, its normal pointing away from the mesh."""
+        centre = rest_positions.mean(axis=0)
+        indices, tractions = [], []
+        for region in self.geometry.region_names:
+            nodes = self.geometry.region_indices(region, rest_positions)
+            outward = rest_positions[nodes] - centre
+            outward /= np.linalg.norm(outward, axis=-1, keepdims=True)
+
+            indices += nodes
+            tractions.append(np.einsum('nij,nj->ni',
+                                       self.stress(rest_positions[nodes]), outward))
+
+        self.point_load = self.node.addObject('ConstantForceField', name='traction',
+                                              template=self.vec_type, indices=indices,
+                                              forces=np.concatenate(tractions))

From 37b9a4ea80bbbccef8e56f63f2e53a3ff7507353 Mon Sep 17 00:00:00 2001
From: Themis Skamagkis <themisskamagkis@gmail.com>
Date: Mon, 12 Oct 2026 09:49:19 +0200
Subject: [PATCH 28/55] [verification] Add the MMS scene

---
 verification/scene/scene.py | 89 +++++++++++++++++++++++++++++++++++++
 1 file changed, 89 insertions(+)
 create mode 100644 verification/scene/scene.py

diff --git a/verification/scene/scene.py b/verification/scene/scene.py
new file mode 100644
index 0000000..1237e82
--- /dev/null
+++ b/verification/scene/scene.py
@@ -0,0 +1,89 @@
+"""Verification suite: an MMS scene whose BCs come from a manufactured solution."""
+
+from common.sofa import SofaScene
+from common.sofa.conventions import VEC_DIM
+from common.sofa.controllers import RegionClamp
+from common.sofa.prefabs import PREFAB_BY_GEOMETRY
+
+from .controllers import ApplyManufacturedSourceTerm, ApplyManufacturedTraction
+
+
+class MMSScene(SofaScene):
+    """A scene set up to solve a modified PDE. Boundary conditions and a body-force load are
+    chosen so that the solution of that modified problem is a pre-manufactured one."""
+
+    plugins = [
+        "Sofa.Component.SolidMechanics.FEM.Elastic",
+        "Sofa.Component.Constraint.Projective",
+        "Sofa.Component.Engine.Select",
+        "Sofa.Component.Mapping.Linear",
+        "Sofa.Component.MechanicalLoad",
+        "Sofa.Component.LinearSolver.Direct",
+        "Sofa.Component.LinearSolver.Iterative",
+        "Sofa.Component.LinearSolver.Preconditioner",
+        "Sofa.Component.LinearSystem",
+        "Sofa.Component.ODESolver.Backward",
+        "Sofa.Component.StateContainer",
+        "Sofa.Component.Topology.Container.Grid",
+        "Sofa.Component.Topology.Container.Dynamic",
+        "Sofa.Component.Topology.Mapping",
+    ]
+
+    def __init__(self, case, cells):
+        """One mesh of the case, `cells` counting its cells per axis."""
+        super().__init__()
+        self.cells = cells
+        self.geometry = case.geometry
+        self.material = case.material
+        self.force_field = case.force_field
+        self.element = case.element
+        self.solvers = case.solvers
+        self.manufactured_problem = case.manufactured_problem
+        self.source_quadrature_degree = case.source_quadrature_degree
+        self.grid_mapping = case.grid_mapping
+
+    @property
+    def grid_resolution(self):
+        """The grid resolution of `RegularGridTopology`, which is its `n` Data."""
+        return [count + 1 for count in self.cells]
+
+    def body(self, root):
+        root.addObject('DefaultAnimationLoop')
+
+        prefab_cls = PREFAB_BY_GEOMETRY[type(self.geometry)]
+        configs = {'material': self.material, 'forceField': self.force_field,
+                   'solvers': self.solvers, 'gridMapping': self.grid_mapping}
+        prefab = root.addChild(prefab_cls(name=prefab_cls.__name__, geometry=self.geometry,
+                                          grid_resolution=self.grid_resolution, element=self.element, configs=configs))
+        mechanical = prefab.mechanical_node()
+
+        dim, spatial_dimensions = self.geometry.dim, self.geometry.spatial_dimensions
+        VEC = VEC_DIM[spatial_dimensions]
+
+        # Load the mesh with the body force the manufactured solution puts on the rhs
+        mechanical.addObject(ApplyManufacturedSourceTerm(
+            node=mechanical, dofs=mechanical.dofs, source=self.manufactured_problem.source,
+            vec_type=VEC, element=self.element,
+            source_quadrature_degree=self.source_quadrature_degree, name='bodyForceCtrl'))
+
+        # Load the whole boundary with the traction that solution exerts there; on the directions
+        # the case left clamped it is absorbed by the constraint.
+        mechanical.addObject(ApplyManufacturedTraction(
+            geometry=self.geometry, node=mechanical, dofs=mechanical.dofs,
+            stress=self.manufactured_problem.stress,
+            vec_type=VEC, element=self.element, spatial_dimensions=spatial_dimensions,
+            source_quadrature_degree=self.source_quadrature_degree, name='tractionCtrl'))
+
+        # Clamp the mesh dofs by region as determined by the manufactured solution BCs
+        mechanical.addObject(RegionClamp(
+            geometry=self.geometry, dofs=mechanical.dofs, node=mechanical,
+            prescribe_displacement_on=self.manufactured_problem.prescribe_displacement_on,
+            vec_type=VEC, displacement=self.manufactured_problem.u, name='clampCtrl'))
+
+        # For an embedded mesh (dim < spatial_dimensions) the out-of-plane dimensions are fixed.
+        if dim < spatial_dimensions:
+            mechanical.addObject('PartialFixedProjectiveConstraint', name='outOfPlane', template=VEC,
+                           fixAll=True,
+                           fixedDirections=[0] * dim + [1] * (spatial_dimensions - dim)) # e.g. [0, 0] + [1] -> [0, 0, 1]
+
+        return prefab

From 019adc4ee764119d538dde1526b452f0dc60c04c Mon Sep 17 00:00:00 2001
From: Themis Skamagkis <themisskamagkis@gmail.com>
Date: Mon, 12 Oct 2026 09:49:19 +0200
Subject: [PATCH 29/55] [verification] Add the error-convergence study

---
 verification/study/study.py | 549 ++++++++++++++++++++++++++++++++++++
 1 file changed, 549 insertions(+)
 create mode 100644 verification/study/study.py

diff --git a/verification/study/study.py b/verification/study/study.py
new file mode 100644
index 0000000..7441a97
--- /dev/null
+++ b/verification/study/study.py
@@ -0,0 +1,549 @@
+"""Runs an error-convergence-rate verification study in SOFA to exercise FEM Force Fields:
+
+1. Solve a manufactured problem on a sequence of increasingly refined meshes.
+2. Measure the discretisation error at each one.
+3. Judge whether the observed order of accuracy matches what the discretisation should achieve.
+
+References
+----------
+Roache, P. J. (1998). Verification and Validation in Computational Science and Engineering.
+    Hermosa Publishers.
+Roy, C. J. (2005). Review of Code and Solution Verification Procedures for Computational
+    Simulation. Journal of Computational Physics, 205(1), 131-156.
+Stern, F., Wilson, R. V., Coleman, H. W. and Paterson, E. G. (2001). Comprehensive Approach to
+    Verification and Validation of CFD Simulations -- Part 1: Methodology and Procedures. Journal
+    of Fluids Engineering, 123(4), 793-802.
+"""
+
+import json
+import pathlib
+import shutil
+import sys
+from dataclasses import dataclass, field
+from enum import StrEnum
+
+import matplotlib
+matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+from matplotlib import colormaps
+from matplotlib.collections import PolyCollection
+from matplotlib.colors import Normalize, to_hex
+from matplotlib.lines import Line2D
+from matplotlib.ticker import NullLocator
+from mpl_toolkits.mplot3d.art3d import Poly3DCollection
+import numpy as np
+
+import Sofa.Core
+import Sofa.Simulation
+
+from .metrics import Measurement, METRICS
+from .quadrature import MeshQuadrature
+from ..scene.scene import MMSScene
+from common.sofa.conventions import ELEMENTS
+
+
+def _newton_diagnostics(newton):
+    """Iterations, residual reduction and stopping status of the Newton solve at one level."""
+    if newton is None:
+        return None
+
+    residuals = []
+    for token in newton.residualGraph.value.split():
+        try:
+            residuals.append(float(token))
+        except ValueError:      # the map key, not one of its residuals
+            continue
+
+    status = newton.status.value
+    # ConvergedEquilibrium means the iterations never started because the residual was already zero.
+    # On a loaded MMS scene that is not a success, it means the load never reached the system.
+    converged = status.startswith("Converged") and status != "ConvergedEquilibrium"
+
+    # Dimensionless reduction r/r0 of the residual. Undefined when it started at zero.
+    reduction = float(np.sqrt(residuals[-1] / residuals[0])) if residuals and residuals[0] else float("nan")
+
+    # Ignore the residual recorded once before iterating; floor at 0 if a crashed solver left none.
+    return {"iterations": max(len(residuals) - 1, 0),
+            "reduction": reduction,
+            "status": "Converged" if converged else status}
+
+
+# Where a run writes what it produced.
+RESULTS_ROOT = pathlib.Path(__file__).parent.parent / "results"
+
+
+def clear_results():
+    """Delete everything previous runs wrote under the results root."""
+    shutil.rmtree(RESULTS_ROOT, ignore_errors=True)
+
+SURFACE, INK, INK_SOFT, GRID = "#fcfcfb", "#0b0b0b", "#52514e", "#e8e7e3"
+PLOT_ACCEPTED, PLOT_DISCARDED = to_hex(colormaps["Greens"](0.95)), to_hex(colormaps["Reds"](0.95))
+METRIC_FILL = {"L2": INK, "H1": to_hex(colormaps["Blues"](0.72)),
+               "EnergyDifference": SURFACE}
+ELEMENT_MARKER = {"edge": "o", "quad": "s", "hexa": "s", "tri": "^", "tet": "^"}
+
+
+# Column widths shared by every table printed here.
+LABEL_WIDTH, SPACING_WIDTH, ELEMENTS_WIDTH = 16, 20, 17
+RATE_WIDTH, ITERATIONS_WIDTH, RESIDUAL_WIDTH = 16, 14, 12
+
+def _paint(cell, accepted):
+    """Colour a cell green when the settling test kept its rate, red when it discarded it."""
+    GREEN, RED, RESET = "\033[32m", "\033[31m", "\033[0m"
+
+    if not cell.strip() or not sys.stdout.isatty():
+        return cell
+    return f"{GREEN if accepted else RED}{cell}{RESET}"
+
+
+def _mesh_spacing(quadrature, dim):
+    """This mesh's step size h: the representative element size, the `dim`-th root of the mean
+    element measure. `dim` is the element's own dimension, so an embedded mesh divides its area,
+    not a volume.
+
+    A mean, not the largest element diameter: the error norms are global, so one worst element
+    should not set the h they are plotted against, and on a swept mesh a diameter does not halve
+    cleanly under refinement even when every cell does. The refinement factor it yields between
+    two levels is (N_coarse / N_fine) raised to 1/dim, which is what the convergence literature
+    uses where there is no single grid spacing to quote. Measured off the built mesh, so a loaded
+    one needs no cell counts to state it.
+    """
+    return float((np.sum(quadrature.weights) / len(quadrature.node_indices)) ** (1.0 / dim))
+
+
+def _mesh_label(cells):
+    """This mesh's cell count per axis as text, e.g. [24, 24] -> "24x24"."""
+    return 'x'.join(str(count) for count in cells)
+
+
+def _print_progress_bar(level, level_count, label):
+    """Print a bar showing how much of the run is done and which mesh it is on."""
+    BAR_WIDTH, LABEL_WIDTH = 20, 14
+
+    if not sys.stdout.isatty():
+        return
+
+    filled = round(BAR_WIDTH * level / level_count)
+    bar = "=" * filled + "-" * (BAR_WIDTH - filled)
+    progress = f"  {level}/{level_count} [{bar}] {label:<{LABEL_WIDTH}}"
+
+    sys.stdout.write(f"\r{progress}\033[K")
+    sys.stdout.flush()
+
+
+def _clear_progress_bar():
+    """Erase the progress bar so the table below starts on a clean line."""
+    if sys.stdout.isatty():
+        sys.stdout.write("\r\033[K")
+        sys.stdout.flush()
+
+
+@dataclass
+class ConvergenceRate:
+    """The order observed between one mesh refinement level and the one below it, or why there is none."""
+
+    class FailureReason(StrEnum):
+        """Why computing a rate for a level pair failed, printed in its place when there is no rate."""
+        FLOOR = "floor"   # one of the two errors has sunk into its own noise; the ratio is round-off, not order
+        OSC = "osc"       # the error overshot then came back down (Stern, -1 < R < 0)
+        DIV = "div"       # the error is diverging, steadily or through growing alternating increments (Stern, |R| >= 1)
+
+    value: float | None
+    reason: FailureReason | None = None
+
+
+@dataclass
+class MeshRefinementLevelResult:
+    """One mesh of the refinement sequence: what was measured on it, and how its solve went."""
+
+    label: str                                   # this mesh's name in the tables, e.g. "24x24"
+    spacing: float                                # this mesh's step size, the h in the order formula
+    elements_count: int                          # element count the mesh mapping actually produced
+    values: dict                                 # metric name -> its measured error on this mesh
+    noise_floor: dict                            # metric name -> the magnitude below which its error is noise
+    newton: dict | None                          # this mesh's Newton solve diagnostics, or None if none ran
+    rates: dict = field(default_factory=dict)    # metric name -> ConvergenceRate, unknown until this level is added
+    # The solved field, kept for drawing.
+    nodes: np.ndarray | None = None              # rest positions of the mesh nodes
+    elements: np.ndarray | None = None           # node indices of each element
+    displacement: np.ndarray | None = None       # u_h at the nodes
+    boundary: np.ndarray | None = None           # node indices of each boundary face, 3D only
+
+
+@dataclass
+class MetricResult:
+    """What the mesh refinement sequence found for one metric, once every level had been solved."""
+
+    order: float | None     # the order at the finest settled mesh, or None if the order never settled
+    accepted_levels: set    # level indices whose rate falls inside the settled range
+    passed: bool            # whether `order` matches the case's expected order, within tolerance
+
+
+class ErrorConvergenceStudy:
+    """Performs an error-convergence-rate study for a case: builds a SOFA scene and solves it for
+    progressively refined meshes.
+    """
+
+    def __init__(self, case):
+        self.name = case.name
+        self.case = case
+        self.asymptotic_tolerance = case.asymptotic_tolerance
+        self.expected_order = case.expected_order
+        self.expected_order_tolerance = case.expected_order_tolerance
+        self.levels = []
+        self.results = {}
+        self.unconverged_levels = []
+        self.verified = False
+        self.plots = {}
+
+    def run(self):
+        """Solve every mesh in the case's refinement sequence, then evaluate and report the results."""
+        mesh_refinement_levels = self.case.levels()
+        try:
+            for index, cells in enumerate(mesh_refinement_levels, start=1):
+                _print_progress_bar(index, len(mesh_refinement_levels), _mesh_label(cells))
+                self.solve(cells)
+        finally:
+            _clear_progress_bar()
+        self.evaluate_metrics()
+        self.plot()
+        self.report()
+
+    def solve(self, cells):
+        """Build and solve one mesh's SOFA scene, then measure and add it as the next level.
+
+        `cells` is that mesh's cell count per axis.
+        """
+        root = Sofa.Core.Node("root")
+        mechanical = MMSScene(self.case, cells).build(root).mechanical_node()
+        Sofa.Simulation.init(root)
+        Sofa.Simulation.animate(root, root.dt.value)
+
+        self.measure(mechanical, _mesh_label(cells))
+
+        Sofa.Simulation.unload(root)
+
+    def measure(self, mechanical, label):
+        """Measure an already-solved mesh, then add it as the next level.
+
+        `mechanical` is that mesh's mechanical node and `label` names it.
+        """
+        element_kind = ELEMENTS[self.case.element]
+        nodes = mechanical.dofs.rest_position.array()
+        u_h = mechanical.dofs.position.array() - nodes
+        # node_indices[e] = the mesh-node indices forming element e (from the topology).
+        node_indices = getattr(mechanical.topology, element_kind.data_name).array()
+
+        quadrature = MeshQuadrature(nodes, node_indices, element_kind.cpp, self.case.error_quadrature_degree)
+        measurement = Measurement(u_h, mechanical.fem.getPotentialEnergy(),
+                                  self.case.manufactured_problem, quadrature)
+
+        level = MeshRefinementLevelResult(
+            label=label, spacing=_mesh_spacing(quadrature, element_kind.dim),
+            # The element count the mapping actually produced, since the case states cells: one row
+            # of node_indices per element, so 6x the cells for tetrahedra and 2x for triangles.
+            elements_count=len(node_indices),
+            values={metric.name: metric.measure(measurement) for metric in METRICS},
+            noise_floor={metric.name: self.case.relative_noise_floor * metric.scale(measurement)
+                         for metric in METRICS},
+            newton=_newton_diagnostics(getattr(mechanical, 'newton', None)),
+            nodes=np.array(nodes), elements=np.array(node_indices), displacement=np.array(u_h))
+        if element_kind.dim == 3:
+            # The faces the traction controller mapped out of the mesh.
+            kind, _ = element_kind.boundary_mappings[-1]
+            level.boundary = np.array(getattr(mechanical.getChild(kind).topology, ELEMENTS[kind].data_name).array())
+        self.levels.append(level)
+        level.rates = {metric.name: self.rate(metric.name) for metric in METRICS}
+
+    def rate(self, metric_name):
+        """The order across the two finest meshes so far, or None and why there is none."""
+        if len(self.levels) < 2:
+            return ConvergenceRate(None)
+
+        current, previous = self.levels[-1], self.levels[-2]
+        error, before = current.values[metric_name], previous.values[metric_name]
+        # The discretization error should not be measured when below the round-off and/or iterative
+        # solver error (Roy, 2005).
+        if error <= current.noise_floor[metric_name] or before <= previous.noise_floor[metric_name]:
+            return ConvergenceRate(None, ConvergenceRate.FailureReason.FLOOR)
+
+        if len(self.levels) >= 3:
+            # Stern's convergence ratio is read from the increments between three consecutive errors.
+            # This also catches a sequence that oscillates while it is still decreasing on this pair.
+            #
+            #   R = (e_i - e_i-1) / (e_i-1 - e_i-2)
+            #   e_i, e_i-1, e_i-2: the error on this level and the two coarser ones before it
+            #   |R| >= 1:   divergence: steadily if R > 0, alternating in sign if R < 0
+            #   -1 < R < 0: still converging, but oscillating
+            #   0 <= R < 1: the normal case, falls through to the order formula below
+            increment = error - before
+            previous_increment = before - self.levels[-3].values[metric_name]
+            ratio = increment / previous_increment if previous_increment else np.inf
+            if abs(ratio) >= 1.0:
+                return ConvergenceRate(None, ConvergenceRate.FailureReason.DIV)
+            if ratio <= 0.0:
+                return ConvergenceRate(None, ConvergenceRate.FailureReason.OSC)
+        elif error >= before:
+            # With 2 levels only the most that can be said is whether the error increased
+            return ConvergenceRate(None, ConvergenceRate.FailureReason.DIV)
+
+        # float(): the numpy scalar would reach the record, which json cannot write.
+        return ConvergenceRate(float(np.log(error / before)
+                                     / np.log(current.spacing / previous.spacing)))
+
+    def evaluate_metrics(self):
+        """Judge the settled order per metric, then conclude self.verified: did this case verify?"""
+        self.results = {}
+        for metric in METRICS:
+            orders = [(index, level.rates[metric.name].value) for index, level in enumerate(self.levels)
+                      if level.rates[metric.name].value is not None]
+
+            # Contiguity matters: a readable order with a discarded pair beneath it does not extend a run.
+            retained = []
+            for index, order in reversed(orders):
+                if retained and (index + 1 != retained[0][0]
+                                 or abs(order - retained[0][1]) > self.asymptotic_tolerance):
+                    break
+                retained.insert(0, (index, order))
+            # One order cannot demonstrate that anything settled -- it takes two to compare.
+            if len(retained) < 2:
+                retained = []
+
+            if retained:
+                order = retained[-1][1]
+                accepted_levels = {index for index, _ in retained}
+                # bool(): the comparison is on a numpy float, so it yields a numpy bool
+                passed = bool(abs(order - self.expected_order[metric.name]) <= self.expected_order_tolerance)
+            else:
+                order, accepted_levels, passed = None, set(), False
+
+            self.results[metric.name] = MetricResult(order=order, accepted_levels=accepted_levels,
+                                                      passed=passed)
+
+        self.unconverged_levels = [level.label for level in self.levels
+                                   if (level.newton["status"] if level.newton else None)
+                                   not in ("Converged", None)]
+
+        # A settled, passing order is not evidence of anything if the solve never converged.
+        self.verified = not self.unconverged_levels and all(result.passed for result in self.results.values())
+
+    def write_results(self):
+        """Save the figures and the record under <results root>/<case path>."""
+        path = RESULTS_ROOT / self.name
+        path.parent.mkdir(parents=True, exist_ok=True)
+
+        for name, figure in self.plots.items():
+            figure.savefig(f"{path}_{name}.png")
+            plt.close(figure)
+
+        pathlib.Path(f"{path}.json").write_text(json.dumps(self.record(), indent=2))
+
+    def record(self):
+        """Record of what this study measured."""
+        return {"case": self.name,
+                "forceField": self.case.force_field["type"],
+                "order": {metric.name: self.results[metric.name].order for metric in METRICS},
+                "levels": [{"label": level.label,
+                            "spacing": level.spacing,
+                            "elements": level.elements_count,
+                            "values": level.values,
+                            "rates": {metric.name: level.rates[metric.name].value
+                                      for metric in METRICS}}
+                           for level in self.levels]}
+
+    def plot(self, axes=None):
+        """Build the convergence-rate figure into self.plots: error against mesh spacing, log-log,
+        one line per metric, its segments coloured by the verdict on the rate they carry.
+
+        Draws on `axes` when given, e.g. one panel of a side-by-side comparison, else on a figure
+        of its own. Returns the axes drawn on.
+        """
+        spacings = [level.spacing for level in self.levels]
+
+        if axes is None:
+            figure, axes = plt.subplots(figsize=(7.2, 4.8))
+        else:
+            figure = axes.figure
+        figure.patch.set_facecolor(SURFACE)
+        axes.set_facecolor(SURFACE)
+        axes.grid(True, color=GRID, linewidth=0.8)
+        axes.set_axisbelow(True)
+        for side in ("top", "right"):
+            axes.spines[side].set_visible(False)
+        for side in ("left", "bottom"):
+            axes.spines[side].set_color(GRID)
+        axes.tick_params(colors=INK_SOFT, labelsize=9)
+        axes.tick_params(which="minor", length=0)
+        axes.set_xscale("log")
+        axes.set_yscale("log")
+        axes.set_title(self.case.equation, fontsize=11, color=INK, pad=12)
+        axes.text(0.02, 0.97, f"Element: {self.case.element.capitalize()}", transform=axes.transAxes,
+                  ha="left", va="top", fontsize=11, color=INK)
+        axes.set_xlabel("mesh spacing h", fontsize=9, color=INK_SOFT)
+        axes.set_ylabel("Norm", fontsize=9, color=INK_SOFT)
+        axes.set_xlim(min(spacings) / 1.2, max(spacings) * 2.4)
+        axes.set_xticks(spacings)
+        axes.set_xticklabels([f"{spacing:.3g}" for spacing in spacings])
+        axes.xaxis.set_minor_locator(NullLocator())
+
+        marker = ELEMENT_MARKER.get(self.case.element, "D")
+        any_discarded = False
+        for metric in METRICS:
+            errors = [level.values[metric.name] for level in self.levels]
+            accepted_levels = self.results[metric.name].accepted_levels
+
+            for index in range(1, len(self.levels)):
+                kept = index - 1 in accepted_levels and index in accepted_levels
+                any_discarded |= not kept
+                axes.plot(spacings[index - 1:index + 1], errors[index - 1:index + 1],
+                          color=PLOT_ACCEPTED if kept else PLOT_DISCARDED, linewidth=2.0,
+                          solid_capstyle="round", zorder=3)
+            axes.plot(spacings, errors, linestyle="none", marker=marker, markersize=6,
+                      markerfacecolor=METRIC_FILL.get(metric.name, INK), markeredgecolor=INK_SOFT,
+                      markeredgewidth=1.2, zorder=4)
+            axes.annotate(metric.name, (spacings[0], errors[0]), textcoords="offset points",
+                          xytext=(9, 0), ha="left", va="center", fontsize=9, color=INK)
+
+            if len(self.levels) >= 2:
+                rate = self.levels[-1].rates[metric.name]
+                text = f"{rate.value:.2f}" if rate.value is not None else (rate.reason or "")
+                if text:
+                    anchor = ((spacings[-2] * spacings[-1]) ** 0.5, (errors[-2] * errors[-1]) ** 0.5)
+                    verdict = PLOT_ACCEPTED if len(self.levels) - 1 in accepted_levels else PLOT_DISCARDED
+                    axes.annotate(text, anchor, textcoords="offset points", xytext=(6, -8),
+                                  ha="left", va="top", fontsize=9, color=verdict)
+
+        if any_discarded:
+            axes.plot([], [], color=PLOT_DISCARDED, linewidth=2.0, label="pre-asymptotic regime")
+            axes.legend(loc="lower right", frameon=False, fontsize=9, handlelength=1.6,
+                        labelcolor=INK_SOFT)
+
+        figure.tight_layout()
+        self.plots = {"convergence": figure}
+        return axes
+
+    def plot_fields(self):
+        """Construct displacement field plots, shown on the deformed mesh, and store into self.plots"""
+        dim = self.case.geometry.dim
+        if dim == 1:
+            figure, axes = plt.subplots(figsize=(7.2, 4.2))
+            bar = np.linspace(self.levels[0].nodes[:, 0].min(), self.levels[0].nodes[:, 0].max(), 400)
+            samples = np.zeros((bar.size, self.levels[0].nodes.shape[1]))
+            samples[:, 0] = bar
+            axes.plot(bar, self.case.manufactured_problem.u(samples)[:, 0], color=INK,
+                      linewidth=2.4, label="exact")
+            colours = colormaps["viridis"](np.linspace(0.0, 0.85, len(self.levels)))
+            for level, colour in zip(self.levels, colours):
+                order = np.argsort(level.nodes[:, 0])
+                # Node markers on fine levels would merge into a thick band.
+                axes.plot(level.nodes[order, 0], level.displacement[order, 0], color=to_hex(colour),
+                          marker="o" if len(order) <= 25 else None, markersize=4,
+                          linewidth=1.4, label=level.label)
+            axes.set_facecolor(SURFACE)
+            axes.grid(True, color=GRID, linewidth=0.8)
+            axes.set_xlabel("x", fontsize=9, color=INK_SOFT)
+            axes.set_ylabel("u", fontsize=9, color=INK_SOFT)
+            axes.tick_params(colors=INK_SOFT, labelsize=9)
+            axes.legend(frameon=False, fontsize=8, labelcolor=INK_SOFT)
+            figure.patch.set_facecolor(SURFACE)
+            figure.tight_layout()
+        else:
+            figure, panels = plt.subplots(1, len(self.levels), figsize=(2.6 * len(self.levels), 3.0),
+                                          squeeze=False, layout="constrained",
+                                          subplot_kw={"projection": "3d"} if dim == 3 else {})
+            magnitudes = [np.linalg.norm(level.displacement, axis=-1) for level in self.levels]
+            # The finest level sets the scale: coarse levels overshoot and would stretch it.
+            norm = Normalize(vmin=magnitudes[-1].min(), vmax=magnitudes[-1].max())
+            for axes, level, magnitude in zip(panels[0], self.levels, magnitudes):
+                deformed = level.nodes + level.displacement
+                # Each face takes the mean |u_h| of its nodes.
+                if dim == 2:
+                    faces = PolyCollection(deformed[level.elements][..., :2],
+                                           array=magnitude[level.elements].mean(axis=1), cmap="viridis",
+                                           norm=norm, edgecolors="none", antialiased=False)
+                    axes.add_collection(faces)
+                    axes.autoscale_view()
+                    axes.set_aspect("equal")
+                else:
+                    faces = Poly3DCollection(deformed[level.boundary],
+                                             array=magnitude[level.boundary].mean(axis=1), cmap="viridis",
+                                             norm=norm, linewidths=0, antialiased=False)
+                    axes.add_collection3d(faces)
+                    low, high = deformed.min(axis=0), deformed.max(axis=0)
+                    axes.set(xlim=(low[0], high[0]), ylim=(low[1], high[1]), zlim=(low[2], high[2]))
+                    axes.set_box_aspect(high - low)
+                axes.set_axis_off()
+                axes.set_title(level.label, fontsize=9, color=INK)
+            figure.colorbar(faces, ax=panels[0].tolist(), shrink=0.8, label="|u_h|", extend="both")
+            figure.suptitle(f"Element: {self.case.element.capitalize()}", fontsize=11, color=INK)
+            figure.patch.set_facecolor(SURFACE)
+
+        self.plots["fields"] = figure
+
+    def report(self):
+        """Print one row per mesh refinement level, then the order each metric settled on."""
+        header = (f"{'# cells/axis':>{LABEL_WIDTH}} {'mesh step size h':>{SPACING_WIDTH}}"
+                  f" {'# of elements':>{ELEMENTS_WIDTH}}")
+        for metric in METRICS:
+            header += f" {metric.name:>{RATE_WIDTH}}"
+        print(header + f" {'# Newton iters':>{ITERATIONS_WIDTH}} {'r/r0':>{RESIDUAL_WIDTH}}"
+              f" Newton-Raphson status")
+
+        for index, level in enumerate(self.levels):
+            row = (f"{level.label:>{LABEL_WIDTH}} {level.spacing:>{SPACING_WIDTH}.5f}"
+                   f" {level.elements_count:>{ELEMENTS_WIDTH}}")
+            for metric in METRICS:
+                rate = level.rates[metric.name]
+                cell = f"{rate.value:.2f}" if rate.value is not None else (rate.reason or "")
+                accepted = index in self.results[metric.name].accepted_levels
+                row += f" {_paint(f'{cell:>{RATE_WIDTH}}', accepted)}"
+
+            newton = level.newton
+            if newton is None:
+                row += f" {'':>{ITERATIONS_WIDTH}} {'':>{RESIDUAL_WIDTH}} no newton"
+            else:
+                row += (f" {newton['iterations']:>{ITERATIONS_WIDTH}}"
+                        f" {newton['reduction']:>{RESIDUAL_WIDTH}.1e}"
+                        f" {_paint(newton['status'], newton['status'] == 'Converged')}")
+            print(row)
+
+        # Summary rows, as (label, one text per metric, pass per metric or None where uncoloured).
+        metric_results = [self.results[metric.name] for metric in METRICS]
+        rows = [
+            ("order p", [f"{r.order:.2f}" if r.order is not None else "not reached" for r in metric_results],
+             [r.passed for r in metric_results]),
+            ("expected", [f"{self.expected_order[metric.name]:.1f}" for metric in METRICS], None),
+        ]
+
+        print()
+        for label, texts, passed in rows:
+            row = f"{label:>{LABEL_WIDTH}} {'':>{SPACING_WIDTH}} {'':>{ELEMENTS_WIDTH}}"
+            for index, text in enumerate(texts):
+                padded = f"{text:>{RATE_WIDTH}}"
+                row += f" {padded if passed is None else _paint(padded, passed[index])}"
+            print(row)
+
+        if self.unconverged_levels:
+            print()
+            print(f"  Newton did not converge at: {', '.join(self.unconverged_levels)}")
+
+
+def overview(results):
+    """Print one line per case: the order each metric settled on, and whether its solves converged."""
+    print()
+    header = f"{'case':<44}"
+    for metric in METRICS:
+        header += f" {metric.name:>9}"
+    print(header + f" {'solver':>9}")
+
+    for name, study in results:
+        row = f"{name:<44}"
+        for metric in METRICS:
+            metric_result = None if study is None else study.results[metric.name]
+            has_order = metric_result is not None and metric_result.order is not None
+            cell = "error" if study is None else f"{metric_result.order:.2f}" if has_order else "--"
+            row += f" {_paint(f'{cell:>9}', study is not None and metric_result.passed)}"
+        solver = "error" if study is None else (
+            "ok" if not study.unconverged_levels else f"{len(study.unconverged_levels)} bad")
+        row += f" {_paint(f'{solver:>9}', solver == 'ok')}"
+        print(row)

From 480b379bbd7ca73c8280dfa10bf7874e5f45b842 Mon Sep 17 00:00:00 2001
From: Themis Skamagkis <themisskamagkis@gmail.com>
Date: Mon, 12 Oct 2026 09:49:19 +0200
Subject: [PATCH 30/55] [verification] Parse and validate cases

---
 verification/case.py | 209 +++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 209 insertions(+)
 create mode 100644 verification/case.py

diff --git a/verification/case.py b/verification/case.py
new file mode 100644
index 0000000..eb06ae9
--- /dev/null
+++ b/verification/case.py
@@ -0,0 +1,209 @@
+"""The case as one object: parsed, validated, and resolved into a config dictionary that a run needs."""
+
+import json
+import pathlib
+
+from common.geometry import Grid, Bar1D, Beam2D, Beam3D, Tire, Tire2D, Tire3D
+from common.sofa.conventions import ELEMENTS
+
+from . import manufactured
+from .manufactured import ManufacturedProblem
+from .manufactured.materials import LinearElastic
+from .study.metrics import METRICS
+
+CASES_ROOT = pathlib.Path(__file__).parent / "cases"
+
+# case "geometry.type" -> geometry class, keyed by the class's own name.
+GEOMETRIES = {cls.__name__: cls for cls in (Grid, Bar1D, Beam2D, Beam3D,
+                                            Tire, Tire2D, Tire3D)}
+
+# case "material.type" -> material class, keyed by the class's own name.
+MATERIALS = {cls.__name__: cls for cls in (LinearElastic,)}
+
+
+class _ManufacturedSolutionRegistry:
+    """(topological dim of the geometry, case "solution.type") -> ManufacturedSolution class.
+
+    Looked up by naming convention: a case's solution type "Quadratic" at dim 2 resolves to
+    manufactured.Quadratic2D. The dimension is the PDE's own, not the space it is solved in.
+    """
+
+    def __getitem__(self, key):
+        dim, type_name = key
+        class_name = f"{type_name}{dim}D"
+        try:
+            return getattr(manufactured, class_name)
+        except AttributeError:
+            raise KeyError(f"no manufactured solution {class_name!r} for dim={dim}, "
+                           f"solution type={type_name!r}")
+
+
+MANUFACTURED_SOLUTIONS = _ManufacturedSolutionRegistry()
+
+# Required keys for all cases. A missing key is an omission, an unknown key is a typo.
+REQUIRED = {"geometry", "element", "solution", "material", "forceField",
+            "errorQuadratureDegree", "sourceQuadratureDegree", "solvers", "mesh",
+            "asymptoticTolerance", "expectedOrder", "expectedOrderTolerance", "relativeNoiseFloor"}
+
+OPTIONAL = {
+        "tractionOn"      # {region: directions} handed to the manufactured traction; the rest stay Dirichlet.
+        , "gridMapping"   # Data for the mapping meshing the grid, where the element kind needs one.
+}
+
+MESH_KEYS = {"cells", "levels", "refinementRatio"}
+
+
+
+def _validate(name, config):
+    """Every required key present, nothing unknown, and one expectation per metric.
+
+    `solvers` and `forceField` pass their own keys through to SOFA Data verbatim, so this module
+    cannot know what's legal inside them and only checks their presence/shape, not their contents.
+    Every problem is collected before raising, so a case is corrected in one pass, not one key per run.
+    """
+    problems = []
+
+    # Top-level keys: exact match against REQUIRED.
+    missing = REQUIRED - set(config)
+    if missing:
+        problems.append(f"missing {sorted(missing)}")
+    unknown = set(config) - REQUIRED - OPTIONAL
+    if unknown:
+        problems.append(f"unknown {sorted(unknown)}")
+
+    # A grid-native element kind is meshed by no mapping, so its Data would go nowhere.
+    element = config.get("element")
+    if "gridMapping" in config and element in ELEMENTS and ELEMENTS[element].grid_mapping is None:
+        problems.append(f"gridMapping is stated but {element!r} is meshed straight from the grid")
+
+    # mesh: exact match against MESH_KEYS.
+    if "mesh" in config and set(config["mesh"]) != MESH_KEYS:
+        problems.append(f"mesh must state exactly {sorted(MESH_KEYS)}, got {sorted(config['mesh'])}")
+
+    # solution: just needs a type and an amplitude; the rest is that solution's own business.
+    solution_config = config.get("solution", {})
+    absent = {"type", "amplitude"} - set(solution_config)
+    if absent:
+        problems.append(f"solution must state {sorted(absent)}")
+
+    # material: must name a registered type.
+    material = config.get("material", {})
+    if "type" not in material:
+        problems.append(f"material must state a type, one of {sorted(MATERIALS)}")
+    elif material["type"] not in MATERIALS:
+        problems.append(f"unknown material {material['type']!r}, expected one of {sorted(MATERIALS)}")
+
+    # forceField: must be an object stating a type.
+    force_field = config.get("forceField")
+    if not isinstance(force_field, dict):
+        problems.append(f"forceField must be an object stating a type, got {force_field!r}")
+        force_field = {}
+    elif "type" not in force_field:
+        problems.append("forceField must state a type")
+    # forceField must not restate what the prefab already owns: name/template/topology follow
+    # from element and the mesh, so a case stating either could contradict what it already stated.
+    owned = {"name", "template", "topology"} & set(force_field)
+    if owned:
+        problems.append(f"forceField must not state {sorted(owned)}: the prefab sets them")
+    # forceField must not restate a constitutive parameter that belongs in material, so the
+    # manufactured source and the forceField under test can't drift apart.
+    clash = (set(force_field) & set(material)) - {"type"}
+    if clash:
+        problems.append(f"forceField and material both state {sorted(clash)}: state them in material")
+
+    # expectedOrder: one entry per registered metric, exactly.
+    names = {metric.name for metric in METRICS}
+    if "expectedOrder" in config and set(config["expectedOrder"]) != names:
+        problems.append(f"expectedOrder must state exactly {sorted(names)}, "
+                        f"got {sorted(config['expectedOrder'])}")
+
+    if problems:
+        raise ValueError(f"{name}: " + "; ".join(problems))
+
+
+def _validate_boundary_conditions(name, config, geometry):
+    """The case's traction regions exist, and state one direction per dimension of the space."""
+    traction_on = config.get("tractionOn", {})
+    problems = []
+
+    unknown = set(traction_on) - set(geometry.region_names)
+    if unknown:
+        problems.append(f"tractionOn names {sorted(unknown)}, not regions of "
+                        f"{type(geometry).__name__}: {sorted(geometry.region_names)}")
+
+    # One direction per dimension of the mesh; the clamp pads the ones it is embedded in.
+    for region, directions in traction_on.items():
+        if (len(directions) != geometry.dim
+                or any(direction not in (0, 1) for direction in directions)):
+            problems.append(f"tractionOn[{region!r}] must state {geometry.dim} directions, "
+                            f"each 0 or 1, got {directions}")
+
+    if problems:
+        raise ValueError(f"{name}: " + "; ".join(problems))
+
+
+class Case:
+    """One case: the study it specifies, with the geometry and the manufactured problem built."""
+
+    def __init__(self, path, config):
+        path = pathlib.Path(path)
+
+        # Extract root path and keep stem as the name of the case.
+        stem = path.resolve().with_suffix("")
+        try:
+            self.name = stem.relative_to(CASES_ROOT.resolve()).as_posix()
+        except ValueError:
+            self.name = f"{path.parent.name}/{path.stem}"
+        _validate(self.name, config)
+
+        # Geometry Class
+        geometry_config = dict(config["geometry"])
+        self.geometry = GEOMETRIES[geometry_config.pop("type")](**geometry_config)
+
+        _validate_boundary_conditions(self.name, config, self.geometry)
+
+        # Keyed on the geometry's topological dimension.
+        self.traction_on = config.get("tractionOn", {})
+        solution = MANUFACTURED_SOLUTIONS[(self.geometry.dim, config["solution"]["type"])](
+            config["solution"], self.geometry, self.traction_on)
+
+        # Nothing fixed anywhere leaves the rigid modes in the system.
+        if not solution.prescribe_displacement_on:
+            raise ValueError(f"{self.name}: tractionOn frees every direction of every region, "
+                             f"which leaves the problem singular")
+
+        # ManufacturedProblem Class
+        material = MATERIALS[config["material"]["type"]]
+        self.manufactured_problem = ManufacturedProblem(
+            solution, material(*material.get_material_parameters(
+                config["material"], self.geometry.spatial_dimensions)),
+            self.geometry.spatial_dimensions)
+
+        self.element = config["element"]
+        self.material = config["material"]
+        self.force_field = config["forceField"]
+        self.solvers = config["solvers"]
+        self.mesh = config["mesh"]
+        self.error_quadrature_degree = config["errorQuadratureDegree"]
+        self.source_quadrature_degree = config["sourceQuadratureDegree"]
+        self.asymptotic_tolerance = config["asymptoticTolerance"]
+        self.expected_order = config["expectedOrder"]
+        self.expected_order_tolerance = config["expectedOrderTolerance"]
+        self.relative_noise_floor = config["relativeNoiseFloor"]
+        self.grid_mapping = config.get("gridMapping", {})
+
+    @classmethod
+    def load(cls, path):
+        with open(path) as f:
+            return cls(path, json.load(f))
+
+    @property
+    def equation(self):
+        """The manufactured solution in math form, for figures and reports."""
+        return self.manufactured_problem.equation
+
+    def levels(self):
+        """Cells per axis at each refinement level, coarsest first."""
+        ratio = self.mesh["refinementRatio"]
+        return [[count * ratio ** level for count in self.mesh["cells"]]
+                for level in range(self.mesh["levels"])]

From 86ca36e839b46619f05dface3be1a49e41eae4f1 Mon Sep 17 00:00:00 2001
From: Themis Skamagkis <themisskamagkis@gmail.com>
Date: Mon, 12 Oct 2026 09:49:19 +0200
Subject: [PATCH 31/55] [verification] Add the verification runner

---
 verification/run.py | 133 ++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 133 insertions(+)
 create mode 100644 verification/run.py

diff --git a/verification/run.py b/verification/run.py
new file mode 100644
index 0000000..58e04a8
--- /dev/null
+++ b/verification/run.py
@@ -0,0 +1,133 @@
+"""Runner for the verification suite: loads a case and runs its verification study."""
+
+import argparse
+import os
+import pathlib
+import sys
+import traceback
+
+# Make the plugin's packages importable when this file is run directly.
+sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
+
+from verification.case import Case
+from verification.scene.scene import MMSScene
+from verification.study.study import ErrorConvergenceStudy, clear_results, overview
+
+CASES_ROOT = pathlib.Path(__file__).parent / "cases"
+
+
+def resolve_path(path):
+    """`path` as given when it is absolute or already exists, else taken from the cases folder."""
+    path = pathlib.Path(path)
+    return path if path.is_absolute() or path.exists() else CASES_ROOT / path
+
+
+def run_all(directory, args):
+    """Run every case under `directory`, each with its own table, then one line per case.
+
+    Returns [(name, study)], the study being None where that case raised.
+    """
+    TESTED_COMPONENTS = ("LinearSmallStrainFEMForceField",)
+
+    results = []
+    for component in TESTED_COMPONENTS:
+        for path in sorted((directory / component).rglob("*.json")):
+            name = path.relative_to(directory).with_suffix("").as_posix()
+            print(f"\n--- {name} ---")
+            try:
+                results.append((name, run(path)))
+            except Exception as error:
+                # One broken case should not stop the others.
+                print(f"  failed: {type(error).__name__}: {error}")
+                # Keep the traceback right below the case that raised it.
+                if args.traceback:
+                    traceback.print_exc(file=sys.stdout)
+                results.append((name, None))
+    overview(results)
+    return results
+
+
+def run(case_path):
+    """Load one case, run its verification study, and write what it produced to disk."""
+    case = Case.load(case_path)
+    study = ErrorConvergenceStudy(case)
+    study.run()
+    study.write_results()
+    return study
+
+
+def createScene(root):
+    """The runSofa entry point, allowing to inspect a single mesh level of a case in the GUI.
+
+        runSofa -l SofaPython3 run.py --argv <case> --argv <level>
+
+    `<level>` indexes the case's refinement sequence, so `0` is its coarsest mesh.
+    """
+    if len(sys.argv) != 3:
+        sys.exit("run.py under runSofa expects: --argv <case> --argv <level>")
+    case = Case.load(resolve_path(sys.argv[1]))
+    cells = case.levels()[int(sys.argv[2])]
+
+    # Add VisualStyle
+    root.addObject('RequiredPlugin', name='visual', pluginName=["Sofa.Component.Visual"])
+    root.addObject('VisualStyle', displayFlags='showForceFields showBehaviorModels')
+
+    # Build the scene exercising the given mesh level with the MMS
+    MMSScene(case, cells).build(root)
+
+
+def parse_arguments():
+    """Parse the command line into the whole run configuration.
+
+    Returns
+    -------
+    argparse.Namespace
+        case : str or None
+            Path to the single case to run. ``None`` when ``--all`` was given instead. Taken
+            relative to the cases folder by ``resolve_path`` when it does not resolve as typed.
+        all : bool
+            Run every case under the cases folder's force-field subdirectories. Mutually
+            exclusive with ``case``.
+        clean : bool
+            Delete what previous runs wrote before running.
+        traceback : bool
+            Whether a case that raises under ``--all`` also prints its full stack trace. False by
+            default.
+    """
+    parser = argparse.ArgumentParser(description=__doc__)
+    # Must choose either --all or pass the case path
+    which_cases = parser.add_mutually_exclusive_group()
+    which_cases.add_argument("case", nargs="?", help="path to a case")
+    which_cases.add_argument("--all", action="store_true", help="every case under the cases folder")
+    parser.add_argument("--clean", action="store_true",
+                        help="delete the results of previous runs before running")
+    parser.add_argument("--traceback", action="store_true",
+                        help="full stack for a case that raises under --all")
+
+    arguments = parser.parse_args()
+    if not (arguments.case or arguments.all or arguments.clean):
+        parser.error("one of case, --all or --clean is required")
+    return arguments
+
+
+if __name__ == "__main__":
+    args = parse_arguments()
+
+    if args.clean:
+        clear_results()
+
+    # Run all the cases in verification
+    if args.all:
+        failed = [name for name, study in run_all(CASES_ROOT, args)
+                  if study is None or not study.verified]
+    # Run the prescribed case
+    elif args.case:
+        study = run(resolve_path(args.case))
+        failed = [] if study.verified else [study.name]
+    # Nothing left to do: --clean was given on its own
+    else:
+        failed = []
+
+    # On failure, exit with error and report the problematic cases
+    if failed:
+        sys.exit(f"\nfailed: {', '.join(failed)}")

From f9cdff8b508d2b56b5fdf91580966f15aa52f7e8 Mon Sep 17 00:00:00 2001
From: Themis Skamagkis <themisskamagkis@gmail.com>
Date: Mon, 12 Oct 2026 09:49:19 +0200
Subject: [PATCH 32/55] [docs] Add the introductory MMS notebook

---
 verification/docs/IntroductionMMS.ipynb | 202 ++++++++++++++++++++++++
 1 file changed, 202 insertions(+)
 create mode 100644 verification/docs/IntroductionMMS.ipynb

diff --git a/verification/docs/IntroductionMMS.ipynb b/verification/docs/IntroductionMMS.ipynb
new file mode 100644
index 0000000..e97ba2a
--- /dev/null
+++ b/verification/docs/IntroductionMMS.ipynb
@@ -0,0 +1,202 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "01-md-title",
+   "metadata": {},
+   "source": [
+    "# Code Verification by the Method of Manufactured Solutions\n",
+    "\n",
+    "This notebook presents the Method of Manufactured Solutions (MMS) used for code verification and provides a walkthrough of the code that implements it.\n",
+    "\n",
+    "Verification and Validation answer two distinct questions [[1]](#ref1)[[7]](#ref7): \n",
+    "\n",
+    "- **Verification** Does the code solve its equations correctly, at the order of accuracy the discretisation claims?\n",
+    "- **Validation**  Are these the right equations for the physical problem?\n",
+    "\n",
+    "In short, verification is *solving the equations right*, and validation is *solving the right equations* [[4]](#ref4).\n",
+    "\n",
+    "Verifying a code means comparing its solution against an exact one. The problem with classical exact solutions is that they exist\n",
+    "only for very simple geometries, loads and material laws, and even then they leave terms of the\n",
+    "equations untested [[1]](#ref1)[[3]](#ref3). A code can pass them and still be wrong. \n",
+    "MMS removes this limitation: it builds the problem around a solution we construct, so an exact\n",
+    "solution is always available."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "02-md-mms-in-a-nutshell",
+   "metadata": {},
+   "source": [
+    "## The MMS in a nutshell\n",
+    "\n",
+    "Suppose we wish to solve the mathematical model $\\mathcal{L}(u) = 0$, where $\\mathcal{L}$ is the\n",
+    "differential operator and $u$ the unknown. The main idea is to choose a manufactured solution\n",
+    "$\\hat{u}$ first, then apply the governing equations to it: what remains, $s = \\mathcal{L}(\\hat{u})$,\n",
+    "is the source term. By construction, $\\hat{u}$ is the exact solution of the modified model\n",
+    "[[5]](#ref5)\n",
+    "\n",
+    "$$\\mathcal{L}(u) = s$$\n",
+    "\n",
+    "The boundary conditions are derived from $\\hat{u}$ as well. Three properties make this practical\n",
+    "[[5]](#ref5):\n",
+    "\n",
+    "- $s$ depends on the coordinates only, not on $u$. It enters the code as a known load, that is, as\n",
+    "  input.\n",
+    "- Since $s$ is a known function, it needs no discrete version of $\\mathcal{L}$.\n",
+    "- The left-hand side $\\mathcal{L}(u)$ is unchanged, so the discretisation under test is not\n",
+    "  modified; only its input is.\n",
+    "\n",
+    "In short, the expectation is that solving the modified problem gives back $\\hat{u}$.\n",
+    "\n",
+    "In practice, the code does not solve $\\mathcal{L}(u) = s$ itself but a discrete version of it: a\n",
+    "finite element code solves its weak form on a mesh. Its solution $u_h$ is then only an\n",
+    "approximation of $\\hat{u}$, and the difference $e = u_h - \\hat{u}$ is the error. Because $\\hat{u}$\n",
+    "is known exactly, this error can be computed rather than estimated.\n",
+    "\n",
+    "### Measuring the error\n",
+    "\n",
+    "The error has several sources [[2]](#ref2):\n",
+    "\n",
+    "$$e = \\varepsilon_h + \\varepsilon_\\text{ro} + \\varepsilon_\\text{it} + \\varepsilon_\\text{bug}$$\n",
+    "\n",
+    "- $\\varepsilon_h$, the *discretisation error*, which comes from the mesh and vanishes as it is\n",
+    "  refined;\n",
+    "- $\\varepsilon_\\text{ro}$, the *round-off error*, from finite-precision arithmetic;\n",
+    "- $\\varepsilon_\\text{it}$, the *iterative error*, from stopping the linear or Newton solver at a\n",
+    "  finite tolerance;\n",
+    "- $\\varepsilon_\\text{bug}$, the error caused by *coding mistakes*.\n",
+    "\n",
+    "Under a correct implementation $\\varepsilon_\\text{bug} = 0$, and only $\\varepsilon_h$ should remain,\n",
+    "provided $\\varepsilon_\\text{ro}$ and $\\varepsilon_\\text{it}$ are kept far below it.\n",
+    "\n",
+    "The error is a field over $\\Omega$, so its size is measured with norms integrated over the domain:\n",
+    "\n",
+    "$$\\|e\\|_{L^2} = \\left(\\int_\\Omega |e|^2\\,d\\Omega\\right)^{1/2}, \\qquad\n",
+    "  |e|_{H^1} = \\left(\\int_\\Omega \\|\\nabla e\\|^2\\,d\\Omega\\right)^{1/2}$$\n",
+    "\n",
+    "The $L^2$ norm measures the error itself, the $H^1$ semi-norm its gradient.\n",
+    "\n",
+    "Order verification assumes that the discretisation error can be expanded in powers of the mesh size\n",
+    "$h$ [[5]](#ref5):\n",
+    "\n",
+    "$$\\|\\varepsilon_h\\| = C\\,h^p + O(h^{p+1})$$\n",
+    "\n",
+    "$p$ is the formal order of accuracy of the discretisation. $C$ does not depend on $h$, but on the\n",
+    "higher derivatives of the manufactured solution [[5]](#ref5). As $h$ decreases, the higher-order\n",
+    "terms vanish faster than the first one, so on fine enough meshes the leading term dominates and\n",
+    "$\\|\\varepsilon_h\\| \\approx C\\,h^p$. For elements whose shape functions are polynomials of degree\n",
+    "$q$, $p = q + 1$ in the $L^2$ norm and $p = q$ in the $H^1$ semi-norm [[8]](#ref8).\n",
+    "\n",
+    "$\\varepsilon_\\text{ro}$ and $\\varepsilon_\\text{it}$ do not decrease with $h$. They set a noise\n",
+    "floor, and an error below it must not be used to measure the order: the ratio of two such errors\n",
+    "reflects round-off and solver tolerance, not the discretisation [[2]](#ref2). A coding mistake\n",
+    "shows up as an error that does not vanish, or that vanishes at the wrong rate.\n",
+    "\n",
+    "The code is judged by the rate at which the error decreases, not by its size. A single error value\n",
+    "on a single mesh does not tell whether the code is correct: how small is small enough depends on the\n",
+    "solution, the units and the mesh, all of which enter $C$, while $p$ depends on the discretisation\n",
+    "alone. The mesh is therefore refined by a factor $r = h_\\text{coarse}/h_\\text{fine} > 1$, and the\n",
+    "ratio of the two errors gives the observed order of accuracy [[5]](#ref5):\n",
+    "\n",
+    "$$\\hat{p} = \\frac{\\ln\\!\\left(\\|e_{rh}\\| / \\|e_{h}\\|\\right)}{\\ln r}$$\n",
+    "\n",
+    "The test passes when $\\hat{p}$ matches $p$. Checking the order is a much stronger criterion than\n",
+    "checking that the error is small, or even that it decreases [[3]](#ref3).\n",
+    "\n",
+    "Everything under test is implemented in SOFA, in C++: the force fields\n",
+    "(`LinearSmallStrainFEMForceField`), the\n",
+    "solvers, and the components that turn a body force into nodal loads (`VectorSourceTerm`,\n",
+    "`FEMSourceTermIntegrator`). The suite does not modify any of them; they are the discrete\n",
+    "$\\mathcal{L}$ whose left-hand side stays unchanged.\n",
+    "\n",
+    "The Python scripts do everything around them:\n",
+    "\n",
+    "- derive the source term $s$ and the boundary conditions from $\\hat{u}$ and the material law;\n",
+    "- apply them to a SOFA scene as loads and constraints;\n",
+    "- solve the scene on a sequence of refined meshes, measure the error on each and read its order.\n",
+    "\n",
+    "### Further reading\n",
+    "\n",
+    "- Roache [[1]](#ref1): a short review of MMS for code verification by one of its originators.\n",
+    "- Roy [[2]](#ref2): a review of code and solution verification, with a six-step procedure\n",
+    "  and guidance on choosing manufactured solutions.\n",
+    "- Salari and Knupp [[3]](#ref3): the most complete treatment, including the blind test in which\n",
+    "  coding mistakes were planted in a code to measure what MMS detects.\n",
+    "- Oberkampf and Roy [[5]](#ref5), Chapters 5 and 6: the textbook reference on code verification\n",
+    "  and MMS.\n",
+    "- Aycock *et al.* [[6]](#ref6): MMS applied to elastostatic solid mechanics in a finite element\n",
+    "  solver, the closest precedent for this suite."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "03-md-references",
+   "metadata": {},
+   "source": [
+    "<a id=\"references\"></a>\n",
+    "## References\n",
+    "\n",
+    "<a id=\"ref1\"></a>\n",
+    "**[1]** P. J. Roache, \"Code verification by the method of manufactured solutions,\"\n",
+    "*Journal of Fluids Engineering*, 124(1):4&ndash;10, 2002.\n",
+    "[doi:10.1115/1.1436090](https://doi.org/10.1115/1.1436090)\n",
+    "\n",
+    "<a id=\"ref2\"></a>\n",
+    "**[2]** C. J. Roy, \"Review of code and solution verification procedures for computational\n",
+    "simulation,\" *Journal of Computational Physics*, 205(1):131&ndash;156, 2005.\n",
+    "[doi:10.1016/j.jcp.2004.10.036](https://doi.org/10.1016/j.jcp.2004.10.036)\n",
+    "\n",
+    "<a id=\"ref3\"></a>\n",
+    "**[3]** K. Salari and P. Knupp, *Code Verification by the Method of Manufactured Solutions*,\n",
+    "SAND2000-1444, Sandia National Laboratories, Albuquerque, NM, June 2000.\n",
+    "\n",
+    "<a id=\"ref4\"></a>\n",
+    "**[4]** S. Étienne, A. Garon and D. Pelletier, \"Some manufactured solutions for verification of\n",
+    "fluid-structure interaction codes,\" *Computers and Structures*, 106&ndash;107:56&ndash;67, 2012.\n",
+    "[doi:10.1016/j.compstruc.2012.04.006](https://doi.org/10.1016/j.compstruc.2012.04.006)\n",
+    "\n",
+    "<a id=\"ref5\"></a>\n",
+    "**[5]** W. L. Oberkampf and C. J. Roy, *Verification and Validation in Scientific Computing*,\n",
+    "Cambridge University Press, 2010. ISBN 978-0-521-11360-1. &mdash; *Chapter 5 (code verification and\n",
+    "order verification) and Chapter 6, §6.3 (MMS).*\n",
+    "\n",
+    "<a id=\"ref6\"></a>\n",
+    "**[6]** K. I. Aycock, N. Rebelo and B. A. Craven, \"Method of manufactured solutions code\n",
+    "verification of elastostatic solid mechanics problems in a commercial finite element solver,\"\n",
+    "*Computers and Structures*, 229:106175, 2020.\n",
+    "[doi:10.1016/j.compstruc.2019.106175](https://doi.org/10.1016/j.compstruc.2019.106175)\n",
+    "\n",
+    "<a id=\"ref7\"></a>\n",
+    "**[7]** W. L. Oberkampf and T. G. Trucano, \"Verification and validation in computational fluid\n",
+    "dynamics,\" *Progress in Aerospace Sciences*, 38(3):209&ndash;272, 2002.\n",
+    "[doi:10.1016/S0376-0421(02)00005-2](https://doi.org/10.1016/S0376-0421%2802%2900005-2)\n",
+    "\n",
+    "<a id=\"ref8\"></a>\n",
+    "**[8]** P. G. Ciarlet, *The Finite Element Method for Elliptic Problems*, North-Holland, Amsterdam,\n",
+    "1978. Reprinted by SIAM, Classics in Applied Mathematics 40, 2002."
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.12.12"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

From e2cbfcc85a021d2d5ab82393cc4a2781486b28dd Mon Sep 17 00:00:00 2001
From: Themis Skamagkis <themisskamagkis@gmail.com>
Date: Mon, 12 Oct 2026 10:13:06 +0200
Subject: [PATCH 33/55] [docs] Update README

---
 README.md | 15 +++++++++++----
 1 file changed, 11 insertions(+), 4 deletions(-)

diff --git a/README.md b/README.md
index 843a7da..17d619f 100644
--- a/README.md
+++ b/README.md
@@ -1,9 +1,6 @@
 # Verification-Validation
 Initially funded by TIRREX, this repository proposes a framework for verification and validation. It relies on both manufactured solutions and comparison with alternative simulation technologies. This work therefore proves the reliability and compliance of the SOFA models &amp; algorithms.
 
-For how the repository is organised, how the suite works and how to extend it, see
-[docs/architecture.md](docs/architecture.md).
-
 ## Verification
 
 **Verification** uses the method of manufactured solutions: an exact displacement is chosen, the
@@ -11,6 +8,8 @@ loads that produce it are derived, and SOFA solves the problem on a sequence of
 The error must shrink at the rate the discretisation predicts. Covered force fields:
 - `LinearSmallStrainFEMForceField`
 
+To get started, read [`verification/docs/IntroductionMMS.ipynb`](verification/docs/IntroductionMMS.ipynb): it presents the method and builds one verification study step by step, from the manufactured solution to the input file the suite runs.
+
 ## Layout
 
 ```
@@ -28,8 +27,16 @@ verification/                            the verification suite
 
 ```bash
 python verification/run.py <path-to-json-file>  # one test case
-python verification/run.py --all                # every test case
+python verification/run.py --all                # every test case, then one summary line per case
 python verification/run.py --clean              # delete previous results
 ```
 
 Results are written to `verification/results/`.
+
+To inspect one mesh of a case in the SOFA GUI, the runner also acts as a SOFA scene file:
+
+```bash
+runSofa -l SofaPython3 verification/run.py --argv <path-to-json-file> --argv 0
+```
+
+The first argument is the path to the .json case file and the second is the mesh refinement level to be built.

From e41fc5d06092111fa4327a10bce216c011ceaf1b Mon Sep 17 00:00:00 2001
From: Themis Skamagkis <themisskamagkis@gmail.com>
Date: Mon, 12 Oct 2026 10:41:11 +0200
Subject: [PATCH 34/55] [docs] Add a concrete trigonometric function example

---
 verification/docs/IntroductionMMS.ipynb | 767 ++++++++++++++++++++++--
 1 file changed, 727 insertions(+), 40 deletions(-)

diff --git a/verification/docs/IntroductionMMS.ipynb b/verification/docs/IntroductionMMS.ipynb
index e97ba2a..8d5c668 100644
--- a/verification/docs/IntroductionMMS.ipynb
+++ b/verification/docs/IntroductionMMS.ipynb
@@ -9,16 +9,16 @@
     "\n",
     "This notebook presents the Method of Manufactured Solutions (MMS) used for code verification and provides a walkthrough of the code that implements it.\n",
     "\n",
-    "Verification and Validation answer two distinct questions [[1]](#ref1)[[7]](#ref7): \n",
+    "Verification and Validation answer two distinct questions [[2]](#ref2)[[6]](#ref6): \n",
     "\n",
     "- **Verification** Does the code solve its equations correctly, at the order of accuracy the discretisation claims?\n",
     "- **Validation**  Are these the right equations for the physical problem?\n",
     "\n",
-    "In short, verification is *solving the equations right*, and validation is *solving the right equations* [[4]](#ref4).\n",
+    "In short, verification is *solving the equations right*, and validation is *solving the right equations* [[7]](#ref7).\n",
     "\n",
     "Verifying a code means comparing its solution against an exact one. The problem with classical exact solutions is that they exist\n",
     "only for very simple geometries, loads and material laws, and even then they leave terms of the\n",
-    "equations untested [[1]](#ref1)[[3]](#ref3). A code can pass them and still be wrong. \n",
+    "equations untested [[2]](#ref2)[[4]](#ref4). A code can pass them and still be wrong. \n",
     "MMS removes this limitation: it builds the problem around a solution we construct, so an exact\n",
     "solution is always available."
    ]
@@ -34,7 +34,7 @@
     "differential operator and $u$ the unknown. The main idea is to choose a manufactured solution\n",
     "$\\hat{u}$ first, then apply the governing equations to it: what remains, $s = \\mathcal{L}(\\hat{u})$,\n",
     "is the source term. By construction, $\\hat{u}$ is the exact solution of the modified model\n",
-    "[[5]](#ref5)\n",
+    "[[1]](#ref1)[[2]](#ref2)[[5]](#ref5)\n",
     "\n",
     "$$\\mathcal{L}(u) = s$$\n",
     "\n",
@@ -43,7 +43,8 @@
     "\n",
     "- $s$ depends on the coordinates only, not on $u$. It enters the code as a known load, that is, as\n",
     "  input.\n",
-    "- Since $s$ is a known function, it needs no discrete version of $\\mathcal{L}$.\n",
+    "- Since $s$ is a known function, no discrete version of $\\mathcal{L}$ is applied to it. It is\n",
+    "  only integrated into nodal loads, by quadrature.\n",
     "- The left-hand side $\\mathcal{L}(u)$ is unchanged, so the discretisation under test is not\n",
     "  modified; only its input is.\n",
     "\n",
@@ -56,7 +57,9 @@
     "\n",
     "### Measuring the error\n",
     "\n",
-    "The error has several sources [[2]](#ref2):\n",
+    "Numerical error has three sources [[3]](#ref3): the discretisation, finite-precision arithmetic\n",
+    "and the tolerance the solvers stop at. A coding mistake adds a fourth contribution, the one the\n",
+    "study is looking for:\n",
     "\n",
     "$$e = \\varepsilon_h + \\varepsilon_\\text{ro} + \\varepsilon_\\text{it} + \\varepsilon_\\text{bug}$$\n",
     "\n",
@@ -76,6 +79,7 @@
     "  |e|_{H^1} = \\left(\\int_\\Omega \\|\\nabla e\\|^2\\,d\\Omega\\right)^{1/2}$$\n",
     "\n",
     "The $L^2$ norm measures the error itself, the $H^1$ semi-norm its gradient.\n",
+    "This suite also compares the strain energy SOFA computes for $u_h$ with the exact one, $E(\\hat{u})$, through $\\sqrt{2\\,|E(\\hat{u}) - E(u_h)|}$.\n",
     "\n",
     "Order verification assumes that the discretisation error can be expanded in powers of the mesh size\n",
     "$h$ [[5]](#ref5):\n",
@@ -85,13 +89,16 @@
     "$p$ is the formal order of accuracy of the discretisation. $C$ does not depend on $h$, but on the\n",
     "higher derivatives of the manufactured solution [[5]](#ref5). As $h$ decreases, the higher-order\n",
     "terms vanish faster than the first one, so on fine enough meshes the leading term dominates and\n",
-    "$\\|\\varepsilon_h\\| \\approx C\\,h^p$. For elements whose shape functions are polynomials of degree\n",
+    "$\\|\\varepsilon_h\\| \\approx C\\,h^p$. Such meshes are in the *asymptotic regime*. On coarser ones, in the pre-asymptotic regime, the higher-order terms still weigh, and the observed order can differ from $p$ even for a correct code. For elements whose shape functions are polynomials of degree\n",
     "$q$, $p = q + 1$ in the $L^2$ norm and $p = q$ in the $H^1$ semi-norm [[8]](#ref8).\n",
     "\n",
-    "$\\varepsilon_\\text{ro}$ and $\\varepsilon_\\text{it}$ do not decrease with $h$. They set a noise\n",
-    "floor, and an error below it must not be used to measure the order: the ratio of two such errors\n",
-    "reflects round-off and solver tolerance, not the discretisation [[2]](#ref2). A coding mistake\n",
-    "shows up as an error that does not vanish, or that vanishes at the wrong rate.\n",
+    "$\\varepsilon_\\text{ro}$ and $\\varepsilon_\\text{it}$ do not follow $h$: the iterative error is set\n",
+    "by the solver tolerance, and round-off grows as the mesh is refined and the system becomes larger\n",
+    "and worse conditioned. They set a noise floor, and an error below it must not be used to measure\n",
+    "the order: the ratio of two such errors reflects round-off and solver tolerance, not the\n",
+    "discretisation. Both must stay well below $\\varepsilon_h$ on the finest mesh [[5]](#ref5). A coding\n",
+    "mistake shows up as an error that does not vanish, or that vanishes at the wrong rate\n",
+    "[[4]](#ref4).\n",
     "\n",
     "The code is judged by the rate at which the error decreases, not by its size. A single error value\n",
     "on a single mesh does not tell whether the code is correct: how small is small enough depends on the\n",
@@ -102,7 +109,7 @@
     "$$\\hat{p} = \\frac{\\ln\\!\\left(\\|e_{rh}\\| / \\|e_{h}\\|\\right)}{\\ln r}$$\n",
     "\n",
     "The test passes when $\\hat{p}$ matches $p$. Checking the order is a much stronger criterion than\n",
-    "checking that the error is small, or even that it decreases [[3]](#ref3).\n",
+    "checking that the error is small, or even that it decreases [[4]](#ref4).\n",
     "\n",
     "Everything under test is implemented in SOFA, in C++: the force fields\n",
     "(`LinearSmallStrainFEMForceField`), the\n",
@@ -118,43 +125,715 @@
     "\n",
     "### Further reading\n",
     "\n",
-    "- Roache [[1]](#ref1): a short review of MMS for code verification by one of its originators.\n",
-    "- Roy [[2]](#ref2): a review of code and solution verification, with a six-step procedure\n",
-    "  and guidance on choosing manufactured solutions.\n",
-    "- Salari and Knupp [[3]](#ref3): the most complete treatment, including the blind test in which\n",
-    "  coding mistakes were planted in a code to measure what MMS detects.\n",
-    "- Oberkampf and Roy [[5]](#ref5), Chapters 5 and 6: the textbook reference on code verification\n",
-    "  and MMS.\n",
-    "- Aycock *et al.* [[6]](#ref6): MMS applied to elastostatic solid mechanics in a finite element\n",
-    "  solver, the closest precedent for this suite."
+    "- Steinberg and Roache [[1]](#ref1): the original presentation of the method.\n",
+    "- Roache [[2]](#ref2): a short review of MMS for code verification.\n",
+    "- Roy [[3]](#ref3): a review of code and solution verification, with guidance on choosing manufactured solutions.\n",
+    "- Salari and Knupp [[4]](#ref4): the most complete treatment, including examples of coding mistakes detected by MSS\n",
+    "- Oberkampf and Roy [[5]](#ref5), Chapters 5 and 6: the textbook reference on code verification and MMS.\n",
+    "- Aycock *et al.* [[9]](#ref9): MMS applied to elastostatic solid mechanics in a finite element solver.\n"
    ]
   },
   {
    "cell_type": "markdown",
-   "id": "03-md-references",
+   "id": "03-md-example-problem",
+   "metadata": {},
+   "source": [
+    "## Example: a trigonometric manufactured solution on a 2D square\n",
+    "\n",
+    "The rest of this notebook walks through one verification study, building each piece of it in turn.\n",
+    "\n",
+    "### The governing equations\n",
+    "\n",
+    "Under small-strain deformations, the displacement $u$ of a linear elastic body in static equilibrium\n",
+    "and without external loads satisfies\n",
+    "\n",
+    "$$\\mathcal{L}(u) = -\\nabla\\cdot\\sigma(u) = 0, \\qquad\n",
+    "  \\sigma = \\lambda\\,\\mathrm{tr}(\\varepsilon)\\,I + 2\\mu\\,\\varepsilon, \\qquad\n",
+    "  \\varepsilon = \\tfrac{1}{2}\\left(\\nabla u + \\nabla u^T\\right)$$\n",
+    "\n",
+    "This is the equation solved by the `LinearSmallStrainFEMForceField` component in SOFA. $\\mu$ and\n",
+    "$\\lambda$ are the Lamé parameters, obtained from Young's modulus $E$ and Poisson's ratio $\\nu$.\n",
+    "\n",
+    "The body $\\Omega$ is a rectangle of length $L$ and width $W$, here a square with $L = W = 1$. Its\n",
+    "boundary is made of four sides: `left` ($x = 0$), `right` ($x = L$), `bottom` ($y = 0$) and `top`\n",
+    "($y = W$).\n",
+    "\n",
+    "### The manufactured solution\n",
+    "\n",
+    "We choose one sine period along each axis:\n",
+    "\n",
+    "$$\\hat{u}(x, y) = A \\begin{bmatrix} \\sin(kx)\\cos(ky) \\\\ \\cos(kx)\\sin(ky) \\end{bmatrix},\n",
+    "  \\qquad k = \\frac{2\\pi}{L} = \\frac{2\\pi}{W}$$\n",
+    "\n",
+    "$\\hat{u}$ is not a physical solution and does not need to be one. It only has to be smooth, to leave\n",
+    "no derivative appearing in the equations identically zero, and to have derivatives small enough that\n",
+    "the asymptotic regime is reached on affordable meshes [[4]](#ref4)[[5]](#ref5). The amplitude $A$\n",
+    "sets the last of these.\n",
+    "\n",
+    "Its gradient\n",
+    "\n",
+    "$$\\hat{\\varepsilon} = \\nabla\\hat{u} = A k \\begin{bmatrix} c & -s \\\\ -s & c \\end{bmatrix},\n",
+    "  \\qquad c = \\cos(kx)\\cos(ky), \\quad s = \\sin(kx)\\sin(ky)$$\n",
+    "\n",
+    "and the stress\n",
+    "\n",
+    "$$\\hat{\\sigma} = 2 A k \\begin{bmatrix} (\\lambda+\\mu)\\,c & -\\mu\\,s \\\\ -\\mu\\,s & (\\lambda+\\mu)\\,c \\end{bmatrix}$$\n",
+    "\n",
+    "Applying $\\mathcal{L}$ to $\\hat{u}$ gives the source term:\n",
+    "\n",
+    "$$s = \\mathcal{L}(\\hat{u}) = -\\nabla\\cdot\\hat{\\sigma} = 2k^2(\\lambda + 2\\mu)\\,\\hat{u}$$\n",
+    "\n",
+    "\n",
+    "### The boundary conditions\n",
+    "\n",
+    "The boundary conditions are taken from $\\hat{u}$ as well. Multiplying $\\mathcal{L}(u) = s$ by a test\n",
+    "function $v$, integrating over $\\Omega$ and integrating by parts gives the weak form\n",
+    "\n",
+    "$$\\int_\\Omega \\sigma : \\nabla v \\, d\\Omega - \\int_{\\partial\\Omega} t \\cdot v \\, d\\Gamma\n",
+    "  = \\int_\\Omega s \\cdot v \\, d\\Omega, \\qquad t = \\sigma\\,n$$\n",
+    "\n",
+    "where $n$ is the outward normal and $t$ the traction. On each side, and in each direction, a\n",
+    "boundary condition prescribes one of two quantities:\n",
+    "\n",
+    "- the displacement, $u = \\hat{u}$ (Dirichlet). The test function vanishes there, and so does the\n",
+    "  boundary term;\n",
+    "- the traction, $t = \\hat{t} = \\hat{\\sigma}\\,n$ (Neumann), which enters the boundary term as a known\n",
+    "  load.\n",
+    "\n",
+    "Since $\\hat{u}$ is known everywhere, both conditions are consistent with it, and which one goes where\n",
+    "is a choice. Here, the displacement is prescribed on the whole boundary. It is not a zero Dirichlet condition, for instance on `left` it is $\\hat{u}(0, y) = (0,\\ A\\sin(ky))$.\n",
+    "\n",
+    "The modified problem is then\n",
+    "\n",
+    "$$-\\nabla\\cdot\\sigma(u) = 2k^2(\\lambda + 2\\mu)\\,\\hat{u} \\quad \\text{in } \\Omega, \\qquad\n",
+    "  u = \\hat{u} \\quad \\text{on } \\partial\\Omega$$\n",
+    "\n",
+    "and $\\hat{u}$ is its exact solution. In this example $A = 0.02$, $E = 10^6$ and $\\nu = 0.3$."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "04-md-problem-to-code",
+   "metadata": {},
+   "source": [
+    "### From the problem to the code\n",
+    "\n",
+    "Setting up an MMS study is systematic, so the problem above splits into separate concepts, each\n",
+    "with its own place in the repository:\n",
+    "\n",
+    "```\n",
+    "Verification-Validation/\n",
+    "├── common/                   code shared by the suites\n",
+    "│   ├── geometry/             the domains Ω and their named boundary regions\n",
+    "│   └── sofa/\n",
+    "│       ├── prefabs/          how SOFA builds an elastic body on each geometry\n",
+    "│       ├── controllers.py    SOFA controllers the scenes need, e.g. the prescribed displacement clamp\n",
+    "│       └── scene.py          the base of every SOFA scene\n",
+    "└── verification/\n",
+    "    ├── manufactured/         modules that implement manufactured solutions and material laws; combining them forms the manufactured problem to be solved\n",
+    "    ├── scene/                the SOFA scene that loads a body with a manufactured problem\n",
+    "    ├── study/                the error measurement and the convergence study\n",
+    "    ├── case.py               reads a case: the study's configuration\n",
+    "    └── run.py                runs a case\n",
+    "```\n",
+    "\n",
+    "Each of these is introduced below, where the walkthrough reaches it."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "05-md-input-file",
+   "metadata": {},
+   "source": [
+    "### The input file\n",
+    "\n",
+    "A study is set up by an input file, a JSON file the suite reads into a dictionary. The walkthrough\n",
+    "builds it up as it goes, one key at a time. The problem above needs three keys:\n",
+    "\n",
+    "- `geometry` names the geometry and gives its size;\n",
+    "- `solution` names the manufactured solution and gives its amplitude $A$.\n",
+    "- `material` names the constitutive law and gives its parameters.\n",
+    "\n",
+    "Each `type` in the input file names a class."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "06-code-input-file",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "input_file = {\n",
+    "    \"geometry\": {\"type\": \"Beam2D\", \"length\": 1.0, \"width\": 1.0},\n",
+    "    \"solution\": {\"type\": \"Trigonometric\", \"amplitude\": 0.02},\n",
+    "    \"material\": {\"type\": \"LinearElastic\", \"youngModulus\": 1e6, \"poissonRatio\": 0.3},\n",
+    "}"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "07-md-geometry",
+   "metadata": {},
+   "source": [
+    "### The manufactured problem\n",
+    "\n",
+    "The geometry is built first, because the solution depends on it."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "08-code-geometry",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import sys\n",
+    "\n",
+    "# Make the repository's packages importable from this folder.\n",
+    "sys.path.insert(0, \"../..\")\n",
+    "\n",
+    "from verification.case import GEOMETRIES\n",
+    "\n",
+    "geometry_parameters = dict(input_file[\"geometry\"])\n",
+    "geometry = GEOMETRIES[geometry_parameters.pop(\"type\")](**geometry_parameters)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "09-md-solution",
+   "metadata": {},
+   "source": [
+    "The solution's class is looked up by the geometry's dimension: `Trigonometric` in 2D is `Trigonometric2D` (`verification/manufactured/grid.py`). Solutions take the geometry as inputs because various parameters, such as $k$, might depend on the dimensions. Also, the displacement is prescribed on the geometry's regions."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "10-code-solution",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from verification.case import MANUFACTURED_SOLUTIONS\n",
+    "\n",
+    "solution = MANUFACTURED_SOLUTIONS[(geometry.dim, input_file[\"solution\"][\"type\"])](input_file[\"solution\"], geometry)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "11-md-material",
+   "metadata": {},
+   "source": [
+    "The material converts the parameters of the input file into the ones its law is written with, which may depend on the number of spatial dimensions. `LinearElastic` (`verification/manufactured/materials.py`) converts $E$ and $\\nu$ to the Lamé parameters $\\mu$ and $\\lambda$; hyperelastic materials return their own."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "12-code-material",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from verification.case import MATERIALS\n",
+    "\n",
+    "Material = MATERIALS[input_file[\"material\"][\"type\"]]\n",
+    "parameters = Material.get_material_parameters(input_file[\"material\"], geometry.spatial_dimensions)\n",
+    "material = Material(*parameters)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "13-md-problem",
+   "metadata": {},
+   "source": [
+    "The solution and the material are paired into the manufactured problem, which differentiates $\\hat{u}$ symbolically, with `sympy`, into the stress $\\hat{\\sigma}$ and the source term $s$. The problem is `ManufacturedProblem` (`verification/manufactured/problem.py`). It also compiles each expression into a numpy function evaluated at any array of points: `problem.u`, `problem.grad_u`, `problem.stress` and `problem.source`. These work as inputs to apply source terms such as body forces, prescribe displacement on boundaries and traction BCs. Additionally, the error is measured against them during the convergence analysis. These are the expressions it derived, with the numbers of the input file substituted."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "14-code-problem",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import sympy as sp\n",
+    "from IPython.display import Math, display\n",
+    "\n",
+    "from verification.manufactured import ManufacturedProblem\n",
+    "\n",
+    "problem = ManufacturedProblem(solution, material, geometry.spatial_dimensions)\n",
+    "\n",
+    "\n",
+    "def rounded(expression):\n",
+    "    \"\"\"The expression with its floating-point numbers cut to 6 significant digits, for display.\"\"\"\n",
+    "    return expression.xreplace({f: sp.Float(f, 6) for f in expression.atoms(sp.Float)})\n",
+    "\n",
+    "\n",
+    "# The symbolic expressions the problem derived.\n",
+    "display(Math(r\"\\hat{u} = \" + sp.latex(rounded(problem.displacement_expression))),\n",
+    "        Math(r\"\\hat{\\sigma} = \" + sp.latex(rounded(problem.stress_expression))),\n",
+    "        Math(r\"s = \" + sp.latex(rounded(problem.source_expression))))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "15-md-plot-manufactured-solution",
+   "metadata": {},
+   "source": [
+    "The plot below uses `problem.u` to draw the exact $\\hat{u}$ over the square. On the right, the grid of the undeformed square is displaced by $\\hat{u}$: the boundary keeps its shape, but its points slide along it. The prescribed displacement is not zero on the boundary."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "16-code-plot-manufactured-solution",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import matplotlib.pyplot as plt\n",
+    "import numpy as np\n",
+    "\n",
+    "# û is small next to the square; magnify it so the deformation is visible.\n",
+    "MAGNIFICATION = 5\n",
+    "\n",
+    "# Evaluate the compiled û on a fine grid of points covering the square.\n",
+    "x, y = np.meshgrid(np.linspace(0, geometry.length, 201), np.linspace(0, geometry.width, 201))\n",
+    "u = problem.u(np.stack([x, y], axis=-1))\n",
+    "\n",
+    "fig, (magnitude, deformed) = plt.subplots(1, 2, figsize=(11, 4.8), layout=\"constrained\")\n",
+    "\n",
+    "# Left: the magnitude of û as colour, its direction as arrows on a coarser grid.\n",
+    "field = magnitude.pcolormesh(x, y, np.linalg.norm(u, axis=-1), shading=\"gouraud\")\n",
+    "fig.colorbar(field, ax=magnitude, label=r\"$|\\hat{u}|$\")\n",
+    "step = slice(None, None, 20)\n",
+    "magnitude.quiver(x[step, step], y[step, step], u[step, step, 0], u[step, step, 1], color=\"white\")\n",
+    "magnitude.set_title(r\"$\\hat{u}$: magnitude and direction\")\n",
+    "\n",
+    "# Right: each line of the undeformed grid, and the same line moved by the magnified û.\n",
+    "# The boundary lines are drawn thicker.\n",
+    "t = np.linspace(0, 1, 201)\n",
+    "for c in np.linspace(0, 1, 11):\n",
+    "    for line in (np.stack([c * geometry.length + 0 * t, t * geometry.width], axis=-1),\n",
+    "                 np.stack([t * geometry.length, c * geometry.width + 0 * t], axis=-1)):\n",
+    "        moved = line + MAGNIFICATION * problem.u(line)\n",
+    "        deformed.plot(*line.T, color=\"0.8\", linewidth=0.6)\n",
+    "        deformed.plot(*moved.T, color=\"black\", linewidth=1.6 if c in (0, 1) else 0.8)\n",
+    "deformed.set_title(rf\"Undeformed grid (grey) displaced by $\\hat{{u}}$, magnified ×{MAGNIFICATION}\")\n",
+    "\n",
+    "for ax in (magnitude, deformed):\n",
+    "    ax.set_aspect(\"equal\")\n",
+    "    ax.set_xlabel(\"$x$\")\n",
+    "    ax.set_ylabel(\"$y$\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "17-md-solving-keys",
+   "metadata": {},
+   "source": [
+    "### Solving on one mesh\n",
+    "\n",
+    "SOFA now solves the modified problem on a mesh. This needs four more keys for the `input_file`:\n",
+    "\n",
+    "- `forceField` is the SOFA component under test. It is the discrete $\\mathcal{L}$, the operator whose order of accuracy the study verifies;\n",
+    "- `element` is the kind of element the geometry is meshed with, here `quad`. SOFA builds the topology container from it;\n",
+    "- `solvers` sets how SOFA solves the discrete system;\n",
+    "- `sourceQuadratureDegree` is the degree of the quadrature that integrates the source term $s$ into nodal loads."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "18-code-solving-keys",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "input_file[\"forceField\"] = {\"type\": \"LinearSmallStrainFEMForceField\"}\n",
+    "input_file[\"element\"] = \"quad\"\n",
+    "input_file[\"solvers\"] = {\n",
+    "    \"integration\":  {\"type\": \"StaticSolver\"},\n",
+    "    \"newton\":       {\"type\": \"NewtonRaphsonSolver\", \"maxNbIterationsNewton\": 1,\n",
+    "                     \"relativeInitialStoppingThreshold\": 1e-6,\n",
+    "                     \"relativeSuccessiveStoppingThreshold\": 0,\n",
+    "                     \"absoluteResidualStoppingThreshold\": 0,\n",
+    "                     \"printLog\": False},\n",
+    "    \"linearSolver\": {\"type\": \"SparseLDLSolver\", \"template\": \"CompressedRowSparseMatrixd\"},\n",
+    "}\n",
+    "input_file[\"sourceQuadratureDegree\"] = 3"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "19-md-scene",
+   "metadata": {},
+   "source": [
+    "The scene is built in three steps:\n",
+    "\n",
+    "1. **The body.** A geometry knows nothing of SOFA; its prefab holds the instructions to build it in SOFA. `ElasticBeam` (`common/sofa/prefabs/grid.py`), the prefab of `Beam2D`, meshes the square with a regular grid of `quad` elements and adds the force field and the solvers.\n",
+    "2. **The source term.** `ApplyManufacturedSourceTerm` (`verification/scene/controllers.py`) evaluates `problem.source`, the $s$ derived above, at the nodes, and integrates it against the shape functions into nodal loads.\n",
+    "3. **The boundary conditions.** `RegionClamp` (`common/sofa/controllers.py`) adds one constraint per region where the solution prescribes the displacement, `problem.prescribe_displacement_on`, here all four sides. It moves their nodes to `problem.u`, the exact $\\hat{u}$, and fixes them there."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "20-code-scene",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import Sofa.Core\n",
+    "import Sofa.Simulation\n",
+    "\n",
+    "from common.sofa.conventions import VEC_DIM\n",
+    "from common.sofa.controllers import RegionClamp\n",
+    "from common.sofa.prefabs import PREFAB_BY_GEOMETRY\n",
+    "from verification.scene.controllers import ApplyManufacturedSourceTerm\n",
+    "from verification.scene.scene import MMSScene\n",
+    "\n",
+    "vec_type = VEC_DIM[geometry.spatial_dimensions]\n",
+    "\n",
+    "\n",
+    "def solve(cells):\n",
+    "    \"\"\"Build the scene for a mesh of `cells` cells per side and solve it. Returns its root and body.\"\"\"\n",
+    "    root = Sofa.Core.Node(\"root\")\n",
+    "    root.addObject(\"RequiredPlugin\", pluginName=MMSScene.plugins)\n",
+    "    root.addObject(\"DefaultAnimationLoop\")\n",
+    "\n",
+    "    # 1.\n",
+    "    Prefab = PREFAB_BY_GEOMETRY[type(geometry)]\n",
+    "    prefab = root.addChild(Prefab(\n",
+    "        name=Prefab.__name__, geometry=geometry, grid_resolution=[count + 1 for count in cells],\n",
+    "        element=input_file[\"element\"],\n",
+    "        configs={\"material\": input_file[\"material\"], \"forceField\": input_file[\"forceField\"],\n",
+    "                 \"solvers\": input_file[\"solvers\"]}))\n",
+    "    body = prefab.mechanical_node()\n",
+    "\n",
+    "    # 2.\n",
+    "    body.addObject(ApplyManufacturedSourceTerm(\n",
+    "        node=body, dofs=body.dofs, source=problem.source, vec_type=vec_type,\n",
+    "        element=input_file[\"element\"], source_quadrature_degree=input_file[\"sourceQuadratureDegree\"],\n",
+    "        name=\"bodyForceCtrl\"))\n",
+    "\n",
+    "    # 3.\n",
+    "    body.addObject(RegionClamp(\n",
+    "        geometry=geometry, dofs=body.dofs, node=body,\n",
+    "        prescribe_displacement_on=problem.prescribe_displacement_on,\n",
+    "        vec_type=vec_type, displacement=problem.u, name=\"clampCtrl\"))\n",
+    "\n",
+    "    # Initialise the scene, which runs the controllers, then solve it: one step of the static solver.\n",
+    "    Sofa.Simulation.init(root)\n",
+    "    Sofa.Simulation.animate(root, root.dt.value)\n",
+    "    return root, body\n",
+    "\n",
+    "\n",
+    "# The mesh: 8 cells along each side of the square.\n",
+    "cells = [8, 8]\n",
+    "root, body = solve(cells)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "21-md-mms-scene",
+   "metadata": {},
+   "source": [
+    "Steps 1 to 3, building the scene for one mesh, are encapsulated in `MMSScene` (`verification/scene/scene.py`). It only builds the scene, initialization and solution are performed separately by the class that handles the mesh convergence study, introduced below. A side can also be given a traction instead of a displacement; the cases in `verification/cases/` show how."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "22-md-plot-sofa-solution",
+   "metadata": {},
+   "source": [
+    "The plot draws SOFA's solution $u_h$, read off the nodes, over the exact $\\hat{u}$, both magnified as before."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "23-code-plot-sofa-solution",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from matplotlib.collections import PolyCollection\n",
+    "\n",
+    "# u_h at the nodes: where SOFA moved them, minus where they started.\n",
+    "nodes = body.dofs.rest_position.array()\n",
+    "u_h = body.dofs.position.array() - nodes\n",
+    "quads = body.topology.quads.array()\n",
+    "\n",
+    "fig, ax = plt.subplots(figsize=(5.5, 5.5), layout=\"constrained\")\n",
+    "\n",
+    "# The exact û, along the lines of the mesh.\n",
+    "for c in np.linspace(0, 1, cells[0] + 1):\n",
+    "    for line in (np.stack([c * geometry.length + 0 * t, t * geometry.width], axis=-1),\n",
+    "                 np.stack([t * geometry.length, c * geometry.width + 0 * t], axis=-1)):\n",
+    "        moved = line + MAGNIFICATION * problem.u(line)\n",
+    "        ax.plot(*moved.T, color=\"0.6\", linewidth=0.8)\n",
+    "\n",
+    "# SOFA's solution u_h.\n",
+    "ax.add_collection(PolyCollection((nodes + MAGNIFICATION * u_h)[quads],\n",
+    "                                 facecolor=\"none\", edgecolor=\"tab:red\", linewidth=1.2))\n",
+    "\n",
+    "ax.set_aspect(\"equal\")\n",
+    "ax.set_xlabel(\"$x$\")\n",
+    "ax.set_ylabel(\"$y$\")\n",
+    "ax.set_title(rf\"$u_h$ on the {cells[0]}×{cells[1]} mesh (red) and $\\hat{{u}}$ (grey), magnified ×{MAGNIFICATION}\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "24-md-error",
+   "metadata": {},
+   "source": [
+    "### Measuring the error\n",
+    "\n",
+    "The norms $\\|e\\|_{L^2}$ and $|e|_{H^1}$ are integrated by quadrature over each element, with `MeshQuadrature` (`verification/study/quadrature.py`). At every quadrature point, $u_h$ and its gradient are interpolated from the nodes with the element's shape functions, whereas $\\hat{u}$ and $\\nabla\\hat{u}$ are evaluated exactly. One more key sets the quadrature: `errorQuadratureDegree`."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "25-code-error",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from common.sofa.conventions import ELEMENTS\n",
+    "from verification.study.quadrature import MeshQuadrature\n",
+    "\n",
+    "input_file[\"errorQuadratureDegree\"] = 2\n",
+    "\n",
+    "\n",
+    "def error_norms(body):\n",
+    "    \"\"\"‖e‖_L² and |e|_H¹ of the solution SOFA found on `body`'s mesh.\"\"\"\n",
+    "    nodes = body.dofs.rest_position.array()\n",
+    "    u_h = body.dofs.position.array() - nodes\n",
+    "\n",
+    "    # The quadrature points and weights of every element of the mesh.\n",
+    "    quadrature = MeshQuadrature(nodes, body.topology.quads.array(),\n",
+    "                                ELEMENTS[input_file[\"element\"]].cpp,\n",
+    "                                input_file[\"errorQuadratureDegree\"])\n",
+    "\n",
+    "    # The error and its gradient at every quadrature point: u_h interpolated, û evaluated exactly.\n",
+    "    e = quadrature.values(u_h) - quadrature.sample(problem.u)\n",
+    "    grad_e = quadrature.grads(u_h) - quadrature.sample(problem.grad_u)\n",
+    "\n",
+    "    # ‖e‖_L² = √(∫ |e|² dΩ) and |e|_H¹ = √(∫ ‖∇e‖² dΩ).\n",
+    "    return (np.sqrt(quadrature.integrate(np.sum(e ** 2, axis=-1))),\n",
+    "            np.sqrt(quadrature.integrate(np.sum(grad_e ** 2, axis=(-2, -1)))))\n",
+    "\n",
+    "\n",
+    "L2, H1 = error_norms(body)\n",
+    "print(f\"‖e‖_L2 = {L2:.3e}, |e|_H1 = {H1:.3e}\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "26-md-metrics",
+   "metadata": {},
+   "source": [
+    "These two numbers alone do not say whether the code is correct. That is decided by how they decrease as the mesh is refined. The computation of the errors is encapsulated in `Measurement`, and each norm in a metric class: `L2`, `H1` and `EnergyDifference`. `METRICS` lists the metrics every mesh is measured with. All are in `verification/study/metrics.py`."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "27-md-refinement",
+   "metadata": {},
+   "source": [
+    "### Refining the mesh\n",
+    "\n",
+    "The scene is now solved on a sequence of meshes, each $r$ times finer than the one before, and the\n",
+    "observed order $\\hat{p}$ is computed between each mesh and the one before it. The shape functions of `quad` elements\n",
+    "have degree $q = 1$, so they should give $p = 2$ in $L^2$ and $p = 1$ in $H^1$, and $p = 1$ for `EnergyDifference`. Three keys set this:\n",
+    "\n",
+    "- `mesh` gives the coarsest mesh's cells per side, the number of meshes, and $r$;\n",
+    "- `expectedOrder` gives $p$ for each metric;\n",
+    "- `expectedOrderTolerance` is how far $\\hat{p}$ on the finest pair of meshes may sit from $p$."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "28-code-refinement",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# 3×3 cells, then 6×6, 12×12, 24×24 and 48×48.\n",
+    "input_file[\"mesh\"] = {\"cells\": [3, 3], \"levels\": 5, \"refinementRatio\": 2}\n",
+    "input_file[\"expectedOrder\"] = {\"L2\": 2.0, \"H1\": 1.0, \"EnergyDifference\": 1.0}\n",
+    "input_file[\"expectedOrderTolerance\"] = 0.15\n",
+    "\n",
+    "mesh = input_file[\"mesh\"]\n",
+    "spacings, norms = [], []\n",
+    "for level in range(mesh[\"levels\"]):\n",
+    "    level_cells = [count * mesh[\"refinementRatio\"] ** level for count in mesh[\"cells\"]]\n",
+    "    level_root, level_body = solve(level_cells)\n",
+    "    spacings.append(geometry.length / level_cells[0])\n",
+    "    norms.append(error_norms(level_body))\n",
+    "    # Free this mesh's scene before building the next one.\n",
+    "    Sofa.Simulation.unload(level_root)\n",
+    "\n",
+    "# One row per mesh, one column per norm: ‖e‖_L² and |e|_H¹.\n",
+    "norms = np.array(norms)\n",
+    "# p̂ = ln(‖e_rh‖ / ‖e_h‖) / ln r, between each mesh and the one before it.\n",
+    "orders = np.log(norms[:-1] / norms[1:]) / np.log(mesh[\"refinementRatio\"])\n",
+    "\n",
+    "print(f\"{'h':>8} {'L2':>10} {'p̂':>6} {'H1':>10} {'p̂':>6}\")\n",
+    "for level, h in enumerate(spacings):\n",
+    "    rates = orders[level - 1] if level else (np.nan, np.nan)\n",
+    "    print(f\"{h:8.4f} {norms[level, 0]:10.3e} {rates[0]:6.2f} {norms[level, 1]:10.3e} {rates[1]:6.2f}\")\n",
+    "\n",
+    "# The test: p̂ on the finest pair of meshes matches p, within the tolerance.\n",
+    "for name, order in zip((\"L2\", \"H1\"), orders[-1]):\n",
+    "    expected = input_file[\"expectedOrder\"][name]\n",
+    "    passed = abs(order - expected) <= input_file[\"expectedOrderTolerance\"]\n",
+    "    print(f\"{name}: p̂ = {order:.2f}, p = {expected}: {'pass' if passed else 'fail'}\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "29-md-convergence-study",
+   "metadata": {},
+   "source": [
+    "This loop is encapsulated in `ErrorConvergenceStudy` (`verification/study/study.py`). It judges the orders more strictly than the single check above:\n",
+    "\n",
+    "- a rate is discarded when either error is below the noise floor, `relativeNoiseFloor` times the size of the exact field. Below it, the round-off and iterative errors $\\varepsilon_\\text{ro}$ and $\\varepsilon_\\text{it}$ from the start dominate, and a rate measured there is meaningless;\n",
+    "- a rate is also discarded when the errors oscillate or diverge;\n",
+    "- the order is the one on the finest meshes, and is kept only once at least two consecutive rates agree within `asymptoticTolerance`: the meshes are then in the asymptotic regime, and rates from coarser meshes are discarded;\n",
+    "- the study fails if a Newton solve did not converge."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "30-code-convergence-study",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "input_file[\"relativeNoiseFloor\"] = 1e-12\n",
+    "input_file[\"asymptoticTolerance\"] = 0.10"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "31-md-full-input-file",
+   "metadata": {},
+   "source": [
+    "## Putting it together\n",
+    "\n",
+    "The input file is now complete. This is all a study needs, as the JSON file read by `run.py`:\n",
+    "\n",
+    "```json\n",
+    "{\n",
+    "  \"geometry\":    { \"type\": \"Beam2D\", \"length\": 1.0, \"width\": 1.0 },\n",
+    "  \"element\":     \"quad\",\n",
+    "  \"material\":    { \"type\": \"LinearElastic\", \"youngModulus\": 1e6, \"poissonRatio\": 0.3 },\n",
+    "  \"forceField\":  { \"type\": \"LinearSmallStrainFEMForceField\" },\n",
+    "  \"solution\":    { \"type\": \"Trigonometric\", \"amplitude\": 0.02 },\n",
+    "  \"errorQuadratureDegree\": 3,\n",
+    "  \"sourceQuadratureDegree\": 3,\n",
+    "  \"solvers\": {\n",
+    "    \"integration\":  { \"type\": \"StaticSolver\" },\n",
+    "    \"newton\":       { \"type\": \"NewtonRaphsonSolver\", \"maxNbIterationsNewton\": 1,\n",
+    "                      \"relativeInitialStoppingThreshold\": 1e-6,\n",
+    "                      \"relativeSuccessiveStoppingThreshold\": 0,\n",
+    "                      \"absoluteResidualStoppingThreshold\": 0,\n",
+    "                      \"printLog\": false },\n",
+    "    \"linearSolver\": { \"type\": \"SparseLDLSolver\", \"template\": \"CompressedRowSparseMatrixd\" }\n",
+    "  },\n",
+    "  \"mesh\":        { \"cells\": [3, 3], \"levels\": 5, \"refinementRatio\": 2 },\n",
+    "  \"asymptoticTolerance\": 0.10,\n",
+    "  \"expectedOrder\": { \"L2\": 2.0, \"H1\": 1.0, \"EnergyDifference\": 1.0 },\n",
+    "  \"expectedOrderTolerance\": 0.15,\n",
+    "  \"relativeNoiseFloor\": 1e-12\n",
+    "}\n",
+    "```"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "32-md-study",
+   "metadata": {},
+   "source": [
+    "Everything the walkthrough built by hand is done by two classes:\n",
+    "\n",
+    "- `Case` (`verification/case.py`) checks the input file, then builds the geometry, the solution, the material and the manufactured problem from it;\n",
+    "- `ErrorConvergenceStudy` builds the scene of every mesh with `MMSScene` and solves it, then computes the errors of that mesh's solution $u_h$ with `Measurement` and their norms with the metrics listed in `METRICS`, and finally judges the order of each metric. It prints one row per mesh with the rate of each metric, or why it was discarded, then the order each metric settled on next to the expected one.\n",
+    "\n",
+    "It also draws two figures: the error of each metric against $h$, its segments green where the rate was kept and red where it was discarded; and $|u_h|$ on each mesh, which approaches $|\\hat{u}|$ as the mesh is refined."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "33-code-study",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from verification.case import Case\n",
+    "from verification.study.study import ErrorConvergenceStudy\n",
+    "\n",
+    "# The path only names the case; the input file is passed as the dictionary built above.\n",
+    "case = Case(\"LinearSmallStrainFEMForceField/Trigonometric2D_Quad_Dirichlet.json\", input_file)\n",
+    "study = ErrorConvergenceStudy(case)\n",
+    "study.run()\n",
+    "print(\"\\nverified:\", study.verified)\n",
+    "\n",
+    "# The convergence plot the study drew while running: the error of each metric against h.\n",
+    "display(study.plots[\"convergence\"])\n",
+    "\n",
+    "# |u_h| on each mesh, coarsest to finest, drawn on the displaced mesh.\n",
+    "study.plot_fields()\n",
+    "display(study.plots[\"fields\"])"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "34-md-running-a-case",
+   "metadata": {},
+   "source": [
+    "Saved as a JSON file under `verification/cases/<force field>/`, the input file is a *case*, and the\n",
+    "runner does all of the above from the command line, writing the results to `verification/results/`:\n",
+    "\n",
+    "```bash\n",
+    "python verification/run.py LinearSmallStrainFEMForceField/Trigonometric2D_Quad_Dirichlet.json\n",
+    "python verification/run.py --all     # every case, then one summary line per case\n",
+    "```\n",
+    "\n",
+    "To inspect one mesh of a case in the SOFA GUI, the runner also acts as a SOFA scene file:\n",
+    "\n",
+    "```bash\n",
+    "runSofa -l SofaPython3 verification/run.py --argv LinearSmallStrainFEMForceField/Trigonometric2D_Quad_Dirichlet.json --argv 0\n",
+    "```\n",
+    "\n",
+    "The first argument is the path to the .json case file and the second is the mesh refinement level to be built."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "35-md-references",
    "metadata": {},
    "source": [
     "<a id=\"references\"></a>\n",
     "## References\n",
     "\n",
     "<a id=\"ref1\"></a>\n",
-    "**[1]** P. J. Roache, \"Code verification by the method of manufactured solutions,\"\n",
+    "**[1]** S. Steinberg and P. J. Roache, \"Symbolic manipulation and computational fluid dynamics,\"\n",
+    "*Journal of Computational Physics*, 57(2):251&ndash;284, 1985. See also P. J. Roache and\n",
+    "S. Steinberg, *AIAA Journal*, 22(10):1390&ndash;1394, 1984.\n",
+    "\n",
+    "<a id=\"ref2\"></a>\n",
+    "**[2]** P. J. Roache, \"Code verification by the method of manufactured solutions,\"\n",
     "*Journal of Fluids Engineering*, 124(1):4&ndash;10, 2002.\n",
     "[doi:10.1115/1.1436090](https://doi.org/10.1115/1.1436090)\n",
     "\n",
-    "<a id=\"ref2\"></a>\n",
-    "**[2]** C. J. Roy, \"Review of code and solution verification procedures for computational\n",
+    "<a id=\"ref3\"></a>\n",
+    "**[3]** C. J. Roy, \"Review of code and solution verification procedures for computational\n",
     "simulation,\" *Journal of Computational Physics*, 205(1):131&ndash;156, 2005.\n",
     "[doi:10.1016/j.jcp.2004.10.036](https://doi.org/10.1016/j.jcp.2004.10.036)\n",
     "\n",
-    "<a id=\"ref3\"></a>\n",
-    "**[3]** K. Salari and P. Knupp, *Code Verification by the Method of Manufactured Solutions*,\n",
-    "SAND2000-1444, Sandia National Laboratories, Albuquerque, NM, June 2000.\n",
-    "\n",
     "<a id=\"ref4\"></a>\n",
-    "**[4]** S. Étienne, A. Garon and D. Pelletier, \"Some manufactured solutions for verification of\n",
-    "fluid-structure interaction codes,\" *Computers and Structures*, 106&ndash;107:56&ndash;67, 2012.\n",
-    "[doi:10.1016/j.compstruc.2012.04.006](https://doi.org/10.1016/j.compstruc.2012.04.006)\n",
+    "**[4]** K. Salari and P. Knupp, *Code Verification by the Method of Manufactured Solutions*,\n",
+    "SAND2000-1444, Sandia National Laboratories, Albuquerque, NM, June 2000.\n",
     "\n",
     "<a id=\"ref5\"></a>\n",
     "**[5]** W. L. Oberkampf and C. J. Roy, *Verification and Validation in Scientific Computing*,\n",
@@ -162,19 +841,27 @@
     "order verification) and Chapter 6, §6.3 (MMS).*\n",
     "\n",
     "<a id=\"ref6\"></a>\n",
-    "**[6]** K. I. Aycock, N. Rebelo and B. A. Craven, \"Method of manufactured solutions code\n",
-    "verification of elastostatic solid mechanics problems in a commercial finite element solver,\"\n",
-    "*Computers and Structures*, 229:106175, 2020.\n",
-    "[doi:10.1016/j.compstruc.2019.106175](https://doi.org/10.1016/j.compstruc.2019.106175)\n",
-    "\n",
-    "<a id=\"ref7\"></a>\n",
-    "**[7]** W. L. Oberkampf and T. G. Trucano, \"Verification and validation in computational fluid\n",
+    "**[6]** W. L. Oberkampf and T. G. Trucano, \"Verification and validation in computational fluid\n",
     "dynamics,\" *Progress in Aerospace Sciences*, 38(3):209&ndash;272, 2002.\n",
     "[doi:10.1016/S0376-0421(02)00005-2](https://doi.org/10.1016/S0376-0421%2802%2900005-2)\n",
     "\n",
+    "<a id=\"ref7\"></a>\n",
+    "**[7]** B. W. Boehm, *Software Engineering Economics*, Prentice-Hall, Englewood Cliffs, NJ, 1981.\n",
+    "\n",
     "<a id=\"ref8\"></a>\n",
     "**[8]** P. G. Ciarlet, *The Finite Element Method for Elliptic Problems*, North-Holland, Amsterdam,\n",
-    "1978. Reprinted by SIAM, Classics in Applied Mathematics 40, 2002."
+    "1978. Reprinted by SIAM, Classics in Applied Mathematics 40, 2002.\n",
+    "\n",
+    "<a id=\"ref9\"></a>\n",
+    "**[9]** K. I. Aycock, N. Rebelo and B. A. Craven, \"Method of manufactured solutions code\n",
+    "verification of elastostatic solid mechanics problems in a commercial finite element solver,\"\n",
+    "*Computers and Structures*, 229:106175, 2020.\n",
+    "[doi:10.1016/j.compstruc.2019.106175](https://doi.org/10.1016/j.compstruc.2019.106175)\n",
+    "\n",
+    "<a id=\"ref10\"></a>\n",
+    "**[10]** S. Étienne, A. Garon and D. Pelletier, \"Some manufactured solutions for verification of\n",
+    "fluid-structure interaction codes,\" *Computers and Structures*, 106&ndash;107:56&ndash;67, 2012.\n",
+    "[doi:10.1016/j.compstruc.2012.04.006](https://doi.org/10.1016/j.compstruc.2012.04.006)"
    ]
   }
  ],

From a55f8d360a2058f71f8184a1f786b1bb8505aa63 Mon Sep 17 00:00:00 2001
From: Themis Skamagkis <themisskamagkis@gmail.com>
Date: Mon, 12 Oct 2026 09:49:19 +0200
Subject: [PATCH 35/55] [docs] Add the LinearSmallStrainFEMForceField notebook

---
 .../docs/LinearSmallStrainFEMForceField.ipynb | 497 ++++++++++++++++++
 1 file changed, 497 insertions(+)
 create mode 100644 verification/docs/LinearSmallStrainFEMForceField.ipynb

diff --git a/verification/docs/LinearSmallStrainFEMForceField.ipynb b/verification/docs/LinearSmallStrainFEMForceField.ipynb
new file mode 100644
index 0000000..38fe22c
--- /dev/null
+++ b/verification/docs/LinearSmallStrainFEMForceField.ipynb
@@ -0,0 +1,497 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "01-md-title",
+   "metadata": {},
+   "source": [
+    "# Verification of `LinearSmallStrainFEMForceField`\n",
+    "\n",
+    "This notebook verifies SOFA's `LinearSmallStrainFEMForceField` with the Method of Manufactured Solutions, presented in [IntroductionMMS.ipynb](IntroductionMMS.ipynb). It presents the governing PDEs and some of the cases used for verification, in `verification/cases/LinearSmallStrainFEMForceField/`."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "02-md-material",
+   "metadata": {},
+   "source": [
+    "## Governing PDE for Linear Elastic Materials\n",
+    "\n",
+    "Under small-strain deformations, the displacement $\\mathbf{u}$ of a linear elastic body satisfies\n",
+    "\n",
+    "$$\\mathcal{L}(\\mathbf{u}) = -\\nabla \\cdot \\boldsymbol{\\sigma}(\\mathbf{u}) = 0$$\n",
+    "\n",
+    "This is the equation solved by the `LinearSmallStrainFEMForceField` component in SOFA. The body occupies a domain $\\Omega$ in $d = 1$, $2$ or $3$ dimensions, and $\\Gamma$ is its boundary.\n",
+    "\n",
+    "The material is defined by its strain energy density $\\psi$, the elastic energy stored per unit volume, and the stress is its derivative with respect to the displacement gradient. For a linear elastic material, $\\psi$ depends on the strain $\\boldsymbol{\\varepsilon}$:\n",
+    "\n",
+    "$$\\boldsymbol{\\varepsilon}(\\mathbf{u}) = \\frac{1}{2}\\left(\\nabla \\mathbf{u} + \\nabla \\mathbf{u}^T\\right), \\qquad\n",
+    "  \\psi(\\boldsymbol{\\varepsilon}) = \\frac{\\lambda}{2}\\,\\mathrm{tr}(\\boldsymbol{\\varepsilon})^2 + \\mu\\,\\boldsymbol{\\varepsilon} : \\boldsymbol{\\varepsilon}$$\n",
+    "\n",
+    "$$\\boldsymbol{\\sigma}(\\mathbf{u}) = \\frac{\\partial \\psi}{\\partial \\nabla \\mathbf{u}} = \\lambda\\,\\mathrm{tr}(\\boldsymbol{\\varepsilon})\\,\\mathbf{I} + 2\\mu\\,\\boldsymbol{\\varepsilon}$$\n",
+    "\n",
+    "The Lamé parameters $\\lambda$ and $\\mu$ are computed from Young's modulus $E$ and the Poisson ratio $\\nu$:\n",
+    "\n",
+    "$$\\mu = \\frac{E}{2(1 + \\nu)}, \\qquad \\lambda = \\frac{E\\nu}{(1 + \\nu)\\,(1 - (d - 1)\\,\\nu)}$$\n",
+    "\n",
+    "### Manufactured Solution and Boundary Conditions\n",
+    "\n",
+    "For a manufactured solution $\\hat{\\mathbf{u}}$, with stress $\\hat{\\boldsymbol{\\sigma}} = \\boldsymbol{\\sigma}(\\hat{\\mathbf{u}})$, the source term is\n",
+    "\n",
+    "$$\\mathbf{s} = \\mathcal{L}(\\hat{\\mathbf{u}}) = -\\nabla \\cdot \\hat{\\boldsymbol{\\sigma}}$$\n",
+    "\n",
+    "The boundary conditions are taken from $\\hat{\\mathbf{u}}$ as well. The weak form shows which quantities a boundary condition can prescribe [[1]](#ref1). Multiplying $\\mathcal{L}(\\mathbf{u}) = \\mathbf{s}$ by a test function $\\mathbf{v}$, integrating over $\\Omega$ and applying the divergence theorem gives\n",
+    "\n",
+    "$$\\int_\\Omega \\boldsymbol{\\sigma} : \\boldsymbol{\\varepsilon}(\\mathbf{v})\\,d\\Omega - \\int_\\Gamma \\mathbf{t} \\cdot \\mathbf{v}\\,d\\Gamma\n",
+    "  = \\int_\\Omega \\mathbf{s} \\cdot \\mathbf{v}\\,d\\Omega$$\n",
+    "\n",
+    "where $\\mathbf{t} = \\boldsymbol{\\sigma} \\cdot \\mathbf{n}$ is the traction and $\\mathbf{n}$ the outward normal of $\\Gamma$.\n",
+    "\n",
+    "At each point of $\\Gamma$ and in each direction, a boundary condition prescribes one of two quantities:\n",
+    "\n",
+    "- the displacement, $\\mathbf{u} = \\hat{\\mathbf{u}}$ (Dirichlet). The test function vanishes there, and so does the boundary term;\n",
+    "- the traction, $\\mathbf{t} = \\hat{\\mathbf{t}} = \\hat{\\boldsymbol{\\sigma}} \\cdot \\mathbf{n}$ (Neumann), which enters the boundary term as a known load.\n",
+    "\n",
+    "Since $\\hat{\\mathbf{u}}$ is known everywhere, both $\\hat{\\mathbf{u}}$ and $\\hat{\\mathbf{t}}$ are known on the whole of $\\Gamma$. Which one goes where is a choice, and a region of $\\Gamma$ can have its displacement prescribed in one direction and its traction in another; enough of the displacement must be prescribed to keep the body from translating and rotating freely.\n",
+    "\n",
+    "The modified problem is then\n",
+    "\n",
+    "$$-\\nabla \\cdot \\boldsymbol{\\sigma}(\\mathbf{u}) = \\mathbf{s} \\quad \\text{in } \\Omega, \\qquad\n",
+    "  \\mathbf{u} = \\hat{\\mathbf{u}} \\quad \\text{on } \\Gamma_D, \\qquad\n",
+    "  \\boldsymbol{\\sigma}(\\mathbf{u}) \\cdot \\mathbf{n} = \\hat{\\mathbf{t}} \\quad \\text{on } \\Gamma_N$$\n",
+    "\n",
+    "with $\\Gamma_D$ and $\\Gamma_N$ the parts of the boundary where the displacement and the traction are prescribed, and $\\hat{\\mathbf{u}}$ is its exact solution. Each case below picks $\\hat{\\mathbf{u}}$, $\\Omega$, $\\Gamma_D$ and $\\Gamma_N$."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "03-md-setup",
+   "metadata": {},
+   "source": [
+    "## Running the Cases\n",
+    "\n",
+    "The cell below loads the suite's runner and prepares the figures. Run it once before any case: every case cell below then runs on its own, in any order."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "04-code-setup",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Make the plugin's packages importable from verification/docs\n",
+    "import sys\n",
+    "sys.path.insert(0, \"../..\")\n",
+    "\n",
+    "# Hide the runner's printed report; figures still show\n",
+    "import os\n",
+    "sys.stdout = open(os.devnull, \"w\")\n",
+    "\n",
+    "# Load the suite's runner\n",
+    "from verification.run import run\n",
+    "\n",
+    "# For showing figures inline\n",
+    "%matplotlib inline\n",
+    "\n",
+    "# For drawing several convergence figures side by side\n",
+    "import matplotlib.pyplot as plt"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "05-md-1d",
+   "metadata": {},
+   "source": [
+    "## Edge Elements on a 1D Bar\n",
+    "\n",
+    "The bar starts at the origin and has length $L$ along $x$. The displacement $u$ is a scalar."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "06-md-quadratic-1d",
+   "metadata": {},
+   "source": [
+    "### Quadratic Function\n",
+    "\n",
+    "$$\\hat{u} = A\\,x^2, \\qquad \\hat{\\sigma} = 2EAx, \\qquad s = -2EA$$\n",
+    "\n",
+    "- `left` ($x = 0$): $\\hat{u} = 0$\n",
+    "- `right` ($x = L$): $\\hat{t} = 2EAL$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "07-code-quadratic-1d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Case: 1D bar, quadratic manufactured solution\n",
+    "study = run(\"../cases/LinearSmallStrainFEMForceField/Quadratic1D.json\")\n",
+    "# u_h along the bar at each level, coarse to fine\n",
+    "study.plot_fields()\n",
+    "study.plots[\"convergence\"]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "08-md-trigonometric-1d",
+   "metadata": {},
+   "source": [
+    "### Trigonometric Function\n",
+    "\n",
+    "$$\\hat{u} = A\\sin(kx), \\qquad \\hat{\\sigma} = EAk\\cos(kx), \\qquad s = EAk^2\\sin(kx)$$\n",
+    "\n",
+    "with $k = \\frac{2\\pi}{L}$.\n",
+    "\n",
+    "- `left` ($x = 0$): $\\hat{u} = 0$\n",
+    "- `right` ($x = L$): $\\hat{t} = EAk$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "09-code-trigonometric-1d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Case: 1D bar, trigonometric manufactured solution\n",
+    "study = run(\"../cases/LinearSmallStrainFEMForceField/Trigonometric1D.json\")\n",
+    "# u_h along the bar at each level, coarse to fine\n",
+    "study.plot_fields()\n",
+    "study.plots[\"convergence\"]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "10-md-2d",
+   "metadata": {},
+   "source": [
+    "## Triangle and Quadrilateral Elements on a 2D Beam\n",
+    "\n",
+    "The beam starts at the origin, with length $L$ along $x$ and width $W$ along $y$."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "11-md-quadratic-2d",
+   "metadata": {},
+   "source": [
+    "### Quadratic Function\n",
+    "\n",
+    "$$\n",
+    "\\hat{\\mathbf{u}} = A\\begin{pmatrix} x^2 + xy \\\\ y^2 + xy \\end{pmatrix}, \n",
+    "\\qquad \\mathbf{s} = -A \\begin{pmatrix} 3\\lambda + 5\\mu \\\\ 3\\lambda + 5\\mu \\end{pmatrix}, \n",
+    "\\qquad \\hat{\\boldsymbol{\\sigma}} = A\\begin{pmatrix} 3\\lambda(x + y) + 2\\mu(2x + y) & \\mu(x + y) \\\\ \\mu(x + y) & 3\\lambda(x + y) + 2\\mu(x + 2y) \\end{pmatrix}\n",
+    "$$\n",
+    "\n",
+    "- `left` ($x = 0$): $\\hat{u}_x = 0$, $\\hat{t}_y = -\\mu A y$\n",
+    "- `bottom` ($y = 0$): $\\hat{u}_y = 0$, $\\hat{t}_x = -\\mu A x$\n",
+    "- `right` ($x = L$): $\\hat{\\mathbf{t}} = A\\,\\big(3\\lambda(L + y) + 2\\mu(2L + y),\\ \\mu(L + y)\\big)$\n",
+    "- `top` ($y = W$): $\\hat{\\mathbf{t}} = A\\,\\big(\\mu(x + W),\\ 3\\lambda(x + W) + 2\\mu(x + 2W)\\big)$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "12-code-quadratic-2d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Case: 2D rectangle, quadratic manufactured solution, triangles (left) and quadrilaterals (right)\n",
+    "tri = run(\"../cases/LinearSmallStrainFEMForceField/Quadratic2D_Tri.json\")\n",
+    "quad = run(\"../cases/LinearSmallStrainFEMForceField/Quadratic2D_Quad.json\")\n",
+    "figure, (triangles, quadrilaterals) = plt.subplots(1, 2, figsize=(14.4, 4.8))\n",
+    "tri.plot(triangles)\n",
+    "quad.plot(quadrilaterals)\n",
+    "\n",
+    "# Deformed shape at each level, coarse to fine\n",
+    "tri.plot_fields()\n",
+    "quad.plot_fields()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "13-md-trigonometric-2d",
+   "metadata": {},
+   "source": [
+    "### Trigonometric Function\n",
+    "\n",
+    "$$\\hat{\\mathbf{u}} = A\\begin{pmatrix} \\sin(kx)\\cos(ky) \\\\ \\cos(kx)\\sin(ky) \\end{pmatrix}, \\qquad \\mathbf{s} = 2k^2(\\lambda + 2\\mu)\\,\\hat{\\mathbf{u}}, \\qquad \\hat{\\boldsymbol{\\sigma}} = 2Ak\\begin{pmatrix} (\\lambda + \\mu)\\cos(kx)\\cos(ky) & -\\mu\\sin(kx)\\sin(ky) \\\\ -\\mu\\sin(kx)\\sin(ky) & (\\lambda + \\mu)\\cos(kx)\\cos(ky) \\end{pmatrix}$$\n",
+    "\n",
+    "with $k = \\frac{2\\pi}{L}$, the same along $x$ and $y$ since the cases use a square, $L = W$.\n",
+    "\n",
+    "- `left` ($x = 0$): $\\hat{u}_x = 0$, $\\hat{t}_y = 0$\n",
+    "- `bottom` ($y = 0$): $\\hat{u}_y = 0$, $\\hat{t}_x = 0$\n",
+    "- `right` ($x = L$): $\\hat{\\mathbf{t}} = 2Ak(\\lambda + \\mu)\\,\\big(\\cos(ky),\\ 0\\big)$\n",
+    "- `top` ($y = L$): $\\hat{\\mathbf{t}} = 2Ak(\\lambda + \\mu)\\,\\big(0,\\ \\cos(kx)\\big)$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "14-code-trigonometric-2d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Case: 2D rectangle, trigonometric manufactured solution, triangles (left) and quadrilaterals (right)\n",
+    "tri = run(\"../cases/LinearSmallStrainFEMForceField/Trigonometric2D_Tri.json\")\n",
+    "quad = run(\"../cases/LinearSmallStrainFEMForceField/Trigonometric2D_Quad.json\")\n",
+    "figure, (triangles, quadrilaterals) = plt.subplots(1, 2, figsize=(14.4, 4.8))\n",
+    "tri.plot(triangles)\n",
+    "quad.plot(quadrilaterals)\n",
+    "\n",
+    "# Deformed shape at each level, coarse to fine\n",
+    "tri.plot_fields()\n",
+    "quad.plot_fields()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "15-md-trigonometric-2d-plane-strain",
+   "metadata": {},
+   "source": [
+    "#### Plane Strain\n",
+    "\n",
+    "The beam is solved in 3D space with $u_z = 0$. $\\lambda$ and $\\mu$ take their 3D values. The displacement is prescribed on all four sides:\n",
+    "\n",
+    "- `left` ($x = 0$) and `right` ($x = L$): $\\hat{\\mathbf{u}} = A\\,\\big(0,\\ \\sin(ky),\\ 0\\big)$\n",
+    "- `bottom` ($y = 0$) and `top` ($y = L$): $\\hat{\\mathbf{u}} = A\\,\\big(\\sin(kx),\\ 0,\\ 0\\big)$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "16-code-trigonometric-2d-plane-strain-tri",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Case: 2D rectangle in plane strain, trigonometric manufactured solution, triangles\n",
+    "study = run(\"../cases/LinearSmallStrainFEMForceField/Trigonometric2D_Tri_PlaneStrain.json\")\n",
+    "# Deformed shape at each level, coarse to fine\n",
+    "study.plot_fields()\n",
+    "study.plots[\"convergence\"]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "17-md-2d-tire",
+   "metadata": {},
+   "source": [
+    "## Triangle and Quadrilateral Elements on a 2D Tire\n",
+    "\n",
+    "The tire is an annulus centred at the origin, with inner radius $R_i$ and outer radius $R_o$, adjusted to 2D from [[2]](#ref2)."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "18-md-quartic-2d",
+   "metadata": {},
+   "source": [
+    "### Quartic Function\n",
+    "\n",
+    "$$\\hat{\\mathbf{u}} = A\\begin{pmatrix} x^4 + \\frac{2}{5}xy \\\\ \\frac{1}{10}y^4 - \\frac{2}{5}xy \\end{pmatrix}, \\qquad \\mathbf{s} = -\\frac{2A}{5}\\begin{pmatrix} 30(\\lambda + 2\\mu)x^2 - (\\lambda + \\mu) \\\\ 3(\\lambda + 2\\mu)y^2 + (\\lambda + \\mu) \\end{pmatrix}, \\qquad \\hat{\\boldsymbol{\\sigma}} = A\\begin{pmatrix} \\lambda\\theta + 2\\mu\\left(4x^3 + \\frac{2}{5}y\\right) & \\frac{2}{5}\\mu(x - y) \\\\ \\frac{2}{5}\\mu(x - y) & \\lambda\\theta + \\frac{4}{5}\\mu\\left(y^3 - x\\right) \\end{pmatrix}$$\n",
+    "\n",
+    "with $\\theta = 4x^3 + \\frac{2}{5}\\left(y^3 + y - x\\right)$.\n",
+    "\n",
+    "The boundary is the inner and outer rims of the annulus:\n",
+    "\n",
+    "- `inner` ($r = R_i$): Dirichlet, $\\mathbf{u} = \\hat{\\mathbf{u}}$\n",
+    "- `outer` ($r = R_o$): Neumann, $\\mathbf{t} = \\hat{\\boldsymbol{\\sigma}} \\cdot \\mathbf{n}$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "19-code-quartic-2d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Case: 2D tire, quartic manufactured solution, triangles (left) and quadrilaterals (right)\n",
+    "tri = run(\"../cases/LinearSmallStrainFEMForceField/Quartic2D_Tri.json\")\n",
+    "quad = run(\"../cases/LinearSmallStrainFEMForceField/Quartic2D_Quad.json\")\n",
+    "figure, (triangles, quadrilaterals) = plt.subplots(1, 2, figsize=(14.4, 4.8))\n",
+    "tri.plot(triangles)\n",
+    "quad.plot(quadrilaterals)\n",
+    "\n",
+    "# Deformed shape at each level, coarse to fine\n",
+    "tri.plot_fields()\n",
+    "quad.plot_fields()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "20-md-3d",
+   "metadata": {},
+   "source": [
+    "## Tetrahedron and Hexahedron Elements on a 3D Beam\n",
+    "\n",
+    "The beam starts at the origin, with length $L$ along $x$, width $W$ along $y$ and height $H$ along $z$."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "21-md-quadratic-3d",
+   "metadata": {},
+   "source": [
+    "### Quadratic Function\n",
+    "\n",
+    "$$\\hat{\\mathbf{u}} = A\\begin{pmatrix} x^2 + xy + xz \\\\ y^2 + yz + xy \\\\ z^2 + zx + yz \\end{pmatrix}, \\qquad \\mathbf{s} = -2A\\begin{pmatrix} 2\\lambda + 3\\mu \\\\ 2\\lambda + 3\\mu \\\\ 2\\lambda + 3\\mu \\end{pmatrix}, \\qquad \\hat{\\boldsymbol{\\sigma}} = A\\begin{pmatrix} 4\\lambda\\theta + 2\\mu(2x + y + z) & \\mu(x + y) & \\mu(x + z) \\\\ \\mu(x + y) & 4\\lambda\\theta + 2\\mu(x + 2y + z) & \\mu(y + z) \\\\ \\mu(x + z) & \\mu(y + z) & 4\\lambda\\theta + 2\\mu(x + y + 2z) \\end{pmatrix}$$\n",
+    "\n",
+    "with $\\theta = x + y + z$.\n",
+    "\n",
+    "- `left` ($x = 0$): $\\hat{u}_x = 0$, $\\hat{t}_y = -\\mu A y$, $\\hat{t}_z = -\\mu A z$\n",
+    "- `bottom` ($y = 0$): $\\hat{u}_y = 0$, $\\hat{t}_x = -\\mu A x$, $\\hat{t}_z = -\\mu A z$\n",
+    "- `front` ($z = 0$): $\\hat{u}_z = 0$, $\\hat{t}_x = -\\mu A x$, $\\hat{t}_y = -\\mu A y$\n",
+    "- `right` ($x = L$): $\\hat{\\mathbf{t}} = A\\,\\big(4\\lambda(L + y + z) + 2\\mu(2L + y + z),\\ \\mu(L + y),\\ \\mu(L + z)\\big)$\n",
+    "- `top` ($y = W$): $\\hat{\\mathbf{t}} = A\\,\\big(\\mu(x + W),\\ 4\\lambda(x + W + z) + 2\\mu(x + 2W + z),\\ \\mu(W + z)\\big)$\n",
+    "- `back` ($z = H$): $\\hat{\\mathbf{t}} = A\\,\\big(\\mu(x + H),\\ \\mu(y + H),\\ 4\\lambda(x + y + H) + 2\\mu(x + y + 2H)\\big)$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "22-code-quadratic-3d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Case: 3D box, quadratic manufactured solution, tetrahedra (left) and hexahedra (right)\n",
+    "tet = run(\"../cases/LinearSmallStrainFEMForceField/Quadratic3D_Tetra.json\")\n",
+    "hexa = run(\"../cases/LinearSmallStrainFEMForceField/Quadratic3D_Hexa.json\")\n",
+    "figure, (tetrahedra, hexahedra) = plt.subplots(1, 2, figsize=(14.4, 4.8))\n",
+    "tet.plot(tetrahedra)\n",
+    "hexa.plot(hexahedra)\n",
+    "\n",
+    "# Deformed shape at each level, coarse to fine\n",
+    "tet.plot_fields()\n",
+    "hexa.plot_fields()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "23-md-trigonometric-3d",
+   "metadata": {},
+   "source": [
+    "### Trigonometric Function\n",
+    "\n",
+    "$$\\hat{\\mathbf{u}} = A\\begin{pmatrix} \\sin(kx)\\cos(ky)\\cos(kz) \\\\ \\cos(kx)\\sin(ky)\\cos(kz) \\\\ \\cos(kx)\\cos(ky)\\sin(kz) \\end{pmatrix}, \\qquad \\mathbf{s} = 3k^2(\\lambda + 2\\mu)\\,\\hat{\\mathbf{u}}, \\qquad \\hat{\\boldsymbol{\\sigma}} = Ak\\begin{pmatrix} (3\\lambda + 2\\mu)\\,\\cos(kx)\\cos(ky)\\cos(kz) & -2\\mu\\,\\sin(kx)\\sin(ky)\\cos(kz) & -2\\mu\\,\\sin(kx)\\cos(ky)\\sin(kz) \\\\ -2\\mu\\,\\sin(kx)\\sin(ky)\\cos(kz) & (3\\lambda + 2\\mu)\\,\\cos(kx)\\cos(ky)\\cos(kz) & -2\\mu\\,\\cos(kx)\\sin(ky)\\sin(kz) \\\\ -2\\mu\\,\\sin(kx)\\cos(ky)\\sin(kz) & -2\\mu\\,\\cos(kx)\\sin(ky)\\sin(kz) & (3\\lambda + 2\\mu)\\,\\cos(kx)\\cos(ky)\\cos(kz) \\end{pmatrix}$$\n",
+    "\n",
+    "with $k = \\frac{2\\pi}{L}$, the same along $x$, $y$ and $z$ since the cases use a cube, $L = W = H$.\n",
+    "\n",
+    "- `left` ($x = 0$): $\\hat{u}_x = 0$, $\\hat{t}_y = \\hat{t}_z = 0$\n",
+    "- `bottom` ($y = 0$): $\\hat{u}_y = 0$, $\\hat{t}_x = \\hat{t}_z = 0$\n",
+    "- `front` ($z = 0$): $\\hat{u}_z = 0$, $\\hat{t}_x = \\hat{t}_y = 0$\n",
+    "- `right` ($x = L$): $\\hat{\\mathbf{t}} = Ak(3\\lambda + 2\\mu)\\,\\big(\\cos(ky)\\cos(kz),\\ 0,\\ 0\\big)$\n",
+    "- `top` ($y = L$): $\\hat{\\mathbf{t}} = Ak(3\\lambda + 2\\mu)\\,\\big(0,\\ \\cos(kx)\\cos(kz),\\ 0\\big)$\n",
+    "- `back` ($z = L$): $\\hat{\\mathbf{t}} = Ak(3\\lambda + 2\\mu)\\,\\big(0,\\ 0,\\ \\cos(kx)\\cos(ky)\\big)$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "24-code-trigonometric-3d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Case: 3D box, trigonometric manufactured solution, tetrahedra (left) and hexahedra (right)\n",
+    "tet = run(\"../cases/LinearSmallStrainFEMForceField/Trigonometric3D_Tetra.json\")\n",
+    "hexa = run(\"../cases/LinearSmallStrainFEMForceField/Trigonometric3D_Hexa.json\")\n",
+    "figure, (tetrahedra, hexahedra) = plt.subplots(1, 2, figsize=(14.4, 4.8))\n",
+    "tet.plot(tetrahedra)\n",
+    "hexa.plot(hexahedra)\n",
+    "\n",
+    "# Deformed shape at each level, coarse to fine\n",
+    "tet.plot_fields()\n",
+    "hexa.plot_fields()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "25-md-3d-tire",
+   "metadata": {},
+   "source": [
+    "## Tetrahedron and Hexahedron Elements on a 3D Tire\n",
+    "\n",
+    "The tire is the annulus of the 2D tire, extruded along $z$ from $0$ to $H$."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "26-md-quartic-3d",
+   "metadata": {},
+   "source": [
+    "### Quartic Function\n",
+    "\n",
+    "The manufactured solution is the one of [[2]](#ref2):\n",
+    "\n",
+    "$$\\hat{\\mathbf{u}} = A\\begin{pmatrix} x^4 - \\frac{2}{5}yz \\\\ -\\frac{1}{10}y^4 - \\frac{2}{5}xyz \\\\ z^4 - \\frac{2}{5}xy \\end{pmatrix}, \\qquad \\mathbf{s} = A\\begin{pmatrix} -12(\\lambda + 2\\mu)x^2 + \\frac{2}{5}(\\lambda + \\mu)z \\\\ \\frac{6}{5}(\\lambda + 2\\mu)y^2 \\\\ -12(\\lambda + 2\\mu)z^2 + \\frac{2}{5}(\\lambda + \\mu)x \\end{pmatrix}, \\qquad \\hat{\\boldsymbol{\\sigma}} = A\\begin{pmatrix} \\lambda\\theta + 8\\mu x^3 & -\\frac{2}{5}\\mu z(1 + y) & -\\frac{4}{5}\\mu y \\\\ -\\frac{2}{5}\\mu z(1 + y) & \\lambda\\theta - \\frac{4}{5}\\mu\\left(y^3 + xz\\right) & -\\frac{2}{5}\\mu x(1 + y) \\\\ -\\frac{4}{5}\\mu y & -\\frac{2}{5}\\mu x(1 + y) & \\lambda\\theta + 8\\mu z^3 \\end{pmatrix}$$\n",
+    "\n",
+    "with $\\theta = 4x^3 + 4z^3 - \\frac{2}{5}y^3 - \\frac{2}{5}xz$.\n",
+    "\n",
+    "The boundary is the inner and outer curved surfaces and the two flat faces:\n",
+    "\n",
+    "- `inner` ($r = R_i$): Dirichlet, $\\mathbf{u} = \\hat{\\mathbf{u}}$\n",
+    "- `outer` ($r = R_o$): Dirichlet, $\\mathbf{u} = \\hat{\\mathbf{u}}$\n",
+    "- `front` ($z = 0$): Neumann, $\\mathbf{t} = \\hat{\\boldsymbol{\\sigma}} \\cdot \\mathbf{n}$\n",
+    "- `back` ($z = H$): Neumann, $\\mathbf{t} = \\hat{\\boldsymbol{\\sigma}} \\cdot \\mathbf{n}$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "27-code-quartic-3d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Case: 3D tire, quartic manufactured solution, tetrahedra (left) and hexahedra (right)\n",
+    "tet = run(\"../cases/LinearSmallStrainFEMForceField/Quartic3D_Tetra.json\")\n",
+    "hexa = run(\"../cases/LinearSmallStrainFEMForceField/Quartic3D_Hexa.json\")\n",
+    "figure, (tetrahedra, hexahedra) = plt.subplots(1, 2, figsize=(14.4, 4.8))\n",
+    "tet.plot(tetrahedra)\n",
+    "hexa.plot(hexahedra)\n",
+    "\n",
+    "# Deformed shape at each level, coarse to fine\n",
+    "tet.plot_fields()\n",
+    "hexa.plot_fields()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "28-md-references",
+   "metadata": {},
+   "source": [
+    "<a id=\"references\"></a>\n",
+    "## References\n",
+    "\n",
+    "<a id=\"ref1\"></a>\n",
+    "**[1]** T. J. R. Hughes, *The Finite Element Method: Linear Static and Dynamic Finite Element Analysis*, Dover, Mineola, NY, 2000. ISBN 978-0-486-41181-1. &mdash; *Chapter 2 (the weak formulation of linear elasticity).*\n",
+    "\n",
+    "<a id=\"ref2\"></a>\n",
+    "**[2]** É. Chamberland, A. Fortin and M. Fortin, \"Comparison of the performance of some finite element discretizations for large deformation elasticity problems,\" *Computers and Structures*, 88(11–12):664–673, 2010."
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.12.12"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

From d46f031522ccaf8282979231afaaccc3f6d8101d Mon Sep 17 00:00:00 2001
From: Themis Skamagkis <themisskamagkis@gmail.com>
Date: Mon, 12 Oct 2026 09:49:19 +0200
Subject: [PATCH 36/55] [verification] Add the
 LinearSmallStrainFEMForceField/Quadratic1D case

---
 .../Quadratic1D.json                          | 24 +++++++++++++++++++
 1 file changed, 24 insertions(+)
 create mode 100644 verification/cases/LinearSmallStrainFEMForceField/Quadratic1D.json

diff --git a/verification/cases/LinearSmallStrainFEMForceField/Quadratic1D.json b/verification/cases/LinearSmallStrainFEMForceField/Quadratic1D.json
new file mode 100644
index 0000000..2cfd0c2
--- /dev/null
+++ b/verification/cases/LinearSmallStrainFEMForceField/Quadratic1D.json
@@ -0,0 +1,24 @@
+{
+  "geometry":    { "type": "Bar1D", "length": 2.0 },
+  "element":     "edge",
+  "material":    { "type": "LinearElastic", "youngModulus": 1e6, "poissonRatio": 0.3 },
+  "forceField":  { "type": "LinearSmallStrainFEMForceField" },
+  "tractionOn":  { "right": [1] },
+  "solution":    { "type": "Quadratic", "amplitude": 0.1 },
+  "errorQuadratureDegree": 3,
+  "sourceQuadratureDegree": 2,
+  "solvers": {
+    "integration":  { "type": "StaticSolver" },
+    "newton":       { "type": "NewtonRaphsonSolver", "maxNbIterationsNewton": 1,
+                      "relativeInitialStoppingThreshold": 1e-6,
+                      "relativeSuccessiveStoppingThreshold": 0,
+                      "absoluteResidualStoppingThreshold": 0,
+                      "printLog": false },
+    "linearSolver": { "type": "SparseLDLSolver", "template": "CompressedRowSparseMatrixd" }
+  },
+  "mesh":        { "cells": [3], "levels": 7, "refinementRatio": 2 },
+  "asymptoticTolerance": 0.10,
+  "expectedOrder": { "L2": 2.0, "H1": 1.0, "EnergyDifference": 1.0 },
+  "expectedOrderTolerance": 0.15,
+  "relativeNoiseFloor": 1e-12
+}

From ce162aeddd1c0a5a0746ba0ea7c6fcef45f0ba1d Mon Sep 17 00:00:00 2001
From: Themis Skamagkis <themisskamagkis@gmail.com>
Date: Mon, 12 Oct 2026 09:49:19 +0200
Subject: [PATCH 37/55] [verification] Add the
 LinearSmallStrainFEMForceField/Quadratic2D_Quad case

---
 .../Quadratic2D_Quad.json                     | 24 +++++++++++++++++++
 1 file changed, 24 insertions(+)
 create mode 100644 verification/cases/LinearSmallStrainFEMForceField/Quadratic2D_Quad.json

diff --git a/verification/cases/LinearSmallStrainFEMForceField/Quadratic2D_Quad.json b/verification/cases/LinearSmallStrainFEMForceField/Quadratic2D_Quad.json
new file mode 100644
index 0000000..6b3824c
--- /dev/null
+++ b/verification/cases/LinearSmallStrainFEMForceField/Quadratic2D_Quad.json
@@ -0,0 +1,24 @@
+{
+  "geometry":    { "type": "Beam2D", "length": 1.0, "width": 1.0 },
+  "element":     "quad",
+  "material":    { "type": "LinearElastic", "youngModulus": 1e6, "poissonRatio": 0.3 },
+  "forceField":  { "type": "LinearSmallStrainFEMForceField" },
+  "tractionOn":  { "left": [0, 1], "bottom": [1, 0], "right": [1, 1], "top": [1, 1] },
+  "solution":    { "type": "Quadratic", "amplitude": 0.1 },
+  "errorQuadratureDegree": 3,
+  "sourceQuadratureDegree": 3,
+  "solvers": {
+    "integration":  { "type": "StaticSolver" },
+    "newton":       { "type": "NewtonRaphsonSolver", "maxNbIterationsNewton": 1,
+                      "relativeInitialStoppingThreshold": 1e-6,
+                      "relativeSuccessiveStoppingThreshold": 0,
+                      "absoluteResidualStoppingThreshold": 0,
+                      "printLog": false },
+    "linearSolver": { "type": "SparseLDLSolver", "template": "CompressedRowSparseMatrixd" }
+  },
+  "mesh":        { "cells": [3, 3], "levels": 6, "refinementRatio": 2 },
+  "asymptoticTolerance": 0.10,
+  "expectedOrder": { "L2": 2.0, "H1": 1.0, "EnergyDifference": 1.0 },
+  "expectedOrderTolerance": 0.15,
+  "relativeNoiseFloor": 1e-12
+}

From e566446f3b86f89aab615c068a1f0b81771d5c5c Mon Sep 17 00:00:00 2001
From: Themis Skamagkis <themisskamagkis@gmail.com>
Date: Mon, 12 Oct 2026 09:49:19 +0200
Subject: [PATCH 38/55] [verification] Add the
 LinearSmallStrainFEMForceField/Quadratic2D_Tri case

---
 .../Quadratic2D_Tri.json                      | 24 +++++++++++++++++++
 1 file changed, 24 insertions(+)
 create mode 100644 verification/cases/LinearSmallStrainFEMForceField/Quadratic2D_Tri.json

diff --git a/verification/cases/LinearSmallStrainFEMForceField/Quadratic2D_Tri.json b/verification/cases/LinearSmallStrainFEMForceField/Quadratic2D_Tri.json
new file mode 100644
index 0000000..a65c2c0
--- /dev/null
+++ b/verification/cases/LinearSmallStrainFEMForceField/Quadratic2D_Tri.json
@@ -0,0 +1,24 @@
+{
+  "geometry":    { "type": "Beam2D", "length": 1.0, "width": 1.0 },
+  "element":     "tri",
+  "material":    { "type": "LinearElastic", "youngModulus": 1e6, "poissonRatio": 0.3 },
+  "forceField":  { "type": "LinearSmallStrainFEMForceField" },
+  "tractionOn":  { "left": [0, 1], "bottom": [1, 0], "right": [1, 1], "top": [1, 1] },
+  "solution":    { "type": "Quadratic", "amplitude": 0.1 },
+  "errorQuadratureDegree": 2,
+  "sourceQuadratureDegree": 2,
+  "solvers": {
+    "integration":  { "type": "StaticSolver" },
+    "newton":       { "type": "NewtonRaphsonSolver", "maxNbIterationsNewton": 1,
+                      "relativeInitialStoppingThreshold": 1e-6,
+                      "relativeSuccessiveStoppingThreshold": 0,
+                      "absoluteResidualStoppingThreshold": 0,
+                      "printLog": false },
+    "linearSolver": { "type": "SparseLDLSolver", "template": "CompressedRowSparseMatrixd" }
+  },
+  "mesh":        { "cells": [3, 3], "levels": 6, "refinementRatio": 2 },
+  "asymptoticTolerance": 0.10,
+  "expectedOrder": { "L2": 2.0, "H1": 1.0, "EnergyDifference": 1.0 },
+  "expectedOrderTolerance": 0.15,
+  "relativeNoiseFloor": 1e-12
+}

From 1849b9e1dbb527aa338e3a3902ea171e5eb23ee5 Mon Sep 17 00:00:00 2001
From: Themis Skamagkis <themisskamagkis@gmail.com>
Date: Mon, 12 Oct 2026 09:49:19 +0200
Subject: [PATCH 39/55] [verification] Add the
 LinearSmallStrainFEMForceField/Quadratic3D_Hexa case

---
 .../Quadratic3D_Hexa.json                     | 31 +++++++++++++++++++
 1 file changed, 31 insertions(+)
 create mode 100644 verification/cases/LinearSmallStrainFEMForceField/Quadratic3D_Hexa.json

diff --git a/verification/cases/LinearSmallStrainFEMForceField/Quadratic3D_Hexa.json b/verification/cases/LinearSmallStrainFEMForceField/Quadratic3D_Hexa.json
new file mode 100644
index 0000000..70df3b7
--- /dev/null
+++ b/verification/cases/LinearSmallStrainFEMForceField/Quadratic3D_Hexa.json
@@ -0,0 +1,31 @@
+{
+  "geometry":    { "type": "Beam3D", "length": 1.0, "width": 1.0, "height": 1.0 },
+  "element":     "hexa",
+  "material":    { "type": "LinearElastic", "youngModulus": 1e6, "poissonRatio": 0.3 },
+  "forceField":  { "type": "LinearSmallStrainFEMForceField" },
+  "tractionOn":  { "left": [0, 1, 1], "bottom": [1, 0, 1], "front": [1, 1, 0], "right": [1, 1, 1], "top": [1, 1, 1], "back": [1, 1, 1] },
+  "solution":    { "type": "Quadratic", "amplitude": 0.1 },
+  "errorQuadratureDegree": 3,
+  "sourceQuadratureDegree": 3,
+  "solvers": {
+    "integration":  { "type": "StaticSolver" },
+    "newton":       { "type": "NewtonRaphsonSolver", "maxNbIterationsNewton": 1,
+                      "relativeInitialStoppingThreshold": 1e-6,
+                      "relativeSuccessiveStoppingThreshold": 0,
+                      "absoluteResidualStoppingThreshold": 0,
+                      "printLog": false },
+    "linearSolver": {
+      "type": "PCGLinearSolver", "iterations": 1000, "tolerance": 1e-24, "printLog": false,
+      "system": { "type": "PreconditionedMatrixFreeSystem" },
+      "preconditioner": {
+        "type": "SSORPreconditioner",
+        "system": { "type": "MatrixLinearSystem", "template": "CompressedRowSparseMatrixd" }
+      }
+    }
+  },
+  "mesh":        { "cells": [3, 3, 3], "levels": 5, "refinementRatio": 2 },
+  "asymptoticTolerance": 0.10,
+  "expectedOrder": { "L2": 2.0, "H1": 1.0, "EnergyDifference": 1.0 },
+  "expectedOrderTolerance": 0.15,
+  "relativeNoiseFloor": 1e-10
+}

From aa44bf4784924a076e9e44d00a1c1dc2fe20539f Mon Sep 17 00:00:00 2001
From: Themis Skamagkis <themisskamagkis@gmail.com>
Date: Mon, 12 Oct 2026 09:49:19 +0200
Subject: [PATCH 40/55] [verification] Add the
 LinearSmallStrainFEMForceField/Quadratic3D_Tetra case

---
 .../Quadratic3D_Tetra.json                    | 32 +++++++++++++++++++
 1 file changed, 32 insertions(+)
 create mode 100644 verification/cases/LinearSmallStrainFEMForceField/Quadratic3D_Tetra.json

diff --git a/verification/cases/LinearSmallStrainFEMForceField/Quadratic3D_Tetra.json b/verification/cases/LinearSmallStrainFEMForceField/Quadratic3D_Tetra.json
new file mode 100644
index 0000000..fa3ae31
--- /dev/null
+++ b/verification/cases/LinearSmallStrainFEMForceField/Quadratic3D_Tetra.json
@@ -0,0 +1,32 @@
+{
+  "geometry":    { "type": "Beam3D", "length": 1.0, "width": 1.0, "height": 1.0 },
+  "element":     "tet",
+  "material":    { "type": "LinearElastic", "youngModulus": 1e6, "poissonRatio": 0.3 },
+  "forceField":  { "type": "LinearSmallStrainFEMForceField" },
+  "tractionOn":  { "left": [0, 1, 1], "bottom": [1, 0, 1], "front": [1, 1, 0], "right": [1, 1, 1], "top": [1, 1, 1], "back": [1, 1, 1] },
+  "solution":    { "type": "Quadratic", "amplitude": 0.1 },
+  "errorQuadratureDegree": 2,
+  "sourceQuadratureDegree": 2,
+  "solvers": {
+    "integration":  { "type": "StaticSolver" },
+    "newton":       { "type": "NewtonRaphsonSolver", "maxNbIterationsNewton": 1,
+                      "relativeInitialStoppingThreshold": 1e-6,
+                      "relativeSuccessiveStoppingThreshold": 0,
+                      "absoluteResidualStoppingThreshold": 0,
+                      "printLog": false },
+    "linearSolver": {
+      "type": "PCGLinearSolver", "iterations": 1000, "tolerance": 1e-24, "printLog": false,
+      "system": { "type": "PreconditionedMatrixFreeSystem" },
+      "preconditioner": {
+        "type": "SSORPreconditioner",
+        "system": { "type": "MatrixLinearSystem", "template": "CompressedRowSparseMatrixd" }
+      }
+    }
+  },
+  "gridMapping": { "swapping": true },
+  "mesh":        { "cells": [4, 4, 4], "levels": 5, "refinementRatio": 2 },
+  "asymptoticTolerance": 0.10,
+  "expectedOrder": { "L2": 2.0, "H1": 1.0, "EnergyDifference": 1.0 },
+  "expectedOrderTolerance": 0.15,
+  "relativeNoiseFloor": 1e-10
+}

From 98f4c958023bef3364695997ae7203db773c804f Mon Sep 17 00:00:00 2001
From: Themis Skamagkis <themisskamagkis@gmail.com>
Date: Mon, 12 Oct 2026 09:49:19 +0200
Subject: [PATCH 41/55] [verification] Add the
 LinearSmallStrainFEMForceField/Quartic2D_Quad case

---
 .../Quartic2D_Quad.json                       | 31 +++++++++++++++++++
 1 file changed, 31 insertions(+)
 create mode 100644 verification/cases/LinearSmallStrainFEMForceField/Quartic2D_Quad.json

diff --git a/verification/cases/LinearSmallStrainFEMForceField/Quartic2D_Quad.json b/verification/cases/LinearSmallStrainFEMForceField/Quartic2D_Quad.json
new file mode 100644
index 0000000..f9a49bb
--- /dev/null
+++ b/verification/cases/LinearSmallStrainFEMForceField/Quartic2D_Quad.json
@@ -0,0 +1,31 @@
+{
+  "geometry":    { "type": "Tire2D", "inner_radius": 0.25, "outer_radius": 0.5 },
+  "element":     "quad",
+  "material":    { "type": "LinearElastic", "youngModulus": 1.0, "poissonRatio": 0.3 },
+  "forceField":  { "type": "LinearSmallStrainFEMForceField" },
+  "tractionOn":  { "outer": [1, 1] },
+  "solution":    { "type": "Quartic", "amplitude": 1.0 },
+  "errorQuadratureDegree": 3,
+  "sourceQuadratureDegree": 3,
+  "solvers": {
+    "integration":  { "type": "StaticSolver" },
+    "newton":       { "type": "NewtonRaphsonSolver", "maxNbIterationsNewton": 1,
+                      "relativeInitialStoppingThreshold": 1e-6,
+                      "relativeSuccessiveStoppingThreshold": 0,
+                      "absoluteResidualStoppingThreshold": 0,
+                      "printLog": false },
+    "linearSolver": {
+      "type": "PCGLinearSolver", "iterations": 1000, "tolerance": 1e-24, "printLog": false,
+      "system": { "type": "PreconditionedMatrixFreeSystem" },
+      "preconditioner": {
+        "type": "SSORPreconditioner",
+        "system": { "type": "MatrixLinearSystem", "template": "CompressedRowSparseMatrixd" }
+      }
+    }
+  },
+  "mesh":        { "cells": [2, 16], "levels": 5, "refinementRatio": 2 },
+  "asymptoticTolerance": 0.10,
+  "expectedOrder": { "L2": 2.0, "H1": 1.0, "EnergyDifference": 1.0 },
+  "expectedOrderTolerance": 0.15,
+  "relativeNoiseFloor": 1e-10
+}

From 408fffb64a187ebfcc121ac61ed3a6f188d35eb8 Mon Sep 17 00:00:00 2001
From: Themis Skamagkis <themisskamagkis@gmail.com>
Date: Mon, 12 Oct 2026 09:49:19 +0200
Subject: [PATCH 42/55] [verification] Add the
 LinearSmallStrainFEMForceField/Quartic2D_Tri case

---
 .../Quartic2D_Tri.json                        | 31 +++++++++++++++++++
 1 file changed, 31 insertions(+)
 create mode 100644 verification/cases/LinearSmallStrainFEMForceField/Quartic2D_Tri.json

diff --git a/verification/cases/LinearSmallStrainFEMForceField/Quartic2D_Tri.json b/verification/cases/LinearSmallStrainFEMForceField/Quartic2D_Tri.json
new file mode 100644
index 0000000..c161f73
--- /dev/null
+++ b/verification/cases/LinearSmallStrainFEMForceField/Quartic2D_Tri.json
@@ -0,0 +1,31 @@
+{
+  "geometry":    { "type": "Tire2D", "inner_radius": 0.25, "outer_radius": 0.5 },
+  "element":     "tri",
+  "material":    { "type": "LinearElastic", "youngModulus": 1.0, "poissonRatio": 0.3 },
+  "forceField":  { "type": "LinearSmallStrainFEMForceField" },
+  "tractionOn":  { "outer": [1, 1] },
+  "solution":    { "type": "Quartic", "amplitude": 1.0 },
+  "errorQuadratureDegree": 2,
+  "sourceQuadratureDegree": 2,
+  "solvers": {
+    "integration":  { "type": "StaticSolver" },
+    "newton":       { "type": "NewtonRaphsonSolver", "maxNbIterationsNewton": 1,
+                      "relativeInitialStoppingThreshold": 1e-6,
+                      "relativeSuccessiveStoppingThreshold": 0,
+                      "absoluteResidualStoppingThreshold": 0,
+                      "printLog": false },
+    "linearSolver": {
+      "type": "PCGLinearSolver", "iterations": 1000, "tolerance": 1e-24, "printLog": false,
+      "system": { "type": "PreconditionedMatrixFreeSystem" },
+      "preconditioner": {
+        "type": "SSORPreconditioner",
+        "system": { "type": "MatrixLinearSystem", "template": "CompressedRowSparseMatrixd" }
+      }
+    }
+  },
+  "mesh":        { "cells": [2, 16], "levels": 5, "refinementRatio": 2 },
+  "asymptoticTolerance": 0.10,
+  "expectedOrder": { "L2": 2.0, "H1": 1.0, "EnergyDifference": 1.0 },
+  "expectedOrderTolerance": 0.15,
+  "relativeNoiseFloor": 1e-10
+}

From c0ef61db2db61f7a13f08d883c4c7e7eeb7bd8bf Mon Sep 17 00:00:00 2001
From: Themis Skamagkis <themisskamagkis@gmail.com>
Date: Mon, 12 Oct 2026 09:49:19 +0200
Subject: [PATCH 43/55] [verification] Add the
 LinearSmallStrainFEMForceField/Quartic3D_Hexa case

---
 .../Quartic3D_Hexa.json                       | 31 +++++++++++++++++++
 1 file changed, 31 insertions(+)
 create mode 100644 verification/cases/LinearSmallStrainFEMForceField/Quartic3D_Hexa.json

diff --git a/verification/cases/LinearSmallStrainFEMForceField/Quartic3D_Hexa.json b/verification/cases/LinearSmallStrainFEMForceField/Quartic3D_Hexa.json
new file mode 100644
index 0000000..7bdde69
--- /dev/null
+++ b/verification/cases/LinearSmallStrainFEMForceField/Quartic3D_Hexa.json
@@ -0,0 +1,31 @@
+{
+  "geometry":    { "type": "Tire3D", "inner_radius": 0.25, "outer_radius": 0.5, "extrusion_length": 0.25 },
+  "element":     "hexa",
+  "material":    { "type": "LinearElastic", "youngModulus": 1.0, "poissonRatio": 0.3 },
+  "forceField":  { "type": "LinearSmallStrainFEMForceField" },
+  "tractionOn":  { "front": [1, 1, 1], "back": [1, 1, 1] },
+  "solution":    { "type": "Quartic", "amplitude": 1.0 },
+  "errorQuadratureDegree": 3,
+  "sourceQuadratureDegree": 3,
+  "solvers": {
+    "integration":  { "type": "StaticSolver" },
+    "newton":       { "type": "NewtonRaphsonSolver", "maxNbIterationsNewton": 1,
+                      "relativeInitialStoppingThreshold": 1e-6,
+                      "relativeSuccessiveStoppingThreshold": 0,
+                      "absoluteResidualStoppingThreshold": 0,
+                      "printLog": false },
+    "linearSolver": {
+      "type": "PCGLinearSolver", "iterations": 1000, "tolerance": 1e-24, "printLog": false,
+      "system": { "type": "PreconditionedMatrixFreeSystem" },
+      "preconditioner": {
+        "type": "SSORPreconditioner",
+        "system": { "type": "MatrixLinearSystem", "template": "CompressedRowSparseMatrixd" }
+      }
+    }
+  },
+  "mesh":        { "cells": [2, 16, 2], "levels": 4, "refinementRatio": 2 },
+  "asymptoticTolerance": 0.10,
+  "expectedOrder": { "L2": 2.0, "H1": 1.0, "EnergyDifference": 1.0 },
+  "expectedOrderTolerance": 0.15,
+  "relativeNoiseFloor": 1e-10
+}

From d606b5b0175ea16424803e975d9ce022ac8fd35b Mon Sep 17 00:00:00 2001
From: Themis Skamagkis <themisskamagkis@gmail.com>
Date: Mon, 12 Oct 2026 09:49:20 +0200
Subject: [PATCH 44/55] [verification] Add the
 LinearSmallStrainFEMForceField/Quartic3D_Tetra case

---
 .../Quartic3D_Tetra.json                      | 32 +++++++++++++++++++
 1 file changed, 32 insertions(+)
 create mode 100644 verification/cases/LinearSmallStrainFEMForceField/Quartic3D_Tetra.json

diff --git a/verification/cases/LinearSmallStrainFEMForceField/Quartic3D_Tetra.json b/verification/cases/LinearSmallStrainFEMForceField/Quartic3D_Tetra.json
new file mode 100644
index 0000000..d99fd10
--- /dev/null
+++ b/verification/cases/LinearSmallStrainFEMForceField/Quartic3D_Tetra.json
@@ -0,0 +1,32 @@
+{
+  "geometry":    { "type": "Tire3D", "inner_radius": 0.25, "outer_radius": 0.5, "extrusion_length": 0.25 },
+  "element":     "tet",
+  "material":    { "type": "LinearElastic", "youngModulus": 1.0, "poissonRatio": 0.3 },
+  "forceField":  { "type": "LinearSmallStrainFEMForceField" },
+  "tractionOn":  { "front": [1, 1, 1], "back": [1, 1, 1] },
+  "solution":    { "type": "Quartic", "amplitude": 1.0 },
+  "errorQuadratureDegree": 2,
+  "sourceQuadratureDegree": 2,
+  "solvers": {
+    "integration":  { "type": "StaticSolver" },
+    "newton":       { "type": "NewtonRaphsonSolver", "maxNbIterationsNewton": 1,
+                      "relativeInitialStoppingThreshold": 1e-6,
+                      "relativeSuccessiveStoppingThreshold": 0,
+                      "absoluteResidualStoppingThreshold": 0,
+                      "printLog": false },
+    "linearSolver": {
+      "type": "PCGLinearSolver", "iterations": 1000, "tolerance": 1e-24, "printLog": false,
+      "system": { "type": "PreconditionedMatrixFreeSystem" },
+      "preconditioner": {
+        "type": "SSORPreconditioner",
+        "system": { "type": "MatrixLinearSystem", "template": "CompressedRowSparseMatrixd" }
+      }
+    }
+  },
+  "gridMapping": { "swapping": true },
+  "mesh":        { "cells": [2, 16, 2], "levels": 4, "refinementRatio": 2 },
+  "asymptoticTolerance": 0.10,
+  "expectedOrder": { "L2": 2.0, "H1": 1.0, "EnergyDifference": 1.0 },
+  "expectedOrderTolerance": 0.15,
+  "relativeNoiseFloor": 1e-10
+}

From 31524d18621071afed4f925167959779bacb6389 Mon Sep 17 00:00:00 2001
From: Themis Skamagkis <themisskamagkis@gmail.com>
Date: Mon, 12 Oct 2026 09:49:20 +0200
Subject: [PATCH 45/55] [verification] Add the
 LinearSmallStrainFEMForceField/Trigonometric1D case

---
 .../Trigonometric1D.json                      | 24 +++++++++++++++++++
 1 file changed, 24 insertions(+)
 create mode 100644 verification/cases/LinearSmallStrainFEMForceField/Trigonometric1D.json

diff --git a/verification/cases/LinearSmallStrainFEMForceField/Trigonometric1D.json b/verification/cases/LinearSmallStrainFEMForceField/Trigonometric1D.json
new file mode 100644
index 0000000..827a2f2
--- /dev/null
+++ b/verification/cases/LinearSmallStrainFEMForceField/Trigonometric1D.json
@@ -0,0 +1,24 @@
+{
+  "geometry":    { "type": "Bar1D", "length": 2.0 },
+  "element":     "edge",
+  "material":    { "type": "LinearElastic", "youngModulus": 1e6, "poissonRatio": 0.3 },
+  "forceField":  { "type": "LinearSmallStrainFEMForceField" },
+  "tractionOn":  { "right": [1] },
+  "solution":    { "type": "Trigonometric", "amplitude": 1 },
+  "errorQuadratureDegree": 3,
+  "sourceQuadratureDegree": 2,
+  "solvers": {
+    "integration":  { "type": "StaticSolver" },
+    "newton":       { "type": "NewtonRaphsonSolver", "maxNbIterationsNewton": 1,
+                      "relativeInitialStoppingThreshold": 1e-6,
+                      "relativeSuccessiveStoppingThreshold": 0,
+                      "absoluteResidualStoppingThreshold": 0,
+                      "printLog": false },
+    "linearSolver": { "type": "SparseLDLSolver", "template": "CompressedRowSparseMatrixd" }
+  },
+  "mesh":        { "cells": [3], "levels": 7, "refinementRatio": 2 },
+  "asymptoticTolerance": 0.10,
+  "expectedOrder": { "L2": 2.0, "H1": 1.0, "EnergyDifference": 1.0 },
+  "expectedOrderTolerance": 0.15,
+  "relativeNoiseFloor": 1e-12
+}

From 8523c5585cf74ca20123f7a431938dfc5a6bd329 Mon Sep 17 00:00:00 2001
From: Themis Skamagkis <themisskamagkis@gmail.com>
Date: Mon, 12 Oct 2026 09:49:20 +0200
Subject: [PATCH 46/55] [verification] Add the
 LinearSmallStrainFEMForceField/Trigonometric2D_Quad case

---
 .../Trigonometric2D_Quad.json                 | 24 +++++++++++++++++++
 1 file changed, 24 insertions(+)
 create mode 100644 verification/cases/LinearSmallStrainFEMForceField/Trigonometric2D_Quad.json

diff --git a/verification/cases/LinearSmallStrainFEMForceField/Trigonometric2D_Quad.json b/verification/cases/LinearSmallStrainFEMForceField/Trigonometric2D_Quad.json
new file mode 100644
index 0000000..f608cf9
--- /dev/null
+++ b/verification/cases/LinearSmallStrainFEMForceField/Trigonometric2D_Quad.json
@@ -0,0 +1,24 @@
+{
+  "geometry":    { "type": "Beam2D", "length": 1.0, "width": 1.0 },
+  "element":     "quad",
+  "material":    { "type": "LinearElastic", "youngModulus": 1e6, "poissonRatio": 0.3 },
+  "forceField":  { "type": "LinearSmallStrainFEMForceField" },
+  "tractionOn":  { "left": [0, 1], "bottom": [1, 0], "right": [1, 1], "top": [1, 1] },
+  "solution":    { "type": "Trigonometric", "amplitude": 0.1 },
+  "errorQuadratureDegree": 3,
+  "sourceQuadratureDegree": 3,
+  "solvers": {
+    "integration":  { "type": "StaticSolver" },
+    "newton":       { "type": "NewtonRaphsonSolver", "maxNbIterationsNewton": 1,
+                      "relativeInitialStoppingThreshold": 1e-6,
+                      "relativeSuccessiveStoppingThreshold": 0,
+                      "absoluteResidualStoppingThreshold": 0,
+                      "printLog": false },
+    "linearSolver": { "type": "SparseLDLSolver", "template": "CompressedRowSparseMatrixd" }
+  },
+  "mesh":        { "cells": [3, 3], "levels": 6, "refinementRatio": 2 },
+  "asymptoticTolerance": 0.10,
+  "expectedOrder": { "L2": 2.0, "H1": 1.0, "EnergyDifference": 1.0 },
+  "expectedOrderTolerance": 0.15,
+  "relativeNoiseFloor": 1e-12
+}

From 50c3f62bcfc7e2236328aa188945de1604ddc74b Mon Sep 17 00:00:00 2001
From: Themis Skamagkis <themisskamagkis@gmail.com>
Date: Mon, 12 Oct 2026 09:49:20 +0200
Subject: [PATCH 47/55] [verification] Add the
 LinearSmallStrainFEMForceField/Trigonometric2D_Quad_Dirichlet case

---
 .../Trigonometric2D_Quad_Dirichlet.json       | 23 +++++++++++++++++++
 1 file changed, 23 insertions(+)
 create mode 100644 verification/cases/LinearSmallStrainFEMForceField/Trigonometric2D_Quad_Dirichlet.json

diff --git a/verification/cases/LinearSmallStrainFEMForceField/Trigonometric2D_Quad_Dirichlet.json b/verification/cases/LinearSmallStrainFEMForceField/Trigonometric2D_Quad_Dirichlet.json
new file mode 100644
index 0000000..a57bd17
--- /dev/null
+++ b/verification/cases/LinearSmallStrainFEMForceField/Trigonometric2D_Quad_Dirichlet.json
@@ -0,0 +1,23 @@
+{
+  "geometry":    { "type": "Beam2D", "length": 1.0, "width": 1.0 },
+  "element":     "quad",
+  "material":    { "type": "LinearElastic", "youngModulus": 1e6, "poissonRatio": 0.3 },
+  "forceField":  { "type": "LinearSmallStrainFEMForceField" },
+  "solution":    { "type": "Trigonometric", "amplitude": 0.02 },
+  "errorQuadratureDegree": 3,
+  "sourceQuadratureDegree": 3,
+  "solvers": {
+    "integration":  { "type": "StaticSolver" },
+    "newton":       { "type": "NewtonRaphsonSolver", "maxNbIterationsNewton": 1,
+                      "relativeInitialStoppingThreshold": 1e-6,
+                      "relativeSuccessiveStoppingThreshold": 0,
+                      "absoluteResidualStoppingThreshold": 0,
+                      "printLog": false },
+    "linearSolver": { "type": "SparseLDLSolver", "template": "CompressedRowSparseMatrixd" }
+  },
+  "mesh":        { "cells": [3, 3], "levels": 5, "refinementRatio": 2 },
+  "asymptoticTolerance": 0.10,
+  "expectedOrder": { "L2": 2.0, "H1": 1.0, "EnergyDifference": 1.0 },
+  "expectedOrderTolerance": 0.15,
+  "relativeNoiseFloor": 1e-12
+}

From 4aef8ea514651bf3b4a8c85662956e8514299e1b Mon Sep 17 00:00:00 2001
From: Themis Skamagkis <themisskamagkis@gmail.com>
Date: Mon, 12 Oct 2026 09:49:20 +0200
Subject: [PATCH 48/55] [verification] Add the
 LinearSmallStrainFEMForceField/Trigonometric2D_Tri case

---
 .../Trigonometric2D_Tri.json                  | 24 +++++++++++++++++++
 1 file changed, 24 insertions(+)
 create mode 100644 verification/cases/LinearSmallStrainFEMForceField/Trigonometric2D_Tri.json

diff --git a/verification/cases/LinearSmallStrainFEMForceField/Trigonometric2D_Tri.json b/verification/cases/LinearSmallStrainFEMForceField/Trigonometric2D_Tri.json
new file mode 100644
index 0000000..5024659
--- /dev/null
+++ b/verification/cases/LinearSmallStrainFEMForceField/Trigonometric2D_Tri.json
@@ -0,0 +1,24 @@
+{
+  "geometry":    { "type": "Beam2D", "length": 1.0, "width": 1.0 },
+  "element":     "tri",
+  "material":    { "type": "LinearElastic", "youngModulus": 1e6, "poissonRatio": 0.3 },
+  "forceField":  { "type": "LinearSmallStrainFEMForceField" },
+  "tractionOn":  { "left": [0, 1], "bottom": [1, 0], "right": [1, 1], "top": [1, 1] },
+  "solution":    { "type": "Trigonometric", "amplitude": 0.1 },
+  "errorQuadratureDegree": 2,
+  "sourceQuadratureDegree": 2,
+  "solvers": {
+    "integration":  { "type": "StaticSolver" },
+    "newton":       { "type": "NewtonRaphsonSolver", "maxNbIterationsNewton": 1,
+                      "relativeInitialStoppingThreshold": 1e-6,
+                      "relativeSuccessiveStoppingThreshold": 0,
+                      "absoluteResidualStoppingThreshold": 0,
+                      "printLog": false },
+    "linearSolver": { "type": "SparseLDLSolver", "template": "CompressedRowSparseMatrixd" }
+  },
+  "mesh":        { "cells": [3, 3], "levels": 6, "refinementRatio": 2 },
+  "asymptoticTolerance": 0.10,
+  "expectedOrder": { "L2": 2.0, "H1": 1.0, "EnergyDifference": 1.0 },
+  "expectedOrderTolerance": 0.15,
+  "relativeNoiseFloor": 1e-12
+}

From e515655f272eb2b34b520cf4a5f09c096a5dbd60 Mon Sep 17 00:00:00 2001
From: Themis Skamagkis <themisskamagkis@gmail.com>
Date: Mon, 12 Oct 2026 09:49:20 +0200
Subject: [PATCH 49/55] [verification] Add the
 LinearSmallStrainFEMForceField/Trigonometric2D_Tri_Dirichlet case

---
 .../Trigonometric2D_Tri_Dirichlet.json        | 23 +++++++++++++++++++
 1 file changed, 23 insertions(+)
 create mode 100644 verification/cases/LinearSmallStrainFEMForceField/Trigonometric2D_Tri_Dirichlet.json

diff --git a/verification/cases/LinearSmallStrainFEMForceField/Trigonometric2D_Tri_Dirichlet.json b/verification/cases/LinearSmallStrainFEMForceField/Trigonometric2D_Tri_Dirichlet.json
new file mode 100644
index 0000000..1efc304
--- /dev/null
+++ b/verification/cases/LinearSmallStrainFEMForceField/Trigonometric2D_Tri_Dirichlet.json
@@ -0,0 +1,23 @@
+{
+  "geometry":    { "type": "Beam2D", "length": 1.0, "width": 1.0 },
+  "element":     "tri",
+  "material":    { "type": "LinearElastic", "youngModulus": 1e6, "poissonRatio": 0.3 },
+  "forceField":  { "type": "LinearSmallStrainFEMForceField" },
+  "solution":    { "type": "Trigonometric", "amplitude": 0.02 },
+  "errorQuadratureDegree": 2,
+  "sourceQuadratureDegree": 2,
+  "solvers": {
+    "integration":  { "type": "StaticSolver" },
+    "newton":       { "type": "NewtonRaphsonSolver", "maxNbIterationsNewton": 1,
+                      "relativeInitialStoppingThreshold": 1e-6,
+                      "relativeSuccessiveStoppingThreshold": 0,
+                      "absoluteResidualStoppingThreshold": 0,
+                      "printLog": false },
+    "linearSolver": { "type": "SparseLDLSolver", "template": "CompressedRowSparseMatrixd" }
+  },
+  "mesh":        { "cells": [3, 3], "levels": 5, "refinementRatio": 2 },
+  "asymptoticTolerance": 0.10,
+  "expectedOrder": { "L2": 2.0, "H1": 1.0, "EnergyDifference": 1.0 },
+  "expectedOrderTolerance": 0.15,
+  "relativeNoiseFloor": 1e-12
+}

From f22a739409e4a8e781841d254a4e367efbf52589 Mon Sep 17 00:00:00 2001
From: Themis Skamagkis <themisskamagkis@gmail.com>
Date: Mon, 12 Oct 2026 09:49:20 +0200
Subject: [PATCH 50/55] [verification] Add the
 LinearSmallStrainFEMForceField/Trigonometric2D_Tri_PlaneStrain case

---
 .../Trigonometric2D_Tri_PlaneStrain.json      | 23 +++++++++++++++++++
 1 file changed, 23 insertions(+)
 create mode 100644 verification/cases/LinearSmallStrainFEMForceField/Trigonometric2D_Tri_PlaneStrain.json

diff --git a/verification/cases/LinearSmallStrainFEMForceField/Trigonometric2D_Tri_PlaneStrain.json b/verification/cases/LinearSmallStrainFEMForceField/Trigonometric2D_Tri_PlaneStrain.json
new file mode 100644
index 0000000..c32ddef
--- /dev/null
+++ b/verification/cases/LinearSmallStrainFEMForceField/Trigonometric2D_Tri_PlaneStrain.json
@@ -0,0 +1,23 @@
+{
+  "geometry":    { "type": "Beam2D", "length": 1.0, "width": 1.0, "spatialDimensions": 3 },
+  "element":     "tri",
+  "material":    { "type": "LinearElastic", "youngModulus": 1e6, "poissonRatio": 0.3 },
+  "forceField":  { "type": "LinearSmallStrainFEMForceField" },
+  "solution":    { "type": "Trigonometric", "amplitude": 0.1 },
+  "errorQuadratureDegree": 2,
+  "sourceQuadratureDegree": 2,
+  "solvers": {
+    "integration":  { "type": "StaticSolver" },
+    "newton":       { "type": "NewtonRaphsonSolver", "maxNbIterationsNewton": 1,
+                      "relativeInitialStoppingThreshold": 1e-6,
+                      "relativeSuccessiveStoppingThreshold": 0,
+                      "absoluteResidualStoppingThreshold": 0,
+                      "printLog": false },
+    "linearSolver": { "type": "SparseLDLSolver", "template": "CompressedRowSparseMatrixd" }
+  },
+  "mesh":        { "cells": [3, 3], "levels": 6, "refinementRatio": 2 },
+  "asymptoticTolerance": 0.10,
+  "expectedOrder": { "L2": 2.0, "H1": 1.0, "EnergyDifference": 1.0 },
+  "expectedOrderTolerance": 0.15,
+  "relativeNoiseFloor": 1e-12
+}

From ddcd1f3f3b3ca697b4534e31feea8426f74bb426 Mon Sep 17 00:00:00 2001
From: Themis Skamagkis <themisskamagkis@gmail.com>
Date: Mon, 12 Oct 2026 09:49:20 +0200
Subject: [PATCH 51/55] [verification] Add the
 LinearSmallStrainFEMForceField/Trigonometric3D_Hexa case

---
 .../Trigonometric3D_Hexa.json                 | 31 +++++++++++++++++++
 1 file changed, 31 insertions(+)
 create mode 100644 verification/cases/LinearSmallStrainFEMForceField/Trigonometric3D_Hexa.json

diff --git a/verification/cases/LinearSmallStrainFEMForceField/Trigonometric3D_Hexa.json b/verification/cases/LinearSmallStrainFEMForceField/Trigonometric3D_Hexa.json
new file mode 100644
index 0000000..2e0c27f
--- /dev/null
+++ b/verification/cases/LinearSmallStrainFEMForceField/Trigonometric3D_Hexa.json
@@ -0,0 +1,31 @@
+{
+  "geometry":    { "type": "Beam3D", "length": 1.0, "width": 1.0, "height": 1.0 },
+  "element":     "hexa",
+  "material":    { "type": "LinearElastic", "youngModulus": 1e6, "poissonRatio": 0.3 },
+  "forceField":  { "type": "LinearSmallStrainFEMForceField" },
+  "tractionOn":  { "left": [0, 1, 1], "bottom": [1, 0, 1], "front": [1, 1, 0], "right": [1, 1, 1], "top": [1, 1, 1], "back": [1, 1, 1] },
+  "solution":    { "type": "Trigonometric", "amplitude": 0.1 },
+  "errorQuadratureDegree": 3,
+  "sourceQuadratureDegree": 3,
+  "solvers": {
+    "integration":  { "type": "StaticSolver" },
+    "newton":       { "type": "NewtonRaphsonSolver", "maxNbIterationsNewton": 1,
+                      "relativeInitialStoppingThreshold": 1e-6,
+                      "relativeSuccessiveStoppingThreshold": 0,
+                      "absoluteResidualStoppingThreshold": 0,
+                      "printLog": false },
+    "linearSolver": {
+      "type": "PCGLinearSolver", "iterations": 1000, "tolerance": 1e-24, "printLog": false,
+      "system": { "type": "PreconditionedMatrixFreeSystem" },
+      "preconditioner": {
+        "type": "SSORPreconditioner",
+        "system": { "type": "MatrixLinearSystem", "template": "CompressedRowSparseMatrixd" }
+      }
+    }
+  },
+  "mesh":        { "cells": [3, 3, 3], "levels": 5, "refinementRatio": 2 },
+  "asymptoticTolerance": 0.10,
+  "expectedOrder": { "L2": 2.0, "H1": 1.0, "EnergyDifference": 1.0 },
+  "expectedOrderTolerance": 0.15,
+  "relativeNoiseFloor": 1e-10
+}

From db61c4fc72f72710594ee64cf19c661b524d6cca Mon Sep 17 00:00:00 2001
From: Themis Skamagkis <themisskamagkis@gmail.com>
Date: Mon, 12 Oct 2026 09:49:20 +0200
Subject: [PATCH 52/55] [verification] Add the
 LinearSmallStrainFEMForceField/Trigonometric3D_Hexa_Dirichlet case

---
 .../Trigonometric3D_Hexa_Dirichlet.json       | 30 +++++++++++++++++++
 1 file changed, 30 insertions(+)
 create mode 100644 verification/cases/LinearSmallStrainFEMForceField/Trigonometric3D_Hexa_Dirichlet.json

diff --git a/verification/cases/LinearSmallStrainFEMForceField/Trigonometric3D_Hexa_Dirichlet.json b/verification/cases/LinearSmallStrainFEMForceField/Trigonometric3D_Hexa_Dirichlet.json
new file mode 100644
index 0000000..f5c0ff8
--- /dev/null
+++ b/verification/cases/LinearSmallStrainFEMForceField/Trigonometric3D_Hexa_Dirichlet.json
@@ -0,0 +1,30 @@
+{
+  "geometry":    { "type": "Beam3D", "length": 1.0, "width": 1.0, "height": 1.0 },
+  "element":     "hexa",
+  "material":    { "type": "LinearElastic", "youngModulus": 1e6, "poissonRatio": 0.3 },
+  "forceField":  { "type": "LinearSmallStrainFEMForceField" },
+  "solution":    { "type": "Trigonometric", "amplitude": 0.02 },
+  "errorQuadratureDegree": 3,
+  "sourceQuadratureDegree": 3,
+  "solvers": {
+    "integration":  { "type": "StaticSolver" },
+    "newton":       { "type": "NewtonRaphsonSolver", "maxNbIterationsNewton": 1,
+                      "relativeInitialStoppingThreshold": 1e-6,
+                      "relativeSuccessiveStoppingThreshold": 0,
+                      "absoluteResidualStoppingThreshold": 0,
+                      "printLog": false },
+    "linearSolver": {
+      "type": "PCGLinearSolver", "iterations": 1000, "tolerance": 1e-24, "printLog": false,
+      "system": { "type": "PreconditionedMatrixFreeSystem" },
+      "preconditioner": {
+        "type": "SSORPreconditioner",
+        "system": { "type": "MatrixLinearSystem", "template": "CompressedRowSparseMatrixd" }
+      }
+    }
+  },
+  "mesh":        { "cells": [3, 3, 3], "levels": 5, "refinementRatio": 2 },
+  "asymptoticTolerance": 0.10,
+  "expectedOrder": { "L2": 2.0, "H1": 1.0, "EnergyDifference": 1.0 },
+  "expectedOrderTolerance": 0.15,
+  "relativeNoiseFloor": 1e-10
+}

From 68c0465804019c37c91a7c64613044f15da53ac2 Mon Sep 17 00:00:00 2001
From: Themis Skamagkis <themisskamagkis@gmail.com>
Date: Mon, 12 Oct 2026 09:49:20 +0200
Subject: [PATCH 53/55] [verification] Add the
 LinearSmallStrainFEMForceField/Trigonometric3D_Tetra case

---
 .../Trigonometric3D_Tetra.json                | 32 +++++++++++++++++++
 1 file changed, 32 insertions(+)
 create mode 100644 verification/cases/LinearSmallStrainFEMForceField/Trigonometric3D_Tetra.json

diff --git a/verification/cases/LinearSmallStrainFEMForceField/Trigonometric3D_Tetra.json b/verification/cases/LinearSmallStrainFEMForceField/Trigonometric3D_Tetra.json
new file mode 100644
index 0000000..26c5075
--- /dev/null
+++ b/verification/cases/LinearSmallStrainFEMForceField/Trigonometric3D_Tetra.json
@@ -0,0 +1,32 @@
+{
+  "geometry":    { "type": "Beam3D", "length": 1.0, "width": 1.0, "height": 1.0 },
+  "element":     "tet",
+  "material":    { "type": "LinearElastic", "youngModulus": 1e6, "poissonRatio": 0.3 },
+  "forceField":  { "type": "LinearSmallStrainFEMForceField" },
+  "tractionOn":  { "left": [0, 1, 1], "bottom": [1, 0, 1], "front": [1, 1, 0], "right": [1, 1, 1], "top": [1, 1, 1], "back": [1, 1, 1] },
+  "solution":    { "type": "Trigonometric", "amplitude": 0.1 },
+  "errorQuadratureDegree": 2,
+  "sourceQuadratureDegree": 2,
+  "solvers": {
+    "integration":  { "type": "StaticSolver" },
+    "newton":       { "type": "NewtonRaphsonSolver", "maxNbIterationsNewton": 1,
+                      "relativeInitialStoppingThreshold": 1e-6,
+                      "relativeSuccessiveStoppingThreshold": 0,
+                      "absoluteResidualStoppingThreshold": 0,
+                      "printLog": false },
+    "linearSolver": {
+      "type": "PCGLinearSolver", "iterations": 1000, "tolerance": 1e-24, "printLog": false,
+      "system": { "type": "PreconditionedMatrixFreeSystem" },
+      "preconditioner": {
+        "type": "SSORPreconditioner",
+        "system": { "type": "MatrixLinearSystem", "template": "CompressedRowSparseMatrixd" }
+      }
+    }
+  },
+  "gridMapping": { "swapping": true },
+  "mesh":        { "cells": [4, 4, 4], "levels": 5, "refinementRatio": 2 },
+  "asymptoticTolerance": 0.10,
+  "expectedOrder": { "L2": 2.0, "H1": 1.0, "EnergyDifference": 1.0 },
+  "expectedOrderTolerance": 0.15,
+  "relativeNoiseFloor": 1e-10
+}

From 37756bb89ec911f43c433700f897580ceaad6ff0 Mon Sep 17 00:00:00 2001
From: Themis Skamagkis <themisskamagkis@gmail.com>
Date: Mon, 12 Oct 2026 09:49:20 +0200
Subject: [PATCH 54/55] [verification] Add the
 LinearSmallStrainFEMForceField/Trigonometric3D_Tetra_Dirichlet case

---
 .../Trigonometric3D_Tetra_Dirichlet.json      | 31 +++++++++++++++++++
 1 file changed, 31 insertions(+)
 create mode 100644 verification/cases/LinearSmallStrainFEMForceField/Trigonometric3D_Tetra_Dirichlet.json

diff --git a/verification/cases/LinearSmallStrainFEMForceField/Trigonometric3D_Tetra_Dirichlet.json b/verification/cases/LinearSmallStrainFEMForceField/Trigonometric3D_Tetra_Dirichlet.json
new file mode 100644
index 0000000..bb6b7ba
--- /dev/null
+++ b/verification/cases/LinearSmallStrainFEMForceField/Trigonometric3D_Tetra_Dirichlet.json
@@ -0,0 +1,31 @@
+{
+  "geometry":    { "type": "Beam3D", "length": 1.0, "width": 1.0, "height": 1.0 },
+  "element":     "tet",
+  "material":    { "type": "LinearElastic", "youngModulus": 1e6, "poissonRatio": 0.3 },
+  "forceField":  { "type": "LinearSmallStrainFEMForceField" },
+  "solution":    { "type": "Trigonometric", "amplitude": 0.02 },
+  "errorQuadratureDegree": 2,
+  "sourceQuadratureDegree": 2,
+  "solvers": {
+    "integration":  { "type": "StaticSolver" },
+    "newton":       { "type": "NewtonRaphsonSolver", "maxNbIterationsNewton": 1,
+                      "relativeInitialStoppingThreshold": 1e-6,
+                      "relativeSuccessiveStoppingThreshold": 0,
+                      "absoluteResidualStoppingThreshold": 0,
+                      "printLog": false },
+    "linearSolver": {
+      "type": "PCGLinearSolver", "iterations": 1000, "tolerance": 1e-24, "printLog": false,
+      "system": { "type": "PreconditionedMatrixFreeSystem" },
+      "preconditioner": {
+        "type": "SSORPreconditioner",
+        "system": { "type": "MatrixLinearSystem", "template": "CompressedRowSparseMatrixd" }
+      }
+    }
+  },
+  "gridMapping": { "swapping": true },
+  "mesh":        { "cells": [4, 4, 4], "levels": 5, "refinementRatio": 2 },
+  "asymptoticTolerance": 0.10,
+  "expectedOrder": { "L2": 2.0, "H1": 1.0, "EnergyDifference": 1.0 },
+  "expectedOrderTolerance": 0.15,
+  "relativeNoiseFloor": 1e-10
+}

From 587a38902ccb824da26b2928ea4bbaf5078b21d4 Mon Sep 17 00:00:00 2001
From: Themis Skamagkis <themisskamagkis@gmail.com>
Date: Mon, 12 Oct 2026 14:35:15 +0200
Subject: [PATCH 55/55] [docs] Update README

---
 README.md | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/README.md b/README.md
index 17d619f..497734e 100644
--- a/README.md
+++ b/README.md
@@ -10,6 +10,9 @@ The error must shrink at the rate the discretisation predicts. Covered force fie
 
 To get started, read [`verification/docs/IntroductionMMS.ipynb`](verification/docs/IntroductionMMS.ipynb): it presents the method and builds one verification study step by step, from the manufactured solution to the input file the suite runs.
 
+The following force field components are currently verified within this repository:
+- `LinearSmallStrainFEMForceField` [`verification/docs/LinearSmallStrainFEMForceField.ipynb`](verification/docs/LinearSmallStrainFEMForceField.ipynb)
+
 ## Layout
 
 ```