Corotational 2D. Rigid rotation only.

examples/Freefem/MMS/1D/bar.py

@@ -53,11 +53,14 @@ 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="LinearSmallStrainFEMForceField"):
     """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
     """
     root.addObject('RequiredPlugin', pluginName=[
         "Elasticity",
@@ -103,7 +106,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 +133,16 @@ 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"])
         return rootNode
     return createScene
 
 
-def solve_bar(mms, E_eff, nx):
+def solve_bar(mms, E_eff, nx, force_field="LinearSmallStrainFEMForceField"):
     """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)
     Sofa.Simulation.init(root)
     Sofa.Simulation.animate(root, root.dt.value)
     Bar   = root.Bar
@@ -180,7 +184,8 @@ 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"])
     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,7 @@
 {
     "length": 2.0,
     "youngModulus": 1e6,
+    "forceField": "LinearSmallStrainFEMForceField",
     "reference": {
         "nx": 10
     },

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

@@ -24,7 +24,8 @@
         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"]),
             h_fn       = lambda nx: 1.0 / (nx - 1),
             error_fns  = {
                 "L2": lambda sol, _m=mms: l2_error_1d(

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

@@ -67,11 +67,13 @@ def compute(nodes, topology):
     return compute
 
 def build_beam_scene(rootNode, mms, element, L=1.0, E=1e6, nu=0.3,
-                     nx=10, ny=10, with_visual=True, dim="2d"):
+                     nx=10, ny=10, with_visual=True, dim="2d",
+                     force_field="LinearSmallStrainFEMForceField"):
     """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
     Returns (dofs, topology). Nodes and connectivity become available
     after `Sofa.Simulation.init(root)` runs, via
     `dofs.rest_position.array()` and `element.read_connectivity(topology)`.
@@ -125,7 +127,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 +152,13 @@ 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, dim="2d",
+               force_field="LinearSmallStrainFEMForceField"):
     """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
     )
     Sofa.Simulation.init(root)
     # Read topology back from SOFA now that init has populated it.
@@ -246,9 +250,10 @@ 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"]
     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, dim=dim, force_field=ff)
     l2  = elem.compute_l2(sol, mms, L)
     h1  = elem.compute_h1(sol, mms, L)
 
@@ -281,6 +286,7 @@ 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"])
         return rootNode
     return createScene

examples/Freefem/MMS/2D/params.json

@@ -1,8 +1,9 @@
 {
     "length": 1.0,
     "youngModulus": 1e6,
+    "forceField": "CorotationalFEMForceField",
     "reference": {
-        "nx": 20,
+        "nx": 5,
         "nu": 0.0,
         "dim": "2d"
     },

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

@@ -19,14 +19,16 @@
 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, dim="2d",
+                      force_field="LinearSmallStrainFEMForceField"):
     """
     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
     """
     hyp = "plane strain" if dim == "3d" else "plane stress"
     print(f"\n  PoissonRatio = {nu}  ({hyp})", flush=True)
@@ -40,7 +42,7 @@ 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, dim=dim, force_field=force_field),
             h_fn       = lambda nx: L / (nx - 1),
             error_fns  = {
                 "L2": lambda sol, _e=elem: _e.compute_l2(sol, mms, L),
@@ -66,6 +68,7 @@ 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"]
     conv = cfg["convergence"]
 
     specs = [
@@ -81,4 +84,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, dim=DIM,
+                                  force_field=ff)

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"))