Corotational 2D. Rigid rotation only.

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,

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
     },

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(

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,16 +120,16 @@ 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(),
                           showObject=with_visual, showObjectScale=0.005 * L)
 
     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

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"
     },

examples/Freefem/MMS/2D/run_convergence.py

@@ -19,14 +19,17 @@
 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
     element on the same axes.
 
     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)