From 681d89de33374110c50a4c8c63af6fa2f7094855 Mon Sep 17 00:00:00 2001 From: fatima Date: Thu, 16 Jul 2026 17:09:47 +0200 Subject: [PATCH] Add 2D beam plane stress, distributed load, SOFA vs FreeFEM --- .../2D/2D_planeStress/beam2d_tri.msh | 330 ++++++++++++++++++ .../comparaison_script_2d_cp.py | 117 +++++++ .../freefem_beam2d_distributed_cp.edp | 49 +++ .../params_beam2d_distributed.json | 7 + .../sofa_beam2d_distributed_cp.py | 222 ++++++++++++ 5 files changed, 725 insertions(+) create mode 100644 examples/Freefem/validation/2D/2D_planeStress/beam2d_tri.msh create mode 100644 examples/Freefem/validation/2D/2D_planeStress/comparaison_script_2d_cp.py create mode 100644 examples/Freefem/validation/2D/2D_planeStress/freefem_beam2d_distributed_cp.edp create mode 100644 examples/Freefem/validation/2D/2D_planeStress/params_beam2d_distributed.json create mode 100644 examples/Freefem/validation/2D/2D_planeStress/sofa_beam2d_distributed_cp.py diff --git a/examples/Freefem/validation/2D/2D_planeStress/beam2d_tri.msh b/examples/Freefem/validation/2D/2D_planeStress/beam2d_tri.msh new file mode 100644 index 00000000..69c0bf53 --- /dev/null +++ b/examples/Freefem/validation/2D/2D_planeStress/beam2d_tri.msh @@ -0,0 +1,330 @@ +$MeshFormat +2.2 0 8 +$EndMeshFormat +$PhysicalNames +5 +1 1 "Fixed" +1 2 "Bottom" +1 3 "Right" +1 4 "Top" +2 5 "Beam" +$EndPhysicalNames +$Nodes +105 +1 0 0 0 +2 1 0 0 +3 1 0.2 0 +4 0 0.2 0 +5 0 0.1500000000001384 0 +6 0 0.100000000000274 0 +7 0 0.05000000000013835 0 +8 0.04999999999989938 0 0 +9 0.09999999999981413 0 0 +10 0.1499999999997028 0 0 +11 0.1999999999995569 0 0 +12 0.2499999999994109 0 0 +13 0.299999999999265 0 0 +14 0.3499999999991191 0 0 +15 0.3999999999989731 0 0 +16 0.4499999999988272 0 0 +17 0.4999999999986921 0 0 +18 0.5499999999988131 0 0 +19 0.599999999998945 0 0 +20 0.6499999999990769 0 0 +21 0.6999999999992088 0 0 +22 0.7499999999993406 0 0 +23 0.7999999999994725 0 0 +24 0.8499999999996044 0 0 +25 0.8999999999997362 0 0 +26 0.9499999999998682 0 0 +27 1 0.04999999999988333 0 +28 1 0.0999999999997371 0 +29 1 0.1499999999998672 0 +30 0.9499999999997918 0.2 0 +31 0.8999999999995836 0.2 0 +32 0.849999999999653 0.2 0 +33 0.8 0.2 0 +34 0.7500000000003471 0.2 0 +35 0.7000000000006938 0.2 0 +36 0.6500000000010406 0.2 0 +37 0.6000000000013874 0.2 0 +38 0.5500000000017343 0.2 0 +39 0.5000000000020595 0.2 0 +40 0.450000000001873 0.2 0 +41 0.400000000001665 0.2 0 +42 0.3500000000014568 0.2 0 +43 0.3000000000012487 0.2 0 +44 0.2500000000010406 0.2 0 +45 0.2000000000008325 0.2 0 +46 0.1500000000006244 0.2 0 +47 0.1000000000004162 0.2 0 +48 0.0500000000002081 0.2 0 +49 0.05000000000013093 0.1500000000001248 0 +50 0.1000000000002657 0.1500000000001113 0 +51 0.150000000000394 0.1500000000000977 0 +52 0.2000000000005136 0.1500000000000841 0 +53 0.2500000000006332 0.1500000000000705 0 +54 0.3000000000007528 0.1500000000000571 0 +55 0.3500000000008724 0.1500000000000435 0 +56 0.4000000000009921 0.1500000000000299 0 +57 0.4500000000011116 0.1500000000000164 0 +58 0.5000000000012177 0.1500000000000028 0 +59 0.5500000000010041 0.1499999999999892 0 +60 0.6000000000007768 0.1499999999999757 0 +61 0.6500000000005497 0.1499999999999621 0 +62 0.7000000000003226 0.1499999999999486 0 +63 0.7500000000000955 0.149999999999935 0 +64 0.7999999999998681 0.1499999999999214 0 +65 0.849999999999641 0.1499999999999079 0 +66 0.8999999999996218 0.1499999999998944 0 +67 0.949999999999811 0.1499999999998808 0 +68 0.05000000000005374 0.1000000000002471 0 +69 0.1000000000001152 0.1000000000002203 0 +70 0.1500000000001636 0.1000000000001934 0 +71 0.2000000000001947 0.1000000000001666 0 +72 0.2500000000002258 0.1000000000001398 0 +73 0.3000000000002568 0.1000000000001129 0 +74 0.350000000000288 0.1000000000000861 0 +75 0.4000000000003191 0.1000000000000592 0 +76 0.4500000000003501 0.1000000000000324 0 +77 0.5000000000003757 0.1000000000000056 0 +78 0.5500000000002736 0.09999999999997869 0 +79 0.6000000000001662 0.09999999999995186 0 +80 0.6500000000000588 0.09999999999992504 0 +81 0.6999999999999513 0.09999999999989816 0 +82 0.7499999999998439 0.09999999999987133 0 +83 0.7999999999997363 0.09999999999984448 0 +84 0.8499999999996287 0.09999999999981764 0 +85 0.89999999999966 0.09999999999979081 0 +86 0.9499999999998299 0.09999999999976396 0 +87 0.04999999999997656 0.0500000000001256 0 +88 0.09999999999996464 0.05000000000011284 0 +89 0.1499999999999332 0.0500000000001001 0 +90 0.1999999999998758 0.05000000000008736 0 +91 0.2499999999998184 0.0500000000000746 0 +92 0.299999999999761 0.05000000000006185 0 +93 0.3499999999997035 0.0500000000000491 0 +94 0.3999999999996461 0.05000000000003635 0 +95 0.4499999999995887 0.0500000000000236 0 +96 0.4999999999995339 0.05000000000001084 0 +97 0.5499999999995435 0.04999999999999809 0 +98 0.5999999999995556 0.04999999999998535 0 +99 0.6499999999995678 0.04999999999997258 0 +100 0.6999999999995801 0.04999999999995983 0 +101 0.7499999999995923 0.04999999999994709 0 +102 0.7999999999996044 0.04999999999993435 0 +103 0.8499999999996166 0.04999999999992158 0 +104 0.8999999999996982 0.04999999999990884 0 +105 0.9499999999998492 0.04999999999989609 0 +$EndNodes +$Elements +208 +1 1 2 1 1 4 5 +2 1 2 1 1 5 6 +3 1 2 1 1 6 7 +4 1 2 1 1 7 1 +5 1 2 2 2 1 8 +6 1 2 2 2 8 9 +7 1 2 2 2 9 10 +8 1 2 2 2 10 11 +9 1 2 2 2 11 12 +10 1 2 2 2 12 13 +11 1 2 2 2 13 14 +12 1 2 2 2 14 15 +13 1 2 2 2 15 16 +14 1 2 2 2 16 17 +15 1 2 2 2 17 18 +16 1 2 2 2 18 19 +17 1 2 2 2 19 20 +18 1 2 2 2 20 21 +19 1 2 2 2 21 22 +20 1 2 2 2 22 23 +21 1 2 2 2 23 24 +22 1 2 2 2 24 25 +23 1 2 2 2 25 26 +24 1 2 2 2 26 2 +25 1 2 3 3 2 27 +26 1 2 3 3 27 28 +27 1 2 3 3 28 29 +28 1 2 3 3 29 3 +29 1 2 4 4 3 30 +30 1 2 4 4 30 31 +31 1 2 4 4 31 32 +32 1 2 4 4 32 33 +33 1 2 4 4 33 34 +34 1 2 4 4 34 35 +35 1 2 4 4 35 36 +36 1 2 4 4 36 37 +37 1 2 4 4 37 38 +38 1 2 4 4 38 39 +39 1 2 4 4 39 40 +40 1 2 4 4 40 41 +41 1 2 4 4 41 42 +42 1 2 4 4 42 43 +43 1 2 4 4 43 44 +44 1 2 4 4 44 45 +45 1 2 4 4 45 46 +46 1 2 4 4 46 47 +47 1 2 4 4 47 48 +48 1 2 4 4 48 4 +49 2 2 5 1 4 5 48 +50 2 2 5 1 48 5 49 +51 2 2 5 1 48 49 47 +52 2 2 5 1 47 49 50 +53 2 2 5 1 47 50 46 +54 2 2 5 1 46 50 51 +55 2 2 5 1 46 51 45 +56 2 2 5 1 45 51 52 +57 2 2 5 1 45 52 44 +58 2 2 5 1 44 52 53 +59 2 2 5 1 44 53 43 +60 2 2 5 1 43 53 54 +61 2 2 5 1 43 54 42 +62 2 2 5 1 42 54 55 +63 2 2 5 1 42 55 41 +64 2 2 5 1 41 55 56 +65 2 2 5 1 41 56 40 +66 2 2 5 1 40 56 57 +67 2 2 5 1 40 57 39 +68 2 2 5 1 39 57 58 +69 2 2 5 1 39 58 38 +70 2 2 5 1 38 58 59 +71 2 2 5 1 38 59 37 +72 2 2 5 1 37 59 60 +73 2 2 5 1 37 60 36 +74 2 2 5 1 36 60 61 +75 2 2 5 1 36 61 35 +76 2 2 5 1 35 61 62 +77 2 2 5 1 35 62 34 +78 2 2 5 1 34 62 63 +79 2 2 5 1 34 63 33 +80 2 2 5 1 33 63 64 +81 2 2 5 1 33 64 32 +82 2 2 5 1 32 64 65 +83 2 2 5 1 32 65 31 +84 2 2 5 1 31 65 66 +85 2 2 5 1 31 66 30 +86 2 2 5 1 30 66 67 +87 2 2 5 1 30 67 3 +88 2 2 5 1 3 67 29 +89 2 2 5 1 5 6 49 +90 2 2 5 1 49 6 68 +91 2 2 5 1 49 68 50 +92 2 2 5 1 50 68 69 +93 2 2 5 1 50 69 51 +94 2 2 5 1 51 69 70 +95 2 2 5 1 51 70 52 +96 2 2 5 1 52 70 71 +97 2 2 5 1 52 71 53 +98 2 2 5 1 53 71 72 +99 2 2 5 1 53 72 54 +100 2 2 5 1 54 72 73 +101 2 2 5 1 54 73 55 +102 2 2 5 1 55 73 74 +103 2 2 5 1 55 74 56 +104 2 2 5 1 56 74 75 +105 2 2 5 1 56 75 57 +106 2 2 5 1 57 75 76 +107 2 2 5 1 57 76 58 +108 2 2 5 1 58 76 77 +109 2 2 5 1 58 77 59 +110 2 2 5 1 59 77 78 +111 2 2 5 1 59 78 60 +112 2 2 5 1 60 78 79 +113 2 2 5 1 60 79 61 +114 2 2 5 1 61 79 80 +115 2 2 5 1 61 80 62 +116 2 2 5 1 62 80 81 +117 2 2 5 1 62 81 63 +118 2 2 5 1 63 81 82 +119 2 2 5 1 63 82 64 +120 2 2 5 1 64 82 83 +121 2 2 5 1 64 83 65 +122 2 2 5 1 65 83 84 +123 2 2 5 1 65 84 66 +124 2 2 5 1 66 84 85 +125 2 2 5 1 66 85 67 +126 2 2 5 1 67 85 86 +127 2 2 5 1 67 86 29 +128 2 2 5 1 29 86 28 +129 2 2 5 1 6 7 68 +130 2 2 5 1 68 7 87 +131 2 2 5 1 68 87 69 +132 2 2 5 1 69 87 88 +133 2 2 5 1 69 88 70 +134 2 2 5 1 70 88 89 +135 2 2 5 1 70 89 71 +136 2 2 5 1 71 89 90 +137 2 2 5 1 71 90 72 +138 2 2 5 1 72 90 91 +139 2 2 5 1 72 91 73 +140 2 2 5 1 73 91 92 +141 2 2 5 1 73 92 74 +142 2 2 5 1 74 92 93 +143 2 2 5 1 74 93 75 +144 2 2 5 1 75 93 94 +145 2 2 5 1 75 94 76 +146 2 2 5 1 76 94 95 +147 2 2 5 1 76 95 77 +148 2 2 5 1 77 95 96 +149 2 2 5 1 77 96 78 +150 2 2 5 1 78 96 97 +151 2 2 5 1 78 97 79 +152 2 2 5 1 79 97 98 +153 2 2 5 1 79 98 80 +154 2 2 5 1 80 98 99 +155 2 2 5 1 80 99 81 +156 2 2 5 1 81 99 100 +157 2 2 5 1 81 100 82 +158 2 2 5 1 82 100 101 +159 2 2 5 1 82 101 83 +160 2 2 5 1 83 101 102 +161 2 2 5 1 83 102 84 +162 2 2 5 1 84 102 103 +163 2 2 5 1 84 103 85 +164 2 2 5 1 85 103 104 +165 2 2 5 1 85 104 86 +166 2 2 5 1 86 104 105 +167 2 2 5 1 86 105 28 +168 2 2 5 1 28 105 27 +169 2 2 5 1 7 1 87 +170 2 2 5 1 87 1 8 +171 2 2 5 1 87 8 88 +172 2 2 5 1 88 8 9 +173 2 2 5 1 88 9 89 +174 2 2 5 1 89 9 10 +175 2 2 5 1 89 10 90 +176 2 2 5 1 90 10 11 +177 2 2 5 1 90 11 91 +178 2 2 5 1 91 11 12 +179 2 2 5 1 91 12 92 +180 2 2 5 1 92 12 13 +181 2 2 5 1 92 13 93 +182 2 2 5 1 93 13 14 +183 2 2 5 1 93 14 94 +184 2 2 5 1 94 14 15 +185 2 2 5 1 94 15 95 +186 2 2 5 1 95 15 16 +187 2 2 5 1 95 16 96 +188 2 2 5 1 96 16 17 +189 2 2 5 1 96 17 97 +190 2 2 5 1 97 17 18 +191 2 2 5 1 97 18 98 +192 2 2 5 1 98 18 19 +193 2 2 5 1 98 19 99 +194 2 2 5 1 99 19 20 +195 2 2 5 1 99 20 100 +196 2 2 5 1 100 20 21 +197 2 2 5 1 100 21 101 +198 2 2 5 1 101 21 22 +199 2 2 5 1 101 22 102 +200 2 2 5 1 102 22 23 +201 2 2 5 1 102 23 103 +202 2 2 5 1 103 23 24 +203 2 2 5 1 103 24 104 +204 2 2 5 1 104 24 25 +205 2 2 5 1 104 25 105 +206 2 2 5 1 105 25 26 +207 2 2 5 1 105 26 27 +208 2 2 5 1 27 26 2 +$EndElements diff --git a/examples/Freefem/validation/2D/2D_planeStress/comparaison_script_2d_cp.py b/examples/Freefem/validation/2D/2D_planeStress/comparaison_script_2d_cp.py new file mode 100644 index 00000000..53bb05e1 --- /dev/null +++ b/examples/Freefem/validation/2D/2D_planeStress/comparaison_script_2d_cp.py @@ -0,0 +1,117 @@ +import json +import os +import sys +import numpy as np +import matplotlib.pyplot as plt +import matplotlib.tri as mtri + +from sofa_beam2d_distributed_cp import sofaRun, read_gmsh_2d +from pyfreefem import FreeFemRunner + + +def _rms(a, b): + """Discrete RMS, normalized by sqrt(n).""" + return np.linalg.norm(a - b) / np.sqrt(a.size) + + +def _default_params_path(): + return os.path.join(os.path.dirname(os.path.abspath(__file__)), "params_beam2d_distributed.json") + + +def _default_mesh_path(): + return os.path.join(os.path.dirname(os.path.abspath(__file__)), "beam2d_tri.msh") + + +def _pair_by_coordinates(x_a, y_a, x_b, y_b, tol=1e-6): + """Return an index array `perm` such that (x_b[perm], y_b[perm]) matches + (x_a, y_a) node-for-node. Both meshes must share the same node set + (here: same beam2d_tri.msh used on both sides).""" + order_a = np.lexsort((y_a, x_a)) + order_b = np.lexsort((y_b, x_b)) + if not (np.allclose(x_a[order_a], x_b[order_b], atol=tol) + and np.allclose(y_a[order_a], y_b[order_b], atol=tol)): + raise ValueError("Node coordinates don't match between SOFA and FreeFEM meshes.") + perm = np.empty_like(order_b) + perm[order_b] = order_a + return perm + + +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) + + q = float(cfg["q"]) + young_modulus = float(cfg["youngModulus"]) + poisson_ratio = float(cfg["poissonRatio"]) + mesh_file = _default_mesh_path() + + # --- Run FreeFEM --- + runner = FreeFemRunner("freefem_beam2d_distributed_cp.edp") + exports = runner.execute({ + 'q': q, + 'youngModulus': young_modulus, + 'poissonRatio': poisson_ratio, + 'meshFile': mesh_file, + }) + x_ff = exports['xcoords'] + y_ff = exports['ycoords'] + ux_ff = exports['ux[]'] + uy_ff = exports['uy[]'] + + # --- Run SOFA --- + pos0_sofa, u_sofa = sofaRun(mesh_file=mesh_file, q=q, + young_modulus=young_modulus, + poisson_ratio=poisson_ratio) + x_sofa, y_sofa = pos0_sofa[:, 0], pos0_sofa[:, 1] + ux_sofa, uy_sofa = u_sofa[:, 0], u_sofa[:, 1] + + + perm = _pair_by_coordinates(x_sofa, y_sofa, x_ff, y_ff) + ux_ff_p = ux_ff[perm] + uy_ff_p = uy_ff[perm] + + rms_ux = _rms(ux_sofa, ux_ff_p) + rms_uy = _rms(uy_sofa, uy_ff_p) + + # --- Write results --- + os.makedirs("results", exist_ok=True) + with open("results/comparison_beam2d_distributed_results.txt", 'w') as f: + header = f"{'x':>10} {'y':>10} {'ux_sofa':>12} {'ux_ff':>12} {'uy_sofa':>12} {'uy_ff':>12}" + f.write(header + "\n") + f.write("-" * len(header) + "\n") + for x, y, uxs, uxf, uys, uyf in zip(x_sofa, y_sofa, ux_sofa, ux_ff_p, uy_sofa, uy_ff_p): + f.write(f"{x:10.4f} {y:10.4f} {uxs:12.6e} {uxf:12.6e} {uys:12.6e} {uyf:12.6e}\n") + + f.write("\n") + f.write("RMS norms (SOFA vs FreeFEM, discrete, sqrt(n)-normalized)\n") + f.write("-" * 55 + "\n") + f.write(f" RMS_ux (SOFA vs FF) = {rms_ux:.6e}\n") + f.write(f" RMS_uy (SOFA vs FF) = {rms_uy:.6e}\n") + + print(f"RMS_ux (SOFA vs FF) = {rms_ux:.6e}") + print(f"RMS_uy (SOFA vs FF) = {rms_uy:.6e}") + + # --- Plots: SOFA / FreeFEM / difference, for ux and uy --- + _, triangles, _ = read_gmsh_2d(mesh_file) + triang = mtri.Triangulation(x_sofa, y_sofa, triangles) + + def _plot_field(ax, vals, title, cmap='viridis'): + tc = ax.tripcolor(triang, vals, shading='gouraud', cmap=cmap) + plt.colorbar(tc, ax=ax) + ax.set_title(title) + ax.set_aspect('equal') + + fig, axes = plt.subplots(2, 3, figsize=(16, 8)) + fig.suptitle("2D Beam — Distributed Load — SOFA vs FreeFEM (plane stress)", fontsize=14) + _plot_field(axes[0, 0], ux_sofa, 'SOFA : u_x') + _plot_field(axes[0, 1], ux_ff_p, 'FreeFEM : u_x') + _plot_field(axes[0, 2], ux_sofa - ux_ff_p, 'Diff u_x', 'RdBu') + _plot_field(axes[1, 0], uy_sofa, 'SOFA : u_y') + _plot_field(axes[1, 1], uy_ff_p, 'FreeFEM : u_y') + _plot_field(axes[1, 2], uy_sofa - uy_ff_p, 'Diff u_y', 'RdBu') + plt.tight_layout() + fig.savefig("results/comparison_beam2d_distributed_fields.png", dpi=150) + plt.close(fig) + \ No newline at end of file diff --git a/examples/Freefem/validation/2D/2D_planeStress/freefem_beam2d_distributed_cp.edp b/examples/Freefem/validation/2D/2D_planeStress/freefem_beam2d_distributed_cp.edp new file mode 100644 index 00000000..5c33e8fa --- /dev/null +++ b/examples/Freefem/validation/2D/2D_planeStress/freefem_beam2d_distributed_cp.edp @@ -0,0 +1,49 @@ +IMPORT "io.edp" +load "gmsh" + +DEFAULT (q, 1000.0) +DEFAULT (youngModulus, 1000.0) +DEFAULT (poissonRatio, 0.3) +DEFAULT (meshFile, "beam2d_tri.msh") + + +real Q = $q; +real E = $youngModulus; +real nu = $poissonRatio; + + +real mu = E / (2.*(1.+nu)); +real lambda = E*nu / (1. - nu^2); // plane stress + +// --- Mesh --- +mesh Th = gmshload("$meshFile"); + + +// --- FE space --- +fespace Vh(Th, P1); +Vh ux, uy, vx, vy; + + +solve Elasticity([ux, uy], [vx, vy]) = + int2d(Th)( + lambda*(dx(ux)+dy(uy))*(dx(vx)+dy(vy)) + + 2.*mu*( dx(ux)*dx(vx) + dy(uy)*dy(vy) ) + + mu*( dy(ux)+dx(uy) )*( dy(vx)+dx(vy) ) + ) + - int2d(Th)( -Q*vy ) + + on(1, ux=0, uy=0); + +Elasticity; + +exportArray(ux[]); +exportArray(uy[]); + + +real[int] xcoords(Th.nv), ycoords(Th.nv); +for (int i = 0; i < Th.nv; i++) { + xcoords[i] = Th(i).x; + ycoords[i] = Th(i).y; +} + +exportArray(xcoords); +exportArray(ycoords); diff --git a/examples/Freefem/validation/2D/2D_planeStress/params_beam2d_distributed.json b/examples/Freefem/validation/2D/2D_planeStress/params_beam2d_distributed.json new file mode 100644 index 00000000..817dc50d --- /dev/null +++ b/examples/Freefem/validation/2D/2D_planeStress/params_beam2d_distributed.json @@ -0,0 +1,7 @@ +{ + +"q": 1000, +"youngModulus": 1000, +"poissonRatio": 0.3 + +} \ No newline at end of file diff --git a/examples/Freefem/validation/2D/2D_planeStress/sofa_beam2d_distributed_cp.py b/examples/Freefem/validation/2D/2D_planeStress/sofa_beam2d_distributed_cp.py new file mode 100644 index 00000000..4d05ed85 --- /dev/null +++ b/examples/Freefem/validation/2D/2D_planeStress/sofa_beam2d_distributed_cp.py @@ -0,0 +1,222 @@ +""" +2D Beam Simulation - Distributed Load - Plane Stress - P1 Triangles +Cross-validation SOFA vs FreeFEM (no analytical solution available for this case). + +Physical labels in beam2d_tri.msh (Gmsh 2.2 format): + 1 = Fixed (x=0, boundary edges) -> Dirichlet + 2 = Bottom (y=0) + 3 = Right (x=L) + 4 = Top (y=H) + 5 = Beam (2D domain, triangles) + +Plane stress is obtained with the plugin's LinearSmallStrainFEMForceField by +using template="Vec2d" (genuine 2D elasticity), as opposed to the plane +strain case which uses template="Vec3d" with z held at 0 (see +sofa_beam2d_distributed.py). This mirrors the switch used in the +incompressibility MMS study (dim="2d" -> Vec2d -> plane stress lame lambda += E*nu/(1-nu**2); dim="3d" -> Vec3d -> plane strain lame lambda = +E*nu/((1+nu)*(1-2*nu))). +""" +import json +import os +import sys +import numpy as np +import Sofa +import Sofa.Core +import Sofa.Simulation + +RESULTS_DIR = "results" + + +def read_gmsh_2d(path): + """Minimal Gmsh 2.2 ASCII reader. + + Returns: + nodes : (N,2) array of x,y coordinates, in file node-id order + (node id 1 -> row 0, etc.) + triangles : (M,3) int array of 0-based node indices, region 5 (Beam) + fixed_idx : sorted list of 0-based node indices on boundary label 1 (Fixed) + """ + with open(path) as f: + lines = f.read().splitlines() + + i = 0 + nodes = None + triangles = [] + fixed_nodes = set() + + while i < len(lines): + line = lines[i].strip() + if line == "$Nodes": + n = int(lines[i + 1]) + nodes = np.zeros((n, 2)) + for k in range(n): + parts = lines[i + 2 + k].split() + node_id = int(parts[0]) + nodes[node_id - 1, 0] = float(parts[1]) + nodes[node_id - 1, 1] = float(parts[2]) + i += 2 + n + elif line == "$Elements": + m = int(lines[i + 1]) + for k in range(m): + parts = lines[i + 2 + k].split() + elm_type = int(parts[1]) + n_tags = int(parts[2]) + physical = int(parts[3]) + node_ids = [int(x) - 1 for x in parts[3 + n_tags:]] + if elm_type == 2: # 3-node triangle + if physical == 5: + triangles.append(node_ids) + elif elm_type == 1: # 2-node line + if physical == 1: # Fixed + fixed_nodes.update(node_ids) + i += 2 + m + else: + i += 1 + + triangles = np.array(triangles, dtype=int) + fixed_idx = sorted(fixed_nodes) + return nodes, triangles, fixed_idx + + +def consistent_nodal_forces(nodes, triangles, q): + """Consistent (lumped, 1/3 per node) nodal forces for a constant body + force (0, -q) per unit volume, integrated over each triangle.""" + N = len(nodes) + F = np.zeros((N, 2)) + for tri in triangles: + pts = nodes[tri, :] + v1 = pts[1] - pts[0] + v2 = pts[2] - pts[0] + area = 0.5 * abs(v1[0] * v2[1] - v1[1] * v2[0]) + for nid in tri: + F[nid, 1] += -q * area / 3.0 + return F + + +def _default_mesh_path(): + return os.path.join(os.path.dirname(os.path.abspath(__file__)), "beam2d_tri.msh") + + +def _default_params_path(): + return os.path.join(os.path.dirname(os.path.abspath(__file__)), "params_beam2d_distributed.json") + + +def create_scene_args(rootNode, mesh_file, q, young_modulus, poisson_ratio): + requiredPlugins = [ + "Elasticity", + "Sofa.Component.Constraint.Projective", + "Sofa.Component.LinearSolver.Direct", + "Sofa.Component.MechanicalLoad", + "Sofa.Component.ODESolver.Backward", + "Sofa.Component.StateContainer", + "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"]) + + nodes, triangles, fixed_idx = read_gmsh_2d(mesh_file) + N = len(nodes) + + # Plane stress -> template Vec2d (genuine 2D elasticity, no z dof) + template = "Vec2d" + positions = nodes.tolist() + tris_list = triangles.tolist() + + F_nodal = consistent_nodal_forces(nodes, triangles, q) + forces_list = F_nodal.tolist() + + with rootNode.addChild('Beam') as Beam: + Beam.addObject('NewtonRaphsonSolver' + , name="newtonSolver" + , printLog=True + , maxNbIterationsNewton=30 + , absoluteResidualStoppingThreshold=1e-12) + Beam.addObject('SparseLDLSolver' + , name="linearSolver" + , template="CompressedRowSparseMatrixd") + Beam.addObject('StaticSolver' + , name="staticSolver" + , newtonSolver="@newtonSolver" + , linearSolver="@linearSolver") + + dofs = Beam.addObject('MechanicalObject' + , name="dofs" + , template=template + , position=positions + , showObject=True + , showObjectScale=0.01) + + Beam.addObject('TriangleSetTopologyContainer' + , name="topology" + , position="@dofs.position" + , triangles=tris_list) + Beam.addObject('TriangleSetTopologyModifier') + + Beam.addObject('LinearSmallStrainFEMForceField' + , name="FEM" + , template=template + , youngModulus=young_modulus + , poissonRatio=poisson_ratio + , topology="@topology") + + Beam.addObject('FixedProjectiveConstraint' + , name="dirichlet" + , indices=fixed_idx) + + Beam.addObject('ConstantForceField' + , name="DistributedLoad" + , indices=list(range(N)) + , forces=forces_list) + + return rootNode, dofs, np.array(positions) + + +def createScene(rootNode): + with open(_default_params_path()) as f: + cfg = json.load(f) + create_scene_args(rootNode + , mesh_file=_default_mesh_path() + , q=float(cfg["q"]) + , young_modulus=float(cfg["youngModulus"]) + , poisson_ratio=float(cfg["poissonRatio"])) + return rootNode + + +def sofaRun(mesh_file, q, young_modulus, poisson_ratio): + root = Sofa.Core.Node("root") + _, dofs, pos0 = create_scene_args(root + , mesh_file=mesh_file + , q=q + , young_modulus=young_modulus + , poisson_ratio=poisson_ratio) + Sofa.Simulation.init(root) + Sofa.Simulation.animate(root, root.dt.value) + + pos_final = np.array(dofs.position.toList()) + u = pos_final - pos0 # (N,2): ux, uy directly, no z dof to drop + x0 = pos0 + + os.makedirs(RESULTS_DIR, exist_ok=True) + out_path = os.path.join(RESULTS_DIR, "sofa_beam2d_distributed_cp_results.txt") + with open(out_path, 'w') as f: + f.write(f"{'x0':>12} {'y0':>12} {'ux':>12} {'uy':>12}\n") + f.write("-" * 54 + "\n") + for (xi, yi), (uxi, uyi) in zip(x0, u): + f.write(f"{xi:12.6f} {yi:12.6f} {uxi:12.6f} {uyi:12.6f}\n") + + return x0, u + + +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(mesh_file=_default_mesh_path() + , q=float(cfg["q"]) + , young_modulus=float(cfg["youngModulus"]) + , poisson_ratio=float(cfg["poissonRatio"])) \ No newline at end of file