From 74c0368d080c9077d5d36ae40bb00937846f9236 Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Mon, 22 Jun 2026 11:11:16 +0200 Subject: [PATCH 1/3] [COR] Allow choosing name for ForceField --- examples/Freefem/MMS/1D/bar.py | 17 +++++++++++------ examples/Freefem/MMS/1D/params.json | 1 + examples/Freefem/MMS/1D/run_convergence.py | 3 ++- examples/Freefem/MMS/2D/beam.py | 18 ++++++++++++------ examples/Freefem/MMS/2D/params.json | 1 + examples/Freefem/MMS/2D/run_convergence.py | 10 +++++++--- examples/Freefem/MMS/3D/params.json | 1 + examples/Freefem/MMS/3D/run_convergence.py | 9 ++++++--- examples/Freefem/MMS/3D/solid.py | 17 +++++++++++------ 9 files changed, 52 insertions(+), 25 deletions(-) diff --git a/examples/Freefem/MMS/1D/bar.py b/examples/Freefem/MMS/1D/bar.py index 3904cb11..1e813acf 100644 --- a/examples/Freefem/MMS/1D/bar.py +++ b/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, diff --git a/examples/Freefem/MMS/1D/params.json b/examples/Freefem/MMS/1D/params.json index c094267f..5dad8e81 100644 --- a/examples/Freefem/MMS/1D/params.json +++ b/examples/Freefem/MMS/1D/params.json @@ -1,6 +1,7 @@ { "length": 2.0, "youngModulus": 1e6, + "forceField": "LinearSmallStrainFEMForceField", "reference": { "nx": 10 }, diff --git a/examples/Freefem/MMS/1D/run_convergence.py b/examples/Freefem/MMS/1D/run_convergence.py index d063d38d..abd59f35 100644 --- a/examples/Freefem/MMS/1D/run_convergence.py +++ b/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( diff --git a/examples/Freefem/MMS/2D/beam.py b/examples/Freefem/MMS/2D/beam.py index 420ccee4..185b71ec 100644 --- a/examples/Freefem/MMS/2D/beam.py +++ b/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 diff --git a/examples/Freefem/MMS/2D/params.json b/examples/Freefem/MMS/2D/params.json index f98c79a1..fb86704f 100644 --- a/examples/Freefem/MMS/2D/params.json +++ b/examples/Freefem/MMS/2D/params.json @@ -1,6 +1,7 @@ { "length": 1.0, "youngModulus": 1e6, + "forceField": "LinearSmallStrainFEMForceField", "reference": { "nx": 20, "nu": 0.0, diff --git a/examples/Freefem/MMS/2D/run_convergence.py b/examples/Freefem/MMS/2D/run_convergence.py index bad05b1f..ad2ad046 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, 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 @@ -27,6 +28,7 @@ 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 """ 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) diff --git a/examples/Freefem/MMS/3D/params.json b/examples/Freefem/MMS/3D/params.json index 6a902e26..cb4522df 100644 --- a/examples/Freefem/MMS/3D/params.json +++ b/examples/Freefem/MMS/3D/params.json @@ -1,6 +1,7 @@ { "length": 1.0, "youngModulus": 1e6, + "forceField": "LinearSmallStrainFEMForceField", "reference": { "nx": 6, "nu": 0.3 diff --git a/examples/Freefem/MMS/3D/run_convergence.py b/examples/Freefem/MMS/3D/run_convergence.py index 44591a83..777ab319 100644 --- a/examples/Freefem/MMS/3D/run_convergence.py +++ b/examples/Freefem/MMS/3D/run_convergence.py @@ -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="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' + force_field : name of the FEM force field to test """ print(f"\n PoissonRatio = {nu}", flush=True) @@ -36,7 +38,7 @@ 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), h_fn = lambda nx: L / (nx - 1), error_fns = { "L2": lambda sol, _e=elem: _e.compute_l2(sol, mms, L), @@ -62,6 +64,7 @@ def convergence_study(elem_specs, mms, L, E, nu, nx_values): cfg = load_params() L = cfg["length"] E = cfg["youngModulus"] + ff = cfg["forceField"] conv = cfg["convergence"] specs = [ @@ -73,4 +76,4 @@ def convergence_study(elem_specs, mms, L, E, nu, nx_values): 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) diff --git a/examples/Freefem/MMS/3D/solid.py b/examples/Freefem/MMS/3D/solid.py index 7fb68bec..b24f9583 100644 --- a/examples/Freefem/MMS/3D/solid.py +++ b/examples/Freefem/MMS/3D/solid.py @@ -56,12 +56,15 @@ def compute(nodes, topology): return compute def build_solid_scene(rootNode, mms, element, L=1.0, E=1e6, nu=0.3, - nx=6, ny=6, nz=6, with_visual=True): + nx=6, ny=6, nz=6, with_visual=True, + force_field="LinearSmallStrainFEMForceField"): """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 """ rootNode.addObject("RequiredPlugin", pluginName=[ "Elasticity", @@ -106,7 +109,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 +133,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="LinearSmallStrainFEMForceField"): """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 ) Sofa.Simulation.init(root) nodes_3d = dofs.rest_position.array().copy() @@ -265,8 +269,9 @@ def run_reference_scene(elem, mms): L, E = cfg["length"], cfg["youngModulus"] nu = ref["nu"] nx = ny = nz = ref["nx"] + ff = cfg["forceField"] - 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) l2 = elem.compute_l2(sol, mms, L) h1 = elem.compute_h1(sol, mms, L) @@ -300,6 +305,6 @@ 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"]) return rootNode return createScene From 4c1fc3d6df72937e4eedea31aab3ebaf3bc6e85c Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Mon, 22 Jun 2026 11:45:49 +0200 Subject: [PATCH 2/3] [2D] Introduce amplitude in 2D trigonometric function --- examples/Freefem/MMS/2D/trigonometric.py | 25 +++++++++++++++--------- 1 file changed, 16 insertions(+), 9 deletions(-) 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() From ca7b3f2406ab2db46dacb2b13b2234bb97c4ae82 Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Mon, 22 Jun 2026 14:35:48 +0200 Subject: [PATCH 3/3] [2D] Step 1: Rigid-body rotation test --- examples/Freefem/MMS/2D/params.json | 4 +- .../Freefem/MMS/2D/test_rigid_rotation.py | 177 ++++++++++++++++++ 2 files changed, 179 insertions(+), 2 deletions(-) create mode 100644 examples/Freefem/MMS/2D/test_rigid_rotation.py diff --git a/examples/Freefem/MMS/2D/params.json b/examples/Freefem/MMS/2D/params.json index fb86704f..5d0bc053 100644 --- a/examples/Freefem/MMS/2D/params.json +++ b/examples/Freefem/MMS/2D/params.json @@ -1,9 +1,9 @@ { "length": 1.0, "youngModulus": 1e6, - "forceField": "LinearSmallStrainFEMForceField", + "forceField": "CorotationalFEMForceField", "reference": { - "nx": 20, + "nx": 5, "nu": 0.0, "dim": "2d" }, 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"))