diff --git a/examples/Freefem/MMS/1D/bar.py b/examples/Freefem/MMS/1D/bar.py
index 3904cb11..6c086be1 100644
--- a/examples/Freefem/MMS/1D/bar.py
+++ b/examples/Freefem/MMS/1D/bar.py
@@ -53,16 +53,20 @@ def compute(nodes, topology):
     return compute
 
 
-def build_bar_scene(root, mms, E_eff, nx):
+def build_bar_scene(root, mms, E_eff, nx, force_field, linear_solver):
     """Populate root with a static 1D bar scene on the non-dimensional domain [0,1].
 
     BodyForce is assembled after init by BodyForceAssembler. BCs: Dirichlet at
     x=0, Neumann `mms.traction_bc(E_eff)` at x=1.
+
+    force_field   : name of the FEM force field to test
+    linear_solver : dict {"type": <name>, "parameters": {...}} of the linear solver
     """
     root.addObject('RequiredPlugin', pluginName=[
         "Elasticity",
         "Sofa.Component.Constraint.Projective",
         "Sofa.Component.LinearSolver.Direct",
+        "Sofa.Component.LinearSolver.Iterative",
         "Sofa.Component.MechanicalLoad",
         "Sofa.Component.ODESolver.Backward",
         "Sofa.Component.StateContainer",
@@ -87,9 +91,9 @@ def build_bar_scene(root, mms, E_eff, nx):
                   maxNbIterationsNewton=10,
                   absoluteResidualStoppingThreshold=1e-10,
                   printLog=False)
-    Bar.addObject('SparseLDLSolver',
+    Bar.addObject(linear_solver["type"],
                   name="linearSolver",
-                  template="CompressedRowSparseMatrixd")
+                  **linear_solver["parameters"])
     Bar.addObject('StaticSolver',
                   name="staticSolver",
                   newtonSolver="@newtonSolver",
@@ -103,7 +107,7 @@ def build_bar_scene(root, mms, E_eff, nx):
                          name="dofs",
                          template="Vec1d")
 
-    Bar.addObject('LinearSmallStrainFEMForceField',
+    Bar.addObject(force_field,
                   name="FEM",
                   template="Vec1d",
                   youngModulus=E_eff,
@@ -130,15 +134,18 @@ def case_scene(mms):
     """Return a `createScene(rootNode)` bound to this MMS case."""
     def createScene(rootNode):
         cfg = load_params()
-        build_bar_scene(rootNode, mms, cfg["E_eff"], cfg["reference"]["nx"])
+        build_bar_scene(rootNode, mms, cfg["E_eff"], cfg["reference"]["nx"],
+                        force_field=cfg["forceField"],
+                        linear_solver=cfg["linearSolver"])
         return rootNode
     return createScene
 
 
-def solve_bar(mms, E_eff, nx):
+def solve_bar(mms, E_eff, nx, force_field, linear_solver):
     """Build, run one static step, and return a BarSolution1D snapshot."""
     root = Sofa.Core.Node("root")
-    build_bar_scene(root, mms, E_eff, nx)
+    build_bar_scene(root, mms, E_eff, nx, force_field=force_field,
+                    linear_solver=linear_solver)
     Sofa.Simulation.init(root)
     Sofa.Simulation.animate(root, root.dt.value)
     Bar   = root.Bar
@@ -180,7 +187,9 @@ def plot_solution(case, x, u_h, u_ex, label_ex):
 def run_reference_scene(mms):
     """Solve one MMS case at the reference mesh, write the solution table and plot."""
     cfg = load_params()
-    sol = solve_bar(mms, cfg["E_eff"], cfg["reference"]["nx"])
+    sol = solve_bar(mms, cfg["E_eff"], cfg["reference"]["nx"],
+                    force_field=cfg["forceField"],
+                    linear_solver=cfg["linearSolver"])
     l2  = l2_error_1d(sol.x0, sol.edges, sol.u_h, mms.u_ex, L2_QUADRATURE_1D)
     h1  = h1_semi_error_1d(sol.x0, sol.edges, sol.u_h, mms.du_ex, H1_QUADRATURE_1D)
     write_solution_table(f"solution_{mms.name}", sol.x0, sol.u_h, mms.u_ex,
diff --git a/examples/Freefem/MMS/1D/params.json b/examples/Freefem/MMS/1D/params.json
index c094267f..166f8e79 100644
--- a/examples/Freefem/MMS/1D/params.json
+++ b/examples/Freefem/MMS/1D/params.json
@@ -1,6 +1,13 @@
 {
     "length": 2.0,
     "youngModulus": 1e6,
+    "forceField": "LinearSmallStrainFEMForceField",
+    "linearSolver": {
+        "type": "SparseLDLSolver",
+        "parameters": {
+            "template": "CompressedRowSparseMatrixd"
+        }
+    },
     "reference": {
         "nx": 10
     },
diff --git a/examples/Freefem/MMS/1D/run_convergence.py b/examples/Freefem/MMS/1D/run_convergence.py
index d063d38d..81ba6a59 100644
--- a/examples/Freefem/MMS/1D/run_convergence.py
+++ b/examples/Freefem/MMS/1D/run_convergence.py
@@ -24,7 +24,9 @@
         stem = f"{mms.name}_convergence"
         hs, errors = run_convergence_series(
             nx_values  = cfg["convergence"]["nx_values"][mms.name],
-            run_fn     = lambda nx, _m=mms: solve_bar(_m, cfg["E_eff"], nx),
+            run_fn     = lambda nx, _m=mms: solve_bar(
+                _m, cfg["E_eff"], nx, force_field=cfg["forceField"],
+                linear_solver=cfg["linearSolver"]),
             h_fn       = lambda nx: 1.0 / (nx - 1),
             error_fns  = {
                 "L2": lambda sol, _m=mms: l2_error_1d(
diff --git a/examples/Freefem/MMS/2D/beam.py b/examples/Freefem/MMS/2D/beam.py
index 420ccee4..d483b39a 100644
--- a/examples/Freefem/MMS/2D/beam.py
+++ b/examples/Freefem/MMS/2D/beam.py
@@ -66,12 +66,15 @@ def compute(nodes, topology):
         return F_xy
     return compute
 
-def build_beam_scene(rootNode, mms, element, L=1.0, E=1e6, nu=0.3,
+def build_beam_scene(rootNode, mms, element, force_field, linear_solver,
+                     L=1.0, E=1e6, nu=0.3,
                      nx=10, ny=10, with_visual=True, dim="2d"):
     """Build a SOFA scene for `mms` on the 2D `element` strategy.
 
     dim : "2d" → Vec2d template / plane stress
           "3d" → Vec3d template / plane strain (z coordinate fixed at 0)
+    force_field   : name of the FEM force field to test
+    linear_solver : dict {"type": <name>, "parameters": {...}} of the linear solver
     Returns (dofs, topology). Nodes and connectivity become available
     after `Sofa.Simulation.init(root)` runs, via
     `dofs.rest_position.array()` and `element.read_connectivity(topology)`.
@@ -83,6 +86,7 @@ def build_beam_scene(rootNode, mms, element, L=1.0, E=1e6, nu=0.3,
         "Sofa.Component.Constraint.Projective",
         "Sofa.Component.Engine.Select",
         "Sofa.Component.LinearSolver.Direct",
+        "Sofa.Component.LinearSolver.Iterative",
         "Sofa.Component.MechanicalLoad",
         "Sofa.Component.ODESolver.Backward",
         "Sofa.Component.StateContainer",
@@ -116,8 +120,8 @@ def build_beam_scene(rootNode, mms, element, L=1.0, E=1e6, nu=0.3,
                    maxNbIterationsNewton=1,
                    absoluteResidualStoppingThreshold=1e-10,
                    printLog=False)
-    Beam.addObject("SparseLDLSolver", name="linearSolver",
-                   template="CompressedRowSparseMatrixd")
+    Beam.addObject(linear_solver["type"], name="linearSolver",
+                   **linear_solver["parameters"])
 
     dofs = Beam.addObject("MechanicalObject", name="dofs", template=tmpl,
                           position=nodes_2d.tolist(),
@@ -125,7 +129,7 @@ def build_beam_scene(rootNode, mms, element, L=1.0, E=1e6, nu=0.3,
 
     topology = element.add_topology(Beam)
 
-    Beam.addObject("LinearSmallStrainFEMForceField", name="FEM", template=tmpl,
+    Beam.addObject(force_field, name="FEM", template=tmpl,
                    youngModulus=E, poissonRatio=nu, topology="@topology")
 
     mms.apply_bcs(Beam, nodes_2d, L, dim)
@@ -150,11 +154,12 @@ def build_beam_scene(rootNode, mms, element, L=1.0, E=1e6, nu=0.3,
 # Simulation runner
 # ─────────────────────────────────────────────────────────────────────────────
 
-def solve_beam(elem, mms, L, E, nu, nx, ny, dim="2d"):
+def solve_beam(elem, mms, L, E, nu, nx, ny, force_field, linear_solver, dim="2d"):
     """Build, init, and run one static step. Returns a BeamSolution2D snapshot."""
     root = Sofa.Core.Node("root")
     dofs, topology = build_beam_scene(
-        root, mms, elem, L=L, E=E, nu=nu, nx=nx, ny=ny, with_visual=False, dim=dim
+        root, mms, elem, L=L, E=E, nu=nu, nx=nx, ny=ny, with_visual=False, dim=dim,
+        force_field=force_field, linear_solver=linear_solver
     )
     Sofa.Simulation.init(root)
     # Read topology back from SOFA now that init has populated it.
@@ -246,9 +251,12 @@ def run_reference_scene(elem, mms):
     L, E    = cfg["length"], cfg["youngModulus"]
     nu, dim = ref["nu"], ref["dim"]
     nx = ny = ref["nx"]
+    ff      = cfg["forceField"]
+    ls      = cfg["linearSolver"]
     hyp     = "plane strain" if dim == "3d" else "plane stress"
 
-    sol = solve_beam(elem, mms, L, E, nu, nx, ny, dim=dim)
+    sol = solve_beam(elem, mms, L, E, nu, nx, ny,
+                     force_field=ff, linear_solver=ls, dim=dim)
     l2  = elem.compute_l2(sol, mms, L)
     h1  = elem.compute_h1(sol, mms, L)
 
@@ -281,6 +289,8 @@ def createScene(rootNode):
         build_beam_scene(rootNode, mms, element,
                          L=cfg["length"], E=cfg["youngModulus"],
                          nu=ref["nu"], nx=ref["nx"], ny=ref["nx"],
-                         with_visual=True, dim=ref["dim"])
+                         with_visual=True, dim=ref["dim"],
+                         force_field=cfg["forceField"],
+                         linear_solver=cfg["linearSolver"])
         return rootNode
     return createScene
diff --git a/examples/Freefem/MMS/2D/params.json b/examples/Freefem/MMS/2D/params.json
index f98c79a1..015e4b44 100644
--- a/examples/Freefem/MMS/2D/params.json
+++ b/examples/Freefem/MMS/2D/params.json
@@ -1,8 +1,15 @@
 {
     "length": 1.0,
     "youngModulus": 1e6,
+    "forceField": "LinearSmallStrainFEMForceField",
+    "linearSolver": {
+        "type": "SparseLDLSolver",
+        "parameters": {
+            "template": "CompressedRowSparseMatrixd"
+        }
+    },
     "reference": {
-        "nx": 20,
+        "nx": 5,
         "nu": 0.0,
         "dim": "2d"
     },
diff --git a/examples/Freefem/MMS/2D/run_convergence.py b/examples/Freefem/MMS/2D/run_convergence.py
index bad05b1f..4b458fd8 100644
--- a/examples/Freefem/MMS/2D/run_convergence.py
+++ b/examples/Freefem/MMS/2D/run_convergence.py
@@ -19,7 +19,8 @@
 from output      import plot_convergence
 
 
-def convergence_study(elem_specs, mms, L, E, nu, nx_values, dim="2d"):
+def convergence_study(elem_specs, mms, L, E, nu, nx_values, force_field,
+                      linear_solver, dim="2d"):
     """
     Run a convergence series for each element type in elem_specs, write a
     per-(element) text table, and one shared plot with L²/H¹ for every
@@ -27,6 +28,8 @@ def convergence_study(elem_specs, mms, L, E, nu, nx_values, dim="2d"):
 
     elem_specs : list of dicts with keys 'elem', 'label', 'l2_style', 'h1_style'
     dim        : "2d" (plane stress) or "3d" (plane strain)
+    force_field   : name of the FEM force field to test
+    linear_solver : dict {"type": <name>, "parameters": {...}} of the linear solver
     """
     hyp = "plane strain" if dim == "3d" else "plane stress"
     print(f"\n  PoissonRatio = {nu}  ({hyp})", flush=True)
@@ -40,7 +43,8 @@ def convergence_study(elem_specs, mms, L, E, nu, nx_values, dim="2d"):
         hs, errors = run_convergence_series(
             nx_values  = nx_values,
             run_fn     = lambda nx, _e=elem: solve_beam(
-                _e, mms, L, E, nu, nx, nx, dim=dim),
+                _e, mms, L, E, nu, nx, nx,
+                force_field=force_field, linear_solver=linear_solver, dim=dim),
             h_fn       = lambda nx: L / (nx - 1),
             error_fns  = {
                 "L2": lambda sol, _e=elem: _e.compute_l2(sol, mms, L),
@@ -66,6 +70,8 @@ def convergence_study(elem_specs, mms, L, E, nu, nx_values, dim="2d"):
     cfg  = load_params()
     L    = cfg["length"]
     E    = cfg["youngModulus"]
+    ff   = cfg["forceField"]
+    ls   = cfg["linearSolver"]
     conv = cfg["convergence"]
 
     specs = [
@@ -81,4 +87,5 @@ def convergence_study(elem_specs, mms, L, E, nu, nx_values, dim="2d"):
         print(f"\n== {mms.name} ==")
         for DIM in conv["dim_values"]:
             for nu in conv["nu_values"]:
-                convergence_study(specs, mms, L, E, nu, nx_vals, dim=DIM)
+                convergence_study(specs, mms, L, E, nu, nx_vals,
+                                  force_field=ff, linear_solver=ls, dim=DIM)
diff --git a/examples/Freefem/MMS/2D/test_rigid_rotation.py b/examples/Freefem/MMS/2D/test_rigid_rotation.py
new file mode 100644
index 00000000..0f43b6fb
--- /dev/null
+++ b/examples/Freefem/MMS/2D/test_rigid_rotation.py
@@ -0,0 +1,177 @@
+"""Rigid-rotation test for the corotational FEM force field.
+
+With a pure rigid-body rotation x = R·X the strain is zero.
+so the internal elastic force must disappear. 
+We expect a corotational force field to give zero force but for a linear small-strain field 
+to yield f != 0, because it sees the rotation as strain.
+
+Run:  python test_rigid_rotation.py
+"""
+
+import os
+import sys
+import numpy as np
+import matplotlib.pyplot as plt
+
+import Sofa
+import Sofa.Core
+import Sofa.Simulation
+import SofaRuntime  # noqa: F401
+
+sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
+from beam import load_params, element_quad, RESULTS_DIR
+
+ANGLE_DEG   = 90.0
+TRANSLATION = [0.0, 0.0, 0.0]
+RATIO_TOL   = 1e-6          # corotational |f| must be < RATIO_TOL * linear |f|
+N_STEPS     = 18            # only used by the commented-out solve-based path
+
+
+def Rz(theta):
+    c, s = np.cos(theta), np.sin(theta)
+    return np.array([[c, -s, 0.0],
+                     [s,  c, 0.0],
+                     [0.0, 0.0, 1.0]])
+
+
+def build_scene(root, L, E, nu, nx, force_field):
+    root.addObject("RequiredPlugin", pluginName=[
+        "Elasticity",
+        "Sofa.Component.Mass",
+        "Sofa.Component.ODESolver.Forward",        # EulerExplicitSolver
+        "Sofa.Component.StateContainer",
+        "Sofa.Component.Topology.Container.Grid",
+        "Sofa.Component.Topology.Container.Dynamic",
+        # TODO: Revisit when testing equilibrium
+        # "Sofa.Component.Constraint.Projective",
+        # "Sofa.Component.Engine.Select",
+        # "Sofa.Component.LinearSolver.Iterative",
+        # "Sofa.Component.ODESolver.Backward",
+    ])
+    root.addObject("DefaultAnimationLoop")
+    root.gravity = [0.0, 0.0, 0.0]
+
+    grid = root.addChild("Grid")
+    grid.addObject("RegularGridTopology", name="grid", nx=nx, ny=nx, nz=1,
+                   min=[0.0, 0.0, 0.0], max=[L, L, 0.0])
+
+    with root.addChild("Beam") as beam:
+        # An ODE solver. TODO: Revisit when testing equilibrium
+        beam.addObject("EulerExplicitSolver", name="odeSolver")
+
+        # TODO: An attempt to demonstrate inconsistencies through a solver 
+        #       Cannot work for now. It should be revisited. 
+        # beam.addObject("StaticSolver", name="staticSolver", printLog=False)
+        # beam.addObject("NewtonRaphsonSolver", name="newton",
+        #                maxNbIterationsNewton=25,
+        #                absoluteResidualStoppingThreshold=1e-10, printLog=False)
+        # beam.addObject("CGLinearSolver", name="linearSolver",
+        #                iterations=5000, tolerance=1e-12, threshold=1e-12)
+        # -------------------------------------------------------------------------
+
+        dofs = beam.addObject("MechanicalObject", name="dofs", template="Vec3d",
+                          position="@../Grid/grid.position")
+        beam.addObject("UniformMass", name="mass", totalMass=1.0)
+        element_quad.add_topology(beam)
+        beam.addObject(force_field, name="FEM", template="Vec3d",
+                   youngModulus=E, poissonRatio=nu, topology="@topology")
+
+        # TODO: An attempt to demonstrate inconsistencies through a solver. These are the BCs.
+        #       Cannot work for now. It should be revisited. 
+        # eps = 1e-5
+        # beam.addObject("BoxROI", name="boundary", template="Vec3d", drawBoxes=False,
+        #                box=[[-eps,   -eps,   -1.0, eps,   L + eps, 1.0],   # x = 0
+        #                     [L - eps, -eps,  -1.0, L + eps, L + eps, 1.0],  # x = L
+        #                     [-eps,   -eps,   -1.0, L + eps, eps,   1.0],    # y = 0
+        #                     [-eps,   L - eps, -1.0, L + eps, L + eps, 1.0]])# y = L
+        # beam.addObject("AffineMovementProjectiveConstraint",
+        #                meshIndices=list(range(nx * nx)),
+        #                indices="@boundary.indices",
+        #                beginConstraintTime=0.0, endConstraintTime=float(N_STEPS),
+        #                rotation=Rz(np.radians(ANGLE_DEG)).tolist(),
+        #                translation=TRANSLATION)
+        return dofs
+
+
+def run(force_field):
+    """Impose x = R·X and return (rest X, rotated config, max internal force magnitude)."""
+    cfg = load_params()
+    L, E   = cfg["length"], cfg["youngModulus"]
+    nu, nx = cfg["reference"]["nu"], cfg["reference"]["nx"]
+
+    root = Sofa.Core.Node("root")
+    dofs = build_scene(root, L, E, nu, nx, force_field)
+    Sofa.Simulation.init(root)
+
+    # Take the rest positions and apply x = R @ X
+    # This will apply the rigid-body rotation to all the dofs
+    X  = dofs.rest_position.array().copy()
+    Xr = (Rz(np.radians(ANGLE_DEG)) @ X.T).T + np.array(TRANSLATION)
+    with dofs.position.writeable() as p:
+        p[:] = Xr
+
+    # Animate the scene once. What this does is call addForce, evaluated at x = R·X at the start of the step. 
+    # force is read right after and is the internal force of the pure rotation.
+    # In the case of LinearSmallStrainFEMForceField, we expect to see an emergent force, resulting 
+    # from the force field understanding the rotation as deformation.
+    # In the case of CorotationalFEMForceField, we expect to see zero force, resulting from the 
+    # force field understanding the displacement as a rigid-body rotation causing zero strain.
+    Sofa.Simulation.animate(root, 1e-8)
+    f = dofs.force.array().copy()
+    Sofa.Simulation.unload(root)
+
+    return X, Xr, float(np.max(np.linalg.norm(f, axis=1)))
+
+
+def write_metrics(forces, path):
+    ''' Write the measured forces to a file'''
+    f_cr  = forces["CorotationalFEMForceField"]
+    f_lin = forces["LinearSmallStrainFEMForceField"]
+    with open(path, "w") as f:
+        f.write(f"Rigid-rotation internal force   angle={ANGLE_DEG} deg   "
+                f"translation={TRANSLATION}\n")
+        f.write("=" * 74 + "\n")
+        for ff, fmax in forces.items():
+            f.write(f"  {ff:32s} max |internal force| = {fmax:.6e}\n")
+        ratio = f_cr / f_lin if f_lin != 0 else float("nan")
+        f.write(f"\n  corotational / linear ratio = {ratio:.3e}  "
+                f"(rigid-body property holds if << 1)\n")
+
+
+def plot_forces(X, Xr, forces, path):
+    ''' Plot the forces throughout the 2D plane'''
+    fig, (ax0, ax1) = plt.subplots(1, 2, figsize=(11, 5))
+
+    ax0.scatter(X[:, 0], X[:, 1], s=12, c="lightgray", label="rest X")
+    ax0.scatter(Xr[:, 0], Xr[:, 1], s=12, c="tab:red", label="imposed R·X")
+    ax0.set_aspect("equal"); ax0.set_xlabel("x"); ax0.set_ylabel("y")
+    ax0.set_title(f"imposed rigid rotation {ANGLE_DEG}°")
+    ax0.legend(loc="best"); ax0.grid(True, alpha=0.3)
+
+    names = list(forces.keys())
+    vals  = [max(forces[n], 1e-30) for n in names]     # floor for log scale
+    ax1.bar(range(len(names)), vals, color=["tab:green", "tab:red"])
+    ax1.set_yscale("log")
+    ax1.set_xticks(range(len(names)))
+    ax1.set_xticklabels([n.replace("FEMForceField", "") for n in names])
+    ax1.set_ylabel("max |internal force|")
+    ax1.set_title("internal force under rigid rotation\n(corotational must be ~0)")
+    ax1.grid(True, axis="y", alpha=0.3)
+
+    fig.tight_layout()
+    fig.savefig(path, dpi=150)
+    plt.close(fig)
+
+
+if __name__ == "__main__":
+    forces = {}
+    geom = None
+    for ff in ("CorotationalFEMForceField", "LinearSmallStrainFEMForceField"):
+        X, Xr, fmax = run(ff)
+        forces[ff] = fmax
+        geom = (X, Xr)
+        print(f"{ff:32s} max |internal force| = {fmax:.6e}")
+
+    os.makedirs(RESULTS_DIR, exist_ok=True)
+    write_metrics(forces, os.path.join(RESULTS_DIR, "rigid_rotation_metrics.txt"))
+    plot_forces(*geom, forces, os.path.join(RESULTS_DIR, "rigid_rotation_force.png"))
diff --git a/examples/Freefem/MMS/2D/trigonometric.py b/examples/Freefem/MMS/2D/trigonometric.py
index a79400d8..73412ce9 100644
--- a/examples/Freefem/MMS/2D/trigonometric.py
+++ b/examples/Freefem/MMS/2D/trigonometric.py
@@ -1,8 +1,8 @@
 """
 Trigonometric 2D MMS on [0,L]^2 with linear-elasticity constitutive law:
 
-    u_ex(x, y) = ( sin(pi x / L) cos(pi y / L),
-                   cos(pi x / L) sin(pi y / L) )
+    u_ex(x, y) = A * ( sin(pi x / L) cos(pi y / L),
+                       cos(pi x / L) sin(pi y / L) )
 
     sigma = lambda tr(eps) I + 2 mu eps,
     with (lambda, mu) selected per dim (plane stress vs plane strain).
@@ -16,6 +16,10 @@
                   quad_q1_rule, tri_p1_rule)
 
 
+# Amplitude of the oscillation functions
+AMPLITUDE = 0.1
+
+
 class Trigonometric(MMSCase2D):
     name       = "trigonometric"
     plot_label = (r"$u_x = \sin(\pi x/L)\cos(\pi y/L),\ "
@@ -25,16 +29,18 @@ class Trigonometric(MMSCase2D):
     source_quadrature_tri  = staticmethod(tri_p1_rule(3))
 
     def u_ex(self, x, y, L):
+        A = AMPLITUDE
         k = np.pi / L
-        return (np.sin(k * x) * np.cos(k * y),
-                np.cos(k * x) * np.sin(k * y))
+        return (A * np.sin(k * x) * np.cos(k * y),
+                A * np.cos(k * x) * np.sin(k * y))
 
     def grad_u_ex(self, x, y, L):
+        A = AMPLITUDE
         k = np.pi / L
-        dux_dx =  k * np.cos(k * x) * np.cos(k * y)
-        dux_dy = -k * np.sin(k * x) * np.sin(k * y)
-        duy_dx = -k * np.sin(k * x) * np.sin(k * y)
-        duy_dy =  k * np.cos(k * x) * np.cos(k * y)
+        dux_dx =  A * k * np.cos(k * x) * np.cos(k * y)
+        dux_dy = -A * k * np.sin(k * x) * np.sin(k * y)
+        duy_dx = -A * k * np.sin(k * x) * np.sin(k * y)
+        duy_dy =  A * k * np.cos(k * x) * np.cos(k * y)
         return np.array([[dux_dx, dux_dy],
                          [duy_dx, duy_dy]])
 
@@ -71,6 +77,7 @@ def apply_bcs(self, Beam, nodes_2d, L, dim):
 
     def source(self, x, y, E, nu, L, dim):
         lam, mu = lame(E, nu, dim)
+        A  = AMPLITUDE
         k  = np.pi / L
         ux = np.sin(k * x) * np.cos(k * y)
         uy = np.cos(k * x) * np.sin(k * y)
@@ -84,7 +91,7 @@ def source(self, x, y, E, nu, L, dim):
              + mu * (d2ux_dyy + d2uy_dxy))
         fy = -(mu * (d2ux_dxy + d2uy_dxx) + lam * d2ux_dxy
              + (lam + 2*mu) * d2uy_dyy)
-        return (fx, fy)
+        return (A * fx, A * fy)
 
 
 mms         = Trigonometric()
diff --git a/examples/Freefem/MMS/3D/params.json b/examples/Freefem/MMS/3D/params.json
index 6a902e26..3e3652e6 100644
--- a/examples/Freefem/MMS/3D/params.json
+++ b/examples/Freefem/MMS/3D/params.json
@@ -1,14 +1,21 @@
 {
     "length": 1.0,
     "youngModulus": 1e6,
+    "forceField": "LinearSmallStrainFEMForceField",
+    "linearSolver": {
+        "type": "SparseLDLSolver",
+        "parameters": {
+            "template": "CompressedRowSparseMatrixd"
+        }
+    },
     "reference": {
         "nx": 6,
         "nu": 0.3
     },
     "convergence": {
-        "nu_values": [0.0, 0.3, 0.49],
+        "nu_values": [0.3, 0.49],
         "nx_values": {
-            "sinus_neumann": [4, 6, 8, 12, 16]
+            "sinusoidal": [10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30]
         }
     }
 }
diff --git a/examples/Freefem/MMS/3D/run_convergence.py b/examples/Freefem/MMS/3D/run_convergence.py
index 44591a83..16a1e6ae 100644
--- a/examples/Freefem/MMS/3D/run_convergence.py
+++ b/examples/Freefem/MMS/3D/run_convergence.py
@@ -5,7 +5,7 @@
 plane-stress / plane-strain `dim` axis (3D has a single constitutive branch).
 """
 
-from sinus_neumann import mms as sinus_neumann_mms
+from sinusoidal import mms as sinusoidal_mms
 
 from solid import (
     RESULTS_DIR,
@@ -17,13 +17,15 @@
 from output      import plot_convergence
 
 
-def convergence_study(elem_specs, mms, L, E, nu, nx_values):
+def convergence_study(elem_specs, mms, L, E, nu, nx_values,
+                      force_field, linear_solver):
     """
     Run a convergence series for each element type in elem_specs, write a
     per-(element) text table, and one shared plot with L²/H¹ for every
     element on the same axes.
 
     elem_specs : list of dicts with keys 'elem', 'label', 'l2_style', 'h1_style'
+    force_field : name of the FEM force field to test
     """
     print(f"\n  PoissonRatio = {nu}", flush=True)
 
@@ -36,7 +38,8 @@ def convergence_study(elem_specs, mms, L, E, nu, nx_values):
         hs, errors = run_convergence_series(
             nx_values  = nx_values,
             run_fn     = lambda nx, _e=elem: solve_solid(
-                _e, mms, L, E, nu, nx, nx, nx),
+                _e, mms, L, E, nu, nx, nx, nx,
+                force_field=force_field, linear_solver=linear_solver),
             h_fn       = lambda nx: L / (nx - 1),
             error_fns  = {
                 "L2": lambda sol, _e=elem: _e.compute_l2(sol, mms, L),
@@ -62,6 +65,8 @@ def convergence_study(elem_specs, mms, L, E, nu, nx_values):
     cfg  = load_params()
     L    = cfg["length"]
     E    = cfg["youngModulus"]
+    ff   = cfg["forceField"]
+    ls   = cfg["linearSolver"]
     conv = cfg["convergence"]
 
     specs = [
@@ -69,8 +74,9 @@ def convergence_study(elem_specs, mms, L, E, nu, nx_values):
          "l2_style": "bo-", "h1_style": "rs--"},
     ]
 
-    for mms in (sinus_neumann_mms,):
+    for mms in (sinusoidal_mms,):
         nx_vals = conv["nx_values"][mms.name]
         print(f"\n== {mms.name} ==")
         for nu in conv["nu_values"]:
-            convergence_study(specs, mms, L, E, nu, nx_vals)
+            convergence_study(specs, mms, L, E, nu, nx_vals,
+                              force_field=ff, linear_solver=ls)
diff --git a/examples/Freefem/MMS/3D/sinus_neumann.py b/examples/Freefem/MMS/3D/sinusoidal.py
similarity index 97%
rename from examples/Freefem/MMS/3D/sinus_neumann.py
rename to examples/Freefem/MMS/3D/sinusoidal.py
index 9cfa947e..f057552e 100644
--- a/examples/Freefem/MMS/3D/sinus_neumann.py
+++ b/examples/Freefem/MMS/3D/sinusoidal.py
@@ -26,7 +26,7 @@
 
 
 class SinusNeumann(MMSCase3D):
-    name       = "sinus_neumann"
+    name       = "sinusoidal"
     plot_label = (r"$u = A(\sin(\pi x/L)\sin(\pi y/L),\ "
                   r"\sin(\pi y/L)\sin(\pi z/L),\ "
                   r"\sin(\pi z/L)\sin(\pi x/L))$")
@@ -89,7 +89,7 @@ def find_corner(pred):
             for k, (xk, yk, zk) in enumerate(xyz):
                 if pred(xk, yk, zk):
                     return k
-            raise RuntimeError("sinus_neumann: BC corner not found")
+            raise RuntimeError("sinusoidal: BC corner not found")
 
         i_origin = find_corner(lambda x, y, z: x < eps     and y < eps     and z < eps)
         i_xL     = find_corner(lambda x, y, z: x > L - eps and y < eps     and z < eps)
diff --git a/examples/Freefem/MMS/3D/solid.py b/examples/Freefem/MMS/3D/solid.py
index 7fb68bec..03c01fcb 100644
--- a/examples/Freefem/MMS/3D/solid.py
+++ b/examples/Freefem/MMS/3D/solid.py
@@ -55,19 +55,24 @@ def compute(nodes, topology):
             nodes, conn, mms, L, E, nu, nx, ny, nz)
     return compute
 
-def build_solid_scene(rootNode, mms, element, L=1.0, E=1e6, nu=0.3,
+def build_solid_scene(rootNode, mms, element, force_field, linear_solver,
+                      L=1.0, E=1e6, nu=0.3,
                       nx=6, ny=6, nz=6, with_visual=True):
     """Build a SOFA scene for `mms` on the 3D `element` strategy.
 
     Returns (dofs, topology). Nodes and connectivity become available
     after `Sofa.Simulation.init(root)` runs, via
     `dofs.rest_position.array()` and `element.read_connectivity(topology)`.
+
+    force_field   : name of the FEM force field to test
+    linear_solver : dict {"type": <name>, "parameters": {...}} of the linear solver
     """
     rootNode.addObject("RequiredPlugin", pluginName=[
         "Elasticity",
         "Sofa.Component.Constraint.Projective",
         "Sofa.Component.Engine.Select",
         "Sofa.Component.LinearSolver.Direct",
+        "Sofa.Component.LinearSolver.Iterative",
         "Sofa.Component.MechanicalLoad",
         "Sofa.Component.ODESolver.Backward",
         "Sofa.Component.StateContainer",
@@ -97,8 +102,8 @@ def build_solid_scene(rootNode, mms, element, L=1.0, E=1e6, nu=0.3,
                     maxNbIterationsNewton=1,
                     absoluteResidualStoppingThreshold=1e-10,
                     printLog=False)
-    Solid.addObject("SparseLDLSolver", name="linearSolver",
-                    template="CompressedRowSparseMatrixd")
+    Solid.addObject(linear_solver["type"], name="linearSolver",
+                    **linear_solver["parameters"])
 
     dofs = Solid.addObject("MechanicalObject", name="dofs", template="Vec3d",
                            position=nodes_3d.tolist(),
@@ -106,7 +111,7 @@ def build_solid_scene(rootNode, mms, element, L=1.0, E=1e6, nu=0.3,
 
     topology = element.add_topology(Solid)
 
-    Solid.addObject("LinearSmallStrainFEMForceField", name="FEM", template="Vec3d",
+    Solid.addObject(force_field, name="FEM", template="Vec3d",
                     youngModulus=E, poissonRatio=nu, topology="@topology")
 
     mms.apply_bcs(Solid, nodes_3d, L)
@@ -130,12 +135,13 @@ def build_solid_scene(rootNode, mms, element, L=1.0, E=1e6, nu=0.3,
 # Simulation runner
 # ─────────────────────────────────────────────────────────────────────────────
 
-def solve_solid(elem, mms, L, E, nu, nx, ny, nz):
+def solve_solid(elem, mms, L, E, nu, nx, ny, nz, force_field, linear_solver):
     """Build, init, and run one static step. Returns a SolidSolution3D snapshot."""
     root = Sofa.Core.Node("root")
     dofs, topology = build_solid_scene(
         root, mms, elem, L=L, E=E, nu=nu,
-        nx=nx, ny=ny, nz=nz, with_visual=False
+        nx=nx, ny=ny, nz=nz, with_visual=False,
+        force_field=force_field, linear_solver=linear_solver
     )
     Sofa.Simulation.init(root)
     nodes_3d = dofs.rest_position.array().copy()
@@ -265,8 +271,11 @@ def run_reference_scene(elem, mms):
     L, E = cfg["length"], cfg["youngModulus"]
     nu   = ref["nu"]
     nx = ny = nz = ref["nx"]
+    ff   = cfg["forceField"]
+    ls   = cfg["linearSolver"]
 
-    sol = solve_solid(elem, mms, L, E, nu, nx, ny, nz)
+    sol = solve_solid(elem, mms, L, E, nu, nx, ny, nz,
+                      force_field=ff, linear_solver=ls)
     l2  = elem.compute_l2(sol, mms, L)
     h1  = elem.compute_h1(sol, mms, L)
 
@@ -291,7 +300,7 @@ def case_scene(mms, element):
     All parameters come from params.json (top-level + `reference` block). Each
     case file exposes:
         createScene = case_scene(mms, element_hex)
-    so that `runSofa sinus_neumann.py` loads the default scene.
+    so that `runSofa sinusoidal.py` loads the default scene.
     """
     def createScene(rootNode):
         cfg = load_params()
@@ -300,6 +309,7 @@ def createScene(rootNode):
                           L=cfg["length"], E=cfg["youngModulus"],
                           nu=ref["nu"],
                           nx=ref["nx"], ny=ref["nx"], nz=ref["nx"],
-                          with_visual=True)
+                          with_visual=True, force_field=cfg["forceField"],
+                          linear_solver=cfg["linearSolver"])
         return rootNode
     return createScene