1D validation: distributed load bar (SOFA vs FreeFEM vs analytical)

examples/Freefem/validation/1D/comparaison_script.py

@@ -0,0 +1,139 @@
+"""
+1D Bar Simulation - Distributed Load - Comparison File
+
+Physical case: bar fixed at x=0, free at x=L, uniform distributed load q.
+Analytical solution (E*u'' + q = 0, u(0)=0, u'(L)=0):
+
+    u(x) = (q/E) * (L*x - x**2/2)
+"""
+import json
+import os
+import sys
+import numpy as np
+import matplotlib.pyplot as plt
+from sofa_bar_distributed import sofaRun
+from pyfreefem import FreeFemRunner
+
+
+def _rms(a, b):
+    return np.linalg.norm(a - b) / np.sqrt(a.size)
+
+
+def _write_freefem_results(path, x, u):
+    x_final = x + u
+    with open(path, 'w') as f:
+        f.write(f"{'x_initial':>12}  {'x_final':>12}  {'u_x':>12}\n")
+        f.write("-" * 42 + "\n")
+        for xi, xf, ui in zip(x, x_final, u):
+            f.write(f"{xi:12.6f}  {xf:12.6f}  {ui:12.6f}\n")
+
+
+def _default_params_path():
+    return os.path.join(os.path.dirname(os.path.abspath(__file__)), "params_distributed.json")
+
+
+def plot_bar_geometry(x0, u_sofa, scale=None):
+
+    x_deformed = x0 + u_sofa
+    if scale is None:
+        span = x0[-1] - x0[0]
+        max_disp = np.max(np.abs(u_sofa)) or 1e-12
+        scale = 0.15 * span / max_disp  
+
+    fig, ax = plt.subplots(figsize=(9, 2.5))
+    y_rest, y_def = 1.0, 0.0
+
+    ax.plot(x0, [y_rest] * len(x0), 'o-', color="gray", label="Rest configuration")
+    ax.plot(x0 + scale * u_sofa, [y_def] * len(x0), 'o-', color="crimson",
+            label=f"Deformed (SOFA, ×{scale:.1f} for visibility)")
+    ax.plot([0, 0], [y_def - 0.3, y_rest + 0.3], 'k-', linewidth=3)  # fixed end marker
+
+    ax.set_yticks([y_def, y_rest])
+    ax.set_yticklabels(["deformed", "rest"])
+    ax.set_xlabel("x")
+    ax.set_title("1D Bar — Distributed Load — Rest vs Deformed (schematic)")
+    ax.legend(loc="upper right", fontsize=8)
+    ax.set_ylim(y_def - 0.6, y_rest + 0.6)
+    fig.tight_layout()
+    return fig
+
+
+if __name__ == "__main__":
+
+    config_file = sys.argv[1] if len(sys.argv) > 1 else _default_params_path()
+    with open(config_file) as f:
+        cfg = json.load(f)
+
+    length        = float(cfg["length"])
+    nx            = int(cfg["nx"])
+    q             = float(cfg["q"])
+    young_modulus = float(cfg["youngModulus"])
+    poisson_ratio = float(cfg["poissonRatio"])
+
+    # === Run FF ===
+    runner = FreeFemRunner("freefem_bar_distributed.edp")
+    exports = runner.execute({
+        'youngModulus': young_modulus,
+        'q': q,
+        'nx': nx,
+        'length': length,
+    })
+    x_ff = exports['xcoords']
+    u_ff = exports['u[]']
+
+    os.makedirs("results", exist_ok=True)
+    _write_freefem_results(os.path.join("results", "freefem_distributed_results.txt"), x_ff, u_ff)
+
+    # === Run SOFA ===
+    x_sofa, u_sofa = sofaRun(length=length
+            , q=q
+            , young_modulus=young_modulus
+            , poisson_ratio=poisson_ratio
+            , nx=nx)
+
+    # ========== analytical sol ===============
+    u_exact = (q / young_modulus) * (length * x_ff - x_ff**2 / 2.0)
+
+
+    rms_sofa_vs_ff    = _rms(u_sofa, u_ff)
+    rms_sofa_vs_exact = _rms(u_sofa, u_exact)
+
+
+    # --- Compare Results ---
+    with open("results/comparison_distributed_results.txt", 'w') as f:
+        header = f"{'x':>10}  {'u_exact':>12}  {'u_ff':>12}  {'u_sofa':>12}"
+        f.write(header + "\n")
+        f.write("-" * len(header) + "\n")
+        for x, ue, uff, us in zip(x_ff, u_exact, u_ff, u_sofa):
+            f.write(f"{x:10.4f}  {ue:12.6f}  {uff:12.6f}  {us:12.6f}\n")
+
+        f.write("\n")
+        f.write("RMS error norms (||a - b||_2 / sqrt(n))\n")
+        f.write("-" * 40 + "\n")
+        f.write(f"  RMS(sofa, ff)     = {rms_sofa_vs_ff:.6e}\n")
+        f.write(f"  RMS(sofa, exact)  = {rms_sofa_vs_exact:.6e}\n")
+
+
+    print("RMS error norms (||a - b||_2 / sqrt(n))")
+    print("-" * 40)
+    print(f"  RMS(sofa, ff)     = {rms_sofa_vs_ff:.6e}")
+    print(f"  RMS(sofa, exact)  = {rms_sofa_vs_exact:.6e}")
+
+
+    # ==== plot 
+    fig, ax = plt.subplots()
+    ax.plot(x_ff,   u_exact, label="Analytical", linestyle="--", color="black")
+    ax.plot(x_ff,   u_ff,    label="FreeFEM",    marker="o", markersize=4, linestyle="none")
+    ax.plot(x_sofa, u_sofa,  label="SOFA",       marker="x", markersize=5, linestyle="none")
+    ax.set_xlabel("x")
+    ax.set_ylabel("Displacement u(x)")
+    ax.set_title("1D Bar — Distributed Load — Displacement Comparison")
+    ax.legend()
+    fig.savefig("results/comparison_distributed_plot.png", dpi=150)
+    plt.close(fig)
+
+
+    fig_geo = plot_bar_geometry(x_sofa, u_sofa)
+    fig_geo.savefig("results/bar_geometry_distributed.png", dpi=150)
+    plt.close(fig_geo)
+

examples/Freefem/validation/1D/freefem_bar_distributed.edp

@@ -0,0 +1,31 @@
+// 1D Bar Simulation - Distributed Load
+IMPORT "io.edp"
+
+DEFAULT (length, 1.0)
+DEFAULT (youngModulus, 1000.0)
+DEFAULT (q, 1000.0)
+DEFAULT (nx, 10)
+
+real E = $youngModulus;
+real Q = $q;
+
+
+meshL Th = segment($nx - 1, [x*$length, 0, 0]);
+
+fespace Vh(Th, P1);
+Vh u, v;
+
+// FV 
+problem Distributed(u,v) =
+    int1d(Th)(E * dx(u) * dx(v))
+    - int1d(Th)(Q * v)
+    + on(1, u=0);
+
+Distributed;
+
+exportArray(u[]);
+
+
+real[int] xcoords(Th.nv);
+for (int i = 0; i < Th.nv; i++) xcoords[i] = Th(i).x;
+exportArray(xcoords);

examples/Freefem/validation/1D/params_distributed.json

@@ -0,0 +1,11 @@
+{
+
+"length":1,
+"nx":10,
+"q": 1000,
+"youngModulus": 1000,
+"poissonRatio":0.3
+
+
+
+}

examples/Freefem/validation/1D/sofa_bar_distributed.py

@@ -0,0 +1,179 @@
+"""
+1D Bar Simulation - Distributed Load - SOFA Scene File
+
+Physical case: bar fixed at x=0 (Dirichlet), free at x=L, subject to a
+uniform distributed load q per unit length (e.g. self-weight).
+
+Consistent nodal forces for a constant q on a uniform mesh of spacing h:
+    F_0     = q*h/2   (absorbed by the Dirichlet reaction, value irrelevant)
+    F_i     = q*h     for interior nodes
+    F_(N-1) = q*h/2   (free end)
+"""
+import json
+import os
+import sys
+import Sofa
+import Sofa.Core
+import Sofa.Simulation
+
+RESULTS_DIR = "results"
+
+
+class DisplacementExporter(Sofa.Core.Controller):
+
+    def __init__(self, dofs_node, output_file, *args, **kwargs):
+        super().__init__(*args, **kwargs)
+        self.dofs_node   = dofs_node
+        self.output_file = output_file
+        self.x_initial   = None
+        self.u_x         = None
+
+    def onSimulationInitDoneEvent(self, event):
+        self.x_initial = self.dofs_node.position.array().flatten().copy()
+
+    def onAnimateEndEvent(self, event):
+        x_final  = self.dofs_node.position.array().flatten()
+        self.u_x = x_final - self.x_initial
+
+        with open(self.output_file, 'w') as f:
+            f.write(f"{'x_initial':>12}  {'x_final':>12}  {'u_x':>12}\n")
+            f.write("-" * 42 + "\n")
+            for xi, xf, ui in zip(self.x_initial, x_final, self.u_x):
+                f.write(f"{xi:12.6f}  {xf:12.6f}  {ui:12.6f}\n")
+
+
+def _consistent_nodal_forces(q, h, nx):
+
+    forces = [q * h] * nx
+    forces[0]     = q * h / 2.0
+    forces[-1]    = q * h / 2.0
+    return [[f] for f in forces]
+
+
+def create_scene_args(rootNode, length, q, young_modulus, poisson_ratio, nx):
+    requiredPlugins = [
+        "Elasticity",
+        "Sofa.Component.Constraint.Projective",
+        "Sofa.Component.LinearSolver.Direct",
+        "Sofa.Component.MechanicalLoad",
+        "Sofa.Component.ODESolver.Backward",
+        "Sofa.Component.StateContainer",
+        "Sofa.Component.Topology.Container.Grid",
+        "Sofa.Component.Topology.Container.Dynamic",
+        "Sofa.Component.Visual",
+        "Sofa.GL.Component.Rendering3D",
+    ]
+
+    rootNode.addObject('RequiredPlugin', pluginName=requiredPlugins)
+    rootNode.addObject('DefaultAnimationLoop')
+    rootNode.addObject('VisualStyle', displayFlags=["showBehaviorModels", "showForceFields"])
+
+    h = length / (nx - 1)
+
+    Grid = rootNode.addChild('Grid')
+    Grid.addObject('RegularGridTopology'
+                    , name="grid"
+                    , nx=nx, ny=1, nz=1
+                    , min=[0., 0., 0.]
+                    , max=[length, 0., 0.])
+
+    with rootNode.addChild('Bar') as Bar:
+        Bar.addObject('NewtonRaphsonSolver'
+                    , name="newtonSolver"
+                    , printLog=True
+                    , warnWhenLineSearchFails=True
+                    , maxNbIterationsNewton=1
+                    , maxNbIterationsLineSearch=1
+                    , lineSearchCoefficient=1
+                    , relativeSuccessiveStoppingThreshold=0
+                    , absoluteResidualStoppingThreshold=1e-7
+                    , absoluteEstimateDifferenceThreshold=1e-12
+                    , relativeInitialStoppingThreshold=1e-12
+                    , relativeEstimateDifferenceThreshold=0
+                    )
+
+        Bar.addObject('SparseLDLSolver'
+                    , name="linearSolver"
+                    , template="CompressedRowSparseMatrixd")
+        Bar.addObject('StaticSolver'
+                    , name="staticSolver"
+                    , newtonSolver="@newtonSolver"
+                    , linearSolver="@linearSolver")
+
+        Bar.addObject('EdgeSetTopologyContainer'
+                    , name="topology"
+                    , edges="@../Grid/grid.edges"
+                    , position="@../Grid/grid.position")
+
+        dofs = Bar.addObject('MechanicalObject'
+                            , name="dofs"
+                            , template="Vec1d"
+                            , showObject=True
+                            , showObjectScale=0.02)
+
+        Bar.addObject('LinearSmallStrainFEMForceField'
+                    , name="FEM"
+                    , template="Vec1d"
+                    , youngModulus=young_modulus
+                    , poissonRatio=poisson_ratio
+                    , topology="@topology")
+
+        Bar.addObject('FixedProjectiveConstraint', indices="0")
+
+        Bar.addObject('ConstantForceField'
+                      , name="DistributedLoad"
+                      , indices=list(range(nx))
+                      , forces=_consistent_nodal_forces(q, h, nx))
+
+    os.makedirs(RESULTS_DIR, exist_ok=True)
+    exporter = rootNode.addObject(
+        DisplacementExporter(
+            dofs_node   = dofs,
+            output_file = os.path.join(RESULTS_DIR, "sofa_distributed_results.txt"),
+            name        = "exportCtrl"
+        )
+    )
+
+    return rootNode, exporter
+
+
+def _default_params_path():
+    """params_distributed.json next to this script, regardless of CWD."""
+    return os.path.join(os.path.dirname(os.path.abspath(__file__)), "params_distributed.json")
+
+
+def createScene(rootNode):
+    with open(_default_params_path()) as f:
+        cfg = json.load(f)
+    create_scene_args(rootNode
+                    , length=float(cfg["length"])
+                    , q=float(cfg["q"])
+                    , young_modulus=float(cfg["youngModulus"])
+                    , poisson_ratio=float(cfg["poissonRatio"])
+                    , nx=int(cfg["nx"]))
+    return rootNode
+
+
+def sofaRun(length, q, young_modulus, poisson_ratio, nx):
+    root = Sofa.Core.Node("root")
+    _, exporter = create_scene_args(root
+                                  , length=length
+                                  , q=q
+                                  , young_modulus=young_modulus
+                                  , poisson_ratio=poisson_ratio
+                                  , nx=nx)
+    Sofa.Simulation.init(root)
+    Sofa.Simulation.animate(root, root.dt.value)
+    return exporter.x_initial, exporter.u_x
+
+
+if __name__ == "__main__":
+    config_file = sys.argv[1] if len(sys.argv) > 1 else _default_params_path()
+    with open(config_file) as f:
+        cfg = json.load(f)
+
+    sofaRun(length=float(cfg["length"])
+            , q=float(cfg["q"])
+            , young_modulus=float(cfg["youngModulus"])
+            , poisson_ratio=float(cfg["poissonRatio"])
+            , nx=int(cfg["nx"]))