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_coro2d.py b/examples/Freefem/MMS/2D/test_coro2d.py
new file mode 100644
index 00000000..41176e49
--- /dev/null
+++ b/examples/Freefem/MMS/2D/test_coro2d.py
@@ -0,0 +1,319 @@
+"""
+u_d-only test : verifies that CorotationalFEMForceField (CR) and
+LinearSmallStrainFEMForceField (LS) produce identical results when the imposed
+displacement field is x = X + u_d(X), with no global rotation (R = I).
+
+Two test cases are considered:
+
+  1. SMOOTH u_d with a non-zero ANTISYMMETRIC gradient
+     (u_x = u_y = A sin(πx/L) sin(πy/L))
+     ===> CR and LS DO NOT converge to each other as h -> 0:
+        the relative difference reaches a finite plateau (expected physical behavior).
+
+  2. SMOOTH u_d with a purely SYMMETRIC, spatially varying gradient
+     (u_x = A sin(kx) cos(ky),  u_y = A cos(kx) sin(ky),  k = π/L)
+     ==> the local rotation is identically zero everywhere, as verified analytically:
+        ∂u_x/∂y = -Ak sin(kx) sin(ky) = ∂u_y/∂x  (equal cross-derivatives)
+     ==> any discrepancy between CR and LS should converge to zero as h -> 0.
+
+Known limitations
+-----------------
+- For nx >= 160, EulerExplicitSolver becomes unstable because the nodal mass
+  tends to zero (totalMass = 1 / N_nodes), causing the force values to diverge
+  after a single animation step. Therefore, the results are only considered
+  reliable for nx <= 80.
+  A future improvement will replace the EulerExplicitSolver + UniformMass
+  combination with a StaticSolver and a finite-difference Jacobian (Kfd)
+  for CorotationalFEMForceField.
+
+- Upstream warning: SOFAElementLinearSmallStrainFEMForceField registers
+  multiple default templates (Vec3d, Quad / Tetrahedron / Hexahedron) in
+  addition to the primary Vec3d, Edge default. This is an issue in the
+  Elasticity plugin itself and is outside the scope of this test.
+
+Run:
+    python test_coro2d.py
+"""
+
+import os
+import sys
+import numpy as np
+import matplotlib.pyplot as plt
+
+import Sofa
+import Sofa.Core
+import Sofa.Simulation
+import SofaRuntime  
+
+sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
+from beam import load_params, element_quad, RESULTS_DIR
+
+NX_LIST = [40, 80]   
+
+
+def build_scene(root, L, E, nu, nx, force_field):
+    root.addObject("RequiredPlugin", name="ScenePlugins", pluginName=[
+        "Elasticity",
+        "Sofa.Component.Mass",
+        "Sofa.Component.ODESolver.Forward",
+        "Sofa.Component.StateContainer",
+        "Sofa.Component.Topology.Container.Grid",
+        "Sofa.Component.Topology.Container.Dynamic",
+    ])
+    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:
+        beam.addObject("EulerExplicitSolver", name="odeSolver")
+        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")
+
+    return dofs
+
+
+# Displacement fields
+
+def u_smooth_rotational(X, amp=0.01, L=1.0, **_):
+    """u_x = u_y = A sin(pi*x/L) sin(pi*y/L)
+
+    Gradient antisymmetric part is non-zero:
+      du_x/dy =  Ak cos(pi*x/L) sin(pi*y/L)
+      du_y/dx =  Ak sin(pi*x/L) cos(pi*y/L)  
+
+    CR and LS are therefore expected NOT to converge toward each other as
+    h ----> 0: the relative discrepancy saturates at a finite value (physical).
+    """
+    k = np.pi / L
+    x, y = X[:, 0], X[:, 1]
+    field = amp * np.sin(k * x) * np.sin(k * y)
+    return np.stack([field, field.copy(), np.zeros_like(x)], axis=1)
+
+
+def u_smooth_irrotational(X, amp=0.01, L=1.0, **_):
+    """u_x =  A sin(kx) cos(ky),   u_y = A cos(kx) sin(ky),   k = pi/L
+
+    Analytically zero local rotation everywhere:
+      du_x/dy = -Ak sin(kx) sin(ky)
+      du_y/dx = -Ak sin(kx) sin(ky)   
+
+    ====> (du_x/dy - du_y/dx) / 2 = 0  exactly.
+
+    Any CR/LS discrepancy is purely from bilinear interpolation and
+    must vanish as h =====> 0.
+    """
+    k = np.pi / L
+    x, y = X[:, 0], X[:, 1]
+    ux = amp * np.sin(k * x) * np.cos(k * y)
+    uy = amp * np.cos(k * x) * np.sin(k * y)
+    return np.stack([ux, uy, np.zeros_like(x)], axis=1)
+ 
+def verify_irrotational(X, amp=0.01, L=1.0):
+    """Print max |du_x/dy - du_y/dx| at the nodal positions.
+
+    Must be 0.0 to machine precision: both partials equal
+    -Ak sin(kx) sin(ky), so their difference is identically zero.
+    """
+    k = np.pi / L
+    x, y = X[:, 0], X[:, 1]
+    dux_dy = -amp * k * np.sin(k * x) * np.sin(k * y)   
+    duy_dx = -amp * k * np.sin(k * x) * np.sin(k * y)   
+    err = np.max(np.abs(dux_dy - duy_dx))
+    print(f"  [verify]  max |∂u_x/∂y - ∂u_y/∂x| = {err:.3e}  "
+          f"(must be 0 — confirms zero continuous local rotation)")
+
+
+def run(force_field, u_field, nx_override=None):
+    cfg = load_params()
+    L, E   = cfg["length"], cfg["youngModulus"]
+    nu, nx = cfg["reference"]["nu"], cfg["reference"]["nx"]
+    if nx_override is not None:
+        nx = nx_override
+
+    root = Sofa.Core.Node("root")
+    dofs = build_scene(root, L, E, nu, nx, force_field)
+    Sofa.Simulation.init(root)
+
+    X  = dofs.rest_position.array().copy()
+    Xd = X + u_field(X, L=L)
+    with dofs.position.writeable() as p:
+        p[:] = Xd
+
+    Sofa.Simulation.animate(root, 1e-8)
+    f = dofs.force.array().copy()
+    Sofa.Simulation.unload(root)
+    return X, Xd, f
+
+ 
+def compare(u_field, nx_list, label):
+    diffs = []
+    for nx in nx_list:
+        _, _, f_cr = run("CorotationalFEMForceField",      u_field, nx)
+        _, _, f_ls = run("LinearSmallStrainFEMForceField", u_field, nx)
+        diff = np.max(np.linalg.norm(f_cr - f_ls, axis=1))
+        ref  = np.max(np.linalg.norm(f_ls,        axis=1))
+        rel  = diff / ref if ref > 0 else diff
+        diffs.append(rel)
+        print(f"  {label:14s} nx={nx:3d}   |f_CR - f_LS| / |f_LS| = {rel:.3e}")
+    return diffs
+
+
+
+
+def write_metrics(results, path):
+    with open(path, "w") as f:
+        f.write("CR vs LS : displacement-only test \n")
+        f.write("=" * 74 + "\n\n")
+        for label, (nx_list, diffs, expected) in results.items():
+            f.write(f"  [{label}]  {expected}\n")
+            for nx, d in zip(nx_list, diffs):
+                f.write(f"    nx={nx:3d}   |f_CR - f_LS| / |f_LS| = {d:.6e}\n")
+            f.write("\n")
+
+
+def plot_forces(results, path):
+    fig, axes = plt.subplots(1, 2, figsize=(14, 5))
+ 
+    ax0 = axes[0]
+    nx_list, diffs, _ = results["rotational"]
+    hs = [1.0 / nx for nx in nx_list]
+    ax0.loglog(hs, diffs,
+               marker="s", linestyle="--", color="tab:red",
+               label=r"rotational  $u_x = u_y$")
+
+    
+    plateau = diffs[-1]
+    ax0.axhline(y=plateau, color="gray", linestyle=":", linewidth=1.0)
+    ax0.text(hs[0] * 1.05, plateau * 1.3,
+             f"plateau ≈ {plateau:.2e}",
+             fontsize=8, color="gray")
+
+    ax0.invert_xaxis()
+    ax0.set_xlabel(r"$h = 1/n_x$  (mesh size, finer $\rightarrow$)")
+    ax0.set_ylabel(r"$\|f_{CR} - f_{LS}\| \;/\; \|f_{LS}\|$")
+    ax0.set_title("Case 1 — rotational field\n"
+                  r"$u_x = u_y = A\sin(\pi x/L)\sin(\pi y/L)$"
+                  "\n→ finite plateau (physical, not a bug)")
+    ax0.legend(loc="best")
+    ax0.grid(True, which="both", alpha=0.3)
+ 
+    ax1 = axes[1]
+    nx_list, diffs, _ = results["irrotational"]
+    hs = [1.0 / nx for nx in nx_list]
+    ax1.loglog(hs, diffs,
+               marker="^", linestyle="-", color="tab:blue",
+               label=r"irrotational  $\partial_y u_x = \partial_x u_y$")
+
+    floor = np.finfo(float).eps          
+    ax1.axhline(y=floor, color="gray", linestyle=":", linewidth=1.0)
+    ax1.text(max(hs) * 0.6, floor * 3.0,
+             r"$\varepsilon_{\rm machine} \approx 2.2\times10^{-16}$",
+             fontsize=8, color="gray")
+
+    ax1.invert_xaxis()
+    ax1.set_xlabel(r"$h = 1/n_x$  (mesh size, finer $\rightarrow$)")
+    ax1.set_ylabel(r"$\|f_{CR} - f_{LS}\| \;/\; \|f_{LS}\|$")
+    ax1.set_title("Case 2 — irrotational field\n"
+                  r"$u_x = A\sin(kx)\cos(ky),\; u_y = A\cos(kx)\sin(ky)$"
+                  "\n→ at machine precision for all $h$")
+    ax1.legend(loc="best")
+    ax1.grid(True, which="both", alpha=0.3)
+
+    fig.suptitle("CR vs LS — pure displacement test (no global rotation)",
+                 fontsize=12, fontweight="bold", y=1.01)
+    fig.tight_layout()
+    fig.savefig(path, dpi=150, bbox_inches="tight")
+    plt.close(fig)
+    print(f"  -> plot saved to {path}")
+
+ 
+
+if __name__ == "__main__":
+    results = {}
+
+    print("== Analytical verification of the irrotational field ==")
+    cfg = load_params()
+    xs = np.linspace(0, cfg["length"], 33)
+    ys = np.linspace(0, cfg["length"], 33)
+    Xg = np.array([[xi, yj, 0.0] for yj in ys for xi in xs])
+    verify_irrotational(Xg, L=cfg["length"])
+
+    # Case 1: rotational 
+    print("\n== 1) u_d smooth with non-zero antisymmetric part ==")
+    print("    (expected: NO convergence, relative diff saturates at a finite value)")
+    results["rotational"] = (
+        NX_LIST,
+        compare(u_smooth_rotational, NX_LIST, "rotational"),
+        "expected: saturation (physical, not a bug)",
+    )
+
+    # Case 2: irrotational 
+    print("\n== 2) u_d smooth, purely symmetric gradient  "
+          "(∂u_x/∂y = ∂u_y/∂x = -Ak sin(kx) sin(ky)) ==")
+    print("    (expected: convergence toward 0 as h -> 0)")
+    results["irrotational"] = (
+        NX_LIST,
+        compare(u_smooth_irrotational, NX_LIST, "irrotational"),
+        "expected: convergence to 0 as h -> 0",
+    )
+
+    os.makedirs(RESULTS_DIR, exist_ok=True)
+    write_metrics(results, os.path.join(RESULTS_DIR, "coro2d_metrics.txt"))
+    plot_forces(results,   os.path.join(RESULTS_DIR, "coro2d_convergence.png"))
+
+ 
+
+VISUALIZE_FIELD = "irrotational"          
+VISUALIZE_FF    = "CorotationalFEMForceField"   
+VISUALIZE_NX    = [4, 8, 16, 40, 80] 
+_FIELDS = {
+    "rotational":   u_smooth_rotational,
+    "irrotational": u_smooth_irrotational,
+}
+
+
+def createScene(root):
+    cfg = load_params()
+    L = cfg["length"]
+
+    u_field = _FIELDS[VISUALIZE_FIELD]
+
+    dofs = build_scene(root, L, cfg["youngModulus"],
+                       cfg["reference"]["nu"], VISUALIZE_NX, VISUALIZE_FF)
+
+    
+    root.addObject("RequiredPlugin", name="VisualPlugins", pluginName=[
+        "Sofa.GL.Component.Rendering3D",
+        "Sofa.GL.Component.Shader",
+        "Sofa.Component.Mapping.Linear",   
+    ])
+
+    class ApplyDisplacement(Sofa.Core.Controller):
+        def __init__(self, *args, **kwargs):
+            super().__init__(*args, **kwargs)
+            self._applied = False
+
+        def onAnimateBeginEvent(self, _):
+            if self._applied:
+                return
+            X = dofs.rest_position.array().copy()
+            Xd = X + u_field(X, L=L)
+            with dofs.position.writeable() as p:
+                p[:] = Xd
+            self._applied = True
+
+    root.addObject(ApplyDisplacement(name="ApplyDisplacement"))
+
+    beam = root.getChild("Beam")
+    visu = beam.addChild("Visual")
+    visu.addObject("OglModel", name="ogl", color="0.2 0.6 1.0 1.0")
+    visu.addObject("IdentityMapping", input="@../dofs", output="@ogl")
+
+    return root
\ No newline at end of file
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