diff --git a/examples/Freefem/sofa/barre-traction.py b/examples/Freefem/sofa/barre-traction.py new file mode 100644 index 00000000..9b609899 --- /dev/null +++ b/examples/Freefem/sofa/barre-traction.py @@ -0,0 +1,151 @@ +""" +1D Bar Simulation - Traction Load +""" +import Sofa +import Sofa.Core +import Sofa.Simulation +import SofaRuntime +import numpy as np + +# Exporter les resultats +class DisplacementExporter(Sofa.Core.Controller): + """Récupère les déplacements après convergence et les écrit dans un .txt""" + + def __init__(self, dofs_node, output_file="sofa_displacement.txt", *args, **kwargs): + super().__init__(*args, **kwargs) + self.dofs_node = dofs_node + self.output_file = output_file + self.exported = False + + def onAnimateEndEvent(self, event): + """Appelé automatiquement par SOFA à la fin de chaque pas de temps""" + if self.exported: + return + + + pos = self.dofs_node.position.array() + + x_deform = pos[:, 0] + u_x = x_deform - self.x_initial + + + order = np.argsort(self.x_initial) + x0_sorted = self.x_initial[order] + ux_sorted = u_x[order] + + with open(self.output_file, 'w') as f: + f.write("x_initial u_x\n") + for xi, ui in zip(x0_sorted, ux_sorted): + f.write(f"{xi:.6f} {ui:.6f}\n") + + print(f"[DisplacementExporter] Déplacements exportés → {self.output_file}") + print(f" u_x(L) = {ux_sorted[-1]:.6f} (attendu : 1.0)") + self.exported = True + + def onSimulationInitDoneEvent(self, event): + """Stocke les positions initiales au démarrage""" + pos = self.dofs_node.position.array() + self.x_initial = pos[:, 0].copy() + print(f"[DisplacementExporter] {len(self.x_initial)} noeuds initialisés.") + + +# creation de la scene +def createScene(rootNode): + L = 1.0 + E_eff = 1000.0 + F = 1000.0 + N = 9 + + rootNode.addObject('DefaultAnimationLoop') + rootNode.addObject('VisualStyle', + displayFlags="showBehaviorModels showForceFields showWireframe") + + plugins = rootNode.addChild('plugins') + for plugin in [ + "Elasticity", + "Sofa.Component.Constraint.Projective", + "Sofa.Component.Engine.Select", + "Sofa.Component.LinearSolver.Direct", + "Sofa.Component.MechanicalLoad", + "Sofa.Component.ODESolver.Backward", + "Sofa.Component.StateContainer", + "Sofa.Component.Topology.Container.Dynamic", + "Sofa.Component.Topology.Container.Grid", + "Sofa.Component.Visual", + "Sofa.GL.Component.Rendering3D", + ]: + plugins.addObject('RequiredPlugin', name=plugin) + + bar = rootNode.addChild('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('RegularGridTopology', + name="grid", + min="0 0 0", + max=f"{L} 0 0", + n=f"{N} 1 1") + + dofs = bar.addObject('MechanicalObject', + template="Vec3d", + name="dofs", + showObject=True, + showObjectScale=0.02) + + edges = bar.addChild('edges') + edges.addObject('EdgeSetTopologyContainer', + name="topology", + position="@../dofs.position", + edges="@../grid.edges") + edges.addObject('LinearSmallStrainFEMForceField', + name="FEM", + youngModulus=E_eff, + poissonRatio=0.0, + topology="@topology") + + bar.addObject('BoxROI', + name="fixed_roi", + box="-0.01 -1 -1 0.01 1 1", + drawBoxes=False) + bar.addObject('FixedProjectiveConstraint', + indices="@fixed_roi.indices") + + bar.addObject('BoxROI', + name="load_roi", + box=f"{L-0.01} -1 -1 {L+0.01} 1 1", + drawBoxes=False) + bar.addObject('ConstantForceField', + indices="@load_roi.indices", + forces=f"{F} 0 0", + showArrowSize=1e-4) + + + rootNode.addObject( + DisplacementExporter( + dofs_node = dofs, + output_file = "sofa_displacement.txt", + name = "exportCtrl" + ) + ) + + return rootNode diff --git a/examples/Freefem/sofa/beam-2d.py b/examples/Freefem/sofa/beam-2d.py new file mode 100644 index 00000000..85b48287 --- /dev/null +++ b/examples/Freefem/sofa/beam-2d.py @@ -0,0 +1,202 @@ +""" +2D Beam Simulation - Plane Deformation Under Gravity +""" + +import Sofa +import Sofa.Core +import numpy as np +import csv +import os +import time + + +# ─── PARAMÈTRES +E = 21.5 +NU = 0.29 +GRAVITY = -0.05 +NX, NY = 141, 36 +TOTAL_MASS = 20.0 +OUTPUT_FILE = "sofa_displacement.txt" + + + +class DisplacementExporter(Sofa.Core.Controller): + """ + Recovers nodal displacements after convergence and writes them in a .txt file + """ + + def __init__(self, *args, **kwargs): + Sofa.Core.Controller.__init__(self, *args, **kwargs) + self.exported = False + self.node = kwargs.get("node", None) + self.iteration = 0 # compter les itérations + + def onAnimateEndEvent(self, event): + self.iteration += 1 + + if self.exported: + return + + if self.iteration < 50: + return + + root = self.getContext().getRootContext() + strain = root.getChild("beam_plane_strain") + if strain is None: + print("[DisplacementExporter] Nœud beam_plane_strain introuvable") + return + + dofs = strain.getObject("dofs") + if dofs is None: + print("[DisplacementExporter] MechanicalObject 'dofs' introuvable") + return + + # positions déformées + pos_def = np.array(dofs.position.value) + + # positions initiales + try: + pos_rest = np.array(dofs.rest_position.value) + except: + pos_rest = _build_initial_grid(NX, NY, y_offset=3.0) + + disp = pos_def - pos_rest + + # Afficher la convergence + max_uy = np.abs(disp[:, 1]).max() + print(f"[CONVERGENCE] NX={NX} NY={NY} nœuds={NX*NY} max|uy|={max_uy:.8f}") + + # export TXT + try: + with open(OUTPUT_FILE, "w") as f: + f.write("x y ux uy\n") + for i in range(len(pos_def)): + x = pos_rest[i, 0] + y = pos_rest[i, 1] - 3.0 # ramène à [0,2] + ux = disp[i, 0] + uy = disp[i, 1] + f.write(f"{x:.6f} {y:.6f} {ux:.8f} {uy:.8f}\n") + + print(f"[DisplacementExporter] ✓ Exporté {len(pos_def)} nœuds → {OUTPUT_FILE}") + self.exported = True + root.animate.value = False + + except Exception as e: + print(f"[DisplacementExporter] Erreur d'écriture: {e}") + + +def _build_initial_grid(nx, ny, y_offset=0.0): + """Reconstruit la grille régulière SOFA si rest_position n'est pas dispo.""" + xs = np.linspace(0, 10, nx) + ys = np.linspace(y_offset, y_offset + 2, ny) + pts = [] + for y in ys: + for x in xs: + pts.append([x, y, 0.0]) + return np.array(pts) + + +def createScene(rootNode): + """Point d'entrée SOFA — construit le graphe de scène.""" + + rootNode.name = "root" + rootNode.gravity = [0, GRAVITY, 0] + rootNode.dt = 1.0 + + # ── Plugins + rootNode.addObject("RequiredPlugin", name="Elasticity") + rootNode.addObject("RequiredPlugin", name="Sofa.Component.Constraint.Projective") + rootNode.addObject("RequiredPlugin", name="Sofa.Component.Engine.Select") + rootNode.addObject("RequiredPlugin", name="Sofa.Component.LinearSolver.Direct") + rootNode.addObject("RequiredPlugin", name="Sofa.Component.Mass") + rootNode.addObject("RequiredPlugin", name="Sofa.Component.ODESolver.Backward") + rootNode.addObject("RequiredPlugin", name="Sofa.Component.StateContainer") + rootNode.addObject("RequiredPlugin", name="Sofa.Component.Topology.Container.Dynamic") + rootNode.addObject("RequiredPlugin", name="Sofa.Component.Topology.Container.Grid") + rootNode.addObject("RequiredPlugin", name="Sofa.Component.Visual") + rootNode.addObject("RequiredPlugin", name="Sofa.GL.Component.Rendering3D") + + rootNode.addObject("DefaultAnimationLoop") + rootNode.addObject("VisualStyle", + displayFlags="showBehaviorModels showForceFields showWireframe") + + # ── Contrôleur d'export + rootNode.addObject(DisplacementExporter(name="exporter", node=rootNode)) + + # ── Nœud beam_plane_strain (Vec3d) + strain = rootNode.addChild("beam_plane_strain") + + # Réduire les logs pour plus de clarté + strain.addObject("NewtonRaphsonSolver", + name="newtonSolver", + printLog=False, + warnWhenLineSearchFails=True, + maxNbIterationsNewton=1, + maxNbIterationsLineSearch=1, + lineSearchCoefficient=1, + relativeSuccessiveStoppingThreshold=0, + absoluteResidualStoppingThreshold=1e-7, + absoluteEstimateDifferenceThreshold=1e-12, + relativeInitialStoppingThreshold=1e-12, + relativeEstimateDifferenceThreshold=0, + ) + strain.addObject("SparseLDLSolver", + name="linearSolver", + template="CompressedRowSparseMatrixd") + strain.addObject("StaticSolver", + name="staticSolver", + newtonSolver="@newtonSolver", + linearSolver="@linearSolver") + + # Grille régulière avec NX × NY points + strain.addObject("RegularGridTopology", + name="grid", + min=[0, 3, 0], + max=[10, 5, 0], + n=[NX, NY, 1]) + strain.addObject("MechanicalObject", + template="Vec3d", + name="dofs") + + # ── Sous-nœud triangles + triangles = strain.addChild("triangles") + triangles.addObject("TriangleSetTopologyContainer", + name="topology", src="@../grid") + triangles.addObject("TriangleSetTopologyModifier") + triangles.addObject("TriangleSetGeometryAlgorithms", + template="Vec3d", drawTriangles=True) + triangles.addObject("MeshMatrixMass", + totalMass=TOTAL_MASS, topology="@topology") + triangles.addObject("LinearSmallStrainFEMForceField", + name="FEM", + youngModulus=E, + poissonRatio=NU, + topology="@topology") + + # ── Condition aux limites : encastrement à gauche + strain.addObject("BoxROI", + name="fixed_roi", + template="Vec3d", + box=[-0.01, 2.99, -1.0, 0.01, 5.01, 1.0], + drawBoxes=True) + strain.addObject("FixedProjectiveConstraint", + template="Vec3d", + indices="@fixed_roi.indices") + + print(f"[createScene] ✓ Scène créée avec NX={NX}, NY={NY}, nœuds={NX*NY}") + + return rootNode + + +start = time.time() + +# simulation + +end = time.time() +print("Temps de calcul :", end - start) + + +if __name__ == "__main__": + # Pour exécution directe (non-SOFA) + print(f"Paramètres: NX={NX}, NY={NY}, output={OUTPUT_FILE}") + print(f"Total nœuds: {NX*NY}")