diff --git a/examples/Freefem/validation/1D/comparaison_script.py b/examples/Freefem/validation/1D/comparaison_script.py new file mode 100644 index 00000000..0af9b342 --- /dev/null +++ b/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) + \ No newline at end of file diff --git a/examples/Freefem/validation/1D/freefem_bar_distributed.edp b/examples/Freefem/validation/1D/freefem_bar_distributed.edp new file mode 100644 index 00000000..a874ac49 --- /dev/null +++ b/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); diff --git a/examples/Freefem/validation/1D/params_distributed.json b/examples/Freefem/validation/1D/params_distributed.json new file mode 100644 index 00000000..e6cf7f68 --- /dev/null +++ b/examples/Freefem/validation/1D/params_distributed.json @@ -0,0 +1,11 @@ +{ + +"length":1, +"nx":10, +"q": 1000, +"youngModulus": 1000, +"poissonRatio":0.3 + + + +} \ No newline at end of file diff --git a/examples/Freefem/validation/1D/sofa_bar_distributed.py b/examples/Freefem/validation/1D/sofa_bar_distributed.py new file mode 100644 index 00000000..30a609d2 --- /dev/null +++ b/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"])) \ No newline at end of file