Add FreeFEM/SOFA validation project

examples/Freefem/projet_validation/cases/bar_1d/params.json

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

examples/Freefem/projet_validation/cases/bar_1d/run.py

@@ -0,0 +1,71 @@
+"""
+Cas 1D — Barre encastrée-libre, charge répartie.
+
+Toute la logique (paramètres, FreeFem, SOFA, solution analytique) est
+encapsulée ici. Le notebook n'appelle que run_case().
+"""
+import json
+import os
+
+from common import env_setup  # noqa: F401 -- corrige PATH pour stdbuf (Windows)
+from pyfreefem import FreeFemRunner
+
+from .sofa_scene import sofaRun
+
+HERE = os.path.dirname(os.path.abspath(__file__))
+
+LABEL = "Barre 1D — charge répartie"
+
+MATH_DESCRIPTION = """
+Équation forte : E·u''(x) + q = 0, avec u(0)=0 (encastrement) et u'(L)=0 (extrémité libre).
+Solution analytique : u(x) = (q/E) · (L·x - x²/2)
+"""
+
+
+def _load_params():
+    with open(os.path.join(HERE, "params.json")) as f:
+        return json.load(f)
+
+
+def _u_exact(x, q, E, L):
+    return (q / E) * (L * x - x**2 / 2.0)
+
+
+def run_case():
+    """
+    Exécute FreeFem et SOFA sur le cas 1D distribué, calcule la solution
+    analytique, et renvoie un dict de résultats prêt pour common.metrics
+    et common.plotting.
+    """
+    cfg = _load_params()
+    length, nx = float(cfg["length"]), int(cfg["nx"])
+    q = float(cfg["q"])
+    young_modulus = float(cfg["youngModulus"])
+    poisson_ratio = float(cfg["poissonRatio"])
+
+    # --- FreeFem++ ---
+    runner = FreeFemRunner(os.path.join(HERE, "freefem_bar_distributed.edp"))
+    exports = runner.execute({
+        "youngModulus": young_modulus,
+        "q": q,
+        "nx": nx,
+        "length": length,
+    }, verbosity=0)
+    x_ff, u_ff = exports["xcoords"], exports["u[]"]
+
+    # --- SOFA ---
+    x_sofa, u_sofa = sofaRun(
+        length=length, q=q,
+        young_modulus=young_modulus, poisson_ratio=poisson_ratio, nx=nx,
+    )
+
+    # --- Analytique ---
+    u_ana = _u_exact(x_ff, q, young_modulus, length)
+
+    return {
+        "label": LABEL,
+        "params": cfg,
+        "x_ff": x_ff, "u_ff": u_ff,
+        "x_sofa": x_sofa, "u_sofa": u_sofa,
+        "u_ana": u_ana,
+    }

examples/Freefem/projet_validation/cases/bar_1d/sofa_scene.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.json next to this script, regardless of CWD."""
+    return os.path.join(os.path.dirname(os.path.abspath(__file__)), "params.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"]))