Add 2D beam simulation: SOFA and FreeFEM files

examples/Freefem/2D-beam/compare_sofa_freefem.py

@@ -0,0 +1,159 @@
+"""
+2D Beam Simulation - Deformed Under Gravity - Comparison File
+"""
+import json
+import os
+import sys
+import numpy as np
+import matplotlib.pyplot as plt
+import matplotlib.tri as mtri
+from gmsh_generate_beam2D import generate_beam2D
+from sofa_beam2d import sofaRun
+from pyfreefem import FreeFemRunner
+
+
+def _l2(a, b):
+    return np.linalg.norm(a - b) / a.size
+
+
+def _write_freefem_results(path, x, y, u_x, u_y):
+    x_final = x + u_x
+    y_final = y + u_y
+    with open(path, 'w') as f:
+        f.write(f"{'x_initial':>12} {'y_initial':>12} {'x_final':>12} {'y_final':>12} {'u_x':>12} {'u_y':>12}\n")
+        f.write("-" * 80 + "\n")
+        for xi, yi, xf, yf, uxi, uyi in zip(x, y, x_final, y_final, u_x, u_y):
+            f.write(f"{xi:12.6f} {yi:12.6f} {xf:12.6f} {yf:12.6f} {uxi:12.6f} {uyi:12.6f}\n")
+
+
+if __name__ == "__main__":
+
+    ## --- Read parameters from JSON file ---
+    config_file = sys.argv[1] if len(sys.argv) > 1 else "params.json"
+    with open(config_file) as f:
+        cfg = json.load(f)
+
+    length          = float(cfg["length"])
+    height          = float(cfg["height"])
+    nx              = int(cfg["nx"])
+    ny              = int(cfg["ny"])
+    force           = float(cfg["force"])
+    rho             = float(cfg["rho"])
+    gravity         = float(cfg["gravity"])
+    young_modulus   = float(cfg["youngModulus"])
+    poisson_ratio   = float(cfg["poissonRatio"])
+    mesh_filename   = cfg["meshfile"]
+
+    # --- Run Gmsh ---
+    path_before_gmsh = os.environ.get('PATH', '')
+    msh_path, x_positions, y_positions = generate_beam2D(
+        length=length,
+        height=height,
+        nx=nx,
+        ny=ny,
+        filename=mesh_filename
+    )
+
+    # --- Run FreeFEM ---
+    # gmsh.initialize() corrupts PATH, breaking pyfreefem's stdbuf call
+    os.environ['PATH'] = path_before_gmsh
+    runner = FreeFemRunner("freefem_beam2d.edp")
+    results = runner.execute({
+        'youngModulus': young_modulus,
+        'poissonRatio': poisson_ratio,
+        'rhoMat': rho,
+        'gravity': gravity,
+        'meshFile': os.path.abspath(msh_path).replace('\\', '/'),
+    })
+    x_ff = results['xcoords']
+    y_ff = results['ycoords']
+    dofs = results['uu[]']  # interleaved [uu_0, vv_0, uu_1, vv_1, ...]
+    u_x_ff = dofs[0::2]
+    u_y_ff = dofs[1::2]
+
+    os.makedirs("results", exist_ok=True)
+    _write_freefem_results(os.path.join("results", "freefem_results.txt"), x_ff, y_ff, u_x_ff, u_y_ff)
+
+    ## --- Run SOFA ---
+    pos_initial, u_sofa = sofaRun(
+        height=height,
+        young_modulus=young_modulus,
+        poisson_ratio=poisson_ratio,
+        rho=rho,
+        gravity=gravity,
+        mesh_file=os.path.abspath(msh_path).replace('\\', '/')
+    )
+    x_sofa   = pos_initial[:, 0]
+    y_sofa   = pos_initial[:, 1]
+    u_x_sofa = u_sofa[:, 0]
+    u_y_sofa = u_sofa[:, 1]
+
+    ## --- Compare Results ---
+    err_ux = u_x_sofa - u_x_ff
+    err_uy = u_y_sofa - u_y_ff
+
+    with open("results/comparison_results.txt", 'w') as f:
+        header = (f"{'x':>12}  {'y':>12}  {'ux_ff':>12}  {'ux_sofa':>12}"
+                  f"  {'uy_ff':>12}  {'uy_sofa':>12}  {'err_ux':>12}  {'err_uy':>12}")
+        f.write(header + "\n")
+        f.write("-" * len(header) + "\n")
+        for vals in zip(x_ff, y_ff, u_x_ff, u_x_sofa, u_y_ff, u_y_sofa, err_ux, err_uy):
+            f.write("  ".join(f"{v:12.6f}" for v in vals) + "\n")
+        f.write("\n")
+        f.write("L2 error norms\n")
+        f.write("-" * 40 + "\n")
+        f.write(f"  ||ux_sofa - ux_ff||_2 = {_l2(u_x_sofa, u_x_ff):.6e}\n")
+        f.write(f"  ||uy_sofa - uy_ff||_2 = {_l2(u_y_sofa, u_y_ff):.6e}\n")
+
+    ## --- Plots ---
+    triang = mtri.Triangulation(x_ff, y_ff)
+
+    # Shared colour ranges so SOFA and FreeFEM panels are directly comparable
+    vmin_ux = min(u_x_ff.min(), u_x_sofa.min())
+    vmax_ux = max(u_x_ff.max(), u_x_sofa.max())
+    vmin_uy = min(u_y_ff.min(), u_y_sofa.min())
+    vmax_uy = max(u_y_ff.max(), u_y_sofa.max())
+
+    # Fig 1 — displacement field: FreeFEM (top row) vs SOFA (bottom row)
+    fig, axes = plt.subplots(2, 2, figsize=(14, 6), constrained_layout=True)
+    fig.suptitle("Displacement field — FreeFEM (top) vs SOFA (bottom)")
+    for row, (label, ux, uy) in enumerate([
+        ("FreeFEM", u_x_ff,   u_y_ff),
+        ("SOFA",    u_x_sofa, u_y_sofa),
+    ]):
+        tcf = axes[row, 0].tricontourf(triang, ux, levels=20, vmin=vmin_ux, vmax=vmax_ux)
+        axes[row, 0].set_title(f"{label} — $u_x$")
+        axes[row, 0].set_aspect('equal')
+        fig.colorbar(tcf, ax=axes[row, 0])
+
+        tcf = axes[row, 1].tricontourf(triang, uy, levels=20, vmin=vmin_uy, vmax=vmax_uy)
+        axes[row, 1].set_title(f"{label} — $u_y$")
+        axes[row, 1].set_aspect('equal')
+        fig.colorbar(tcf, ax=axes[row, 1])
+    fig.savefig("results/plot_displacement_field.png", dpi=150)
+    plt.close(fig)
+
+    # Fig 2 — deformed beams overlaid
+    fig, ax = plt.subplots(figsize=(14, 4))
+    ax.scatter(x_ff   + u_x_ff,   y_ff   + u_y_ff,   s=4, label="FreeFEM", alpha=0.7)
+    ax.scatter(x_sofa + u_x_sofa, y_sofa + u_y_sofa, s=4, label="SOFA",    alpha=0.7, marker='x')
+    ax.set_aspect('equal')
+    ax.set_title("Deformed beam — FreeFEM vs SOFA")
+    ax.set_xlabel("x")
+    ax.set_ylabel("y")
+    ax.legend()
+    fig.savefig("results/plot_deformed_beams.png", dpi=150)
+    plt.close(fig)
+
+    # Figs 3 & 4 — nodal error in x and y
+    for comp, err in [("ux", err_ux), ("uy", err_uy)]:
+        fig, ax = plt.subplots(figsize=(14, 4))
+        tcf = ax.tricontourf(triang, err, levels=20, cmap='RdBu_r')
+        fig.colorbar(tcf, ax=ax, label=f"$u_{comp[1]}$ error (SOFA − FreeFEM)")
+        ax.set_title(f"Nodal error — $u_{comp[1]}$")
+        ax.set_aspect('equal')
+        ax.set_xlabel("x")
+        ax.set_ylabel("y")
+        fig.savefig(f"results/plot_error_{comp}.png", dpi=150)
+        plt.close(fig)
+

examples/Freefem/2D-beam/freefem_beam2d.edp

@@ -0,0 +1,46 @@
+IMPORT "io.edp"
+load "gmsh"
+
+DEFAULT (youngModulus, "210000")
+DEFAULT (rhoMat, "7.85e-6")
+DEFAULT (gravity, "9810")
+DEFAULT (poissonRatio, "0.3")
+DEFAULT (meshFile, "")
+
+real E       = $youngModulus;
+real sigma   = $poissonRatio;
+real rhoMat  = $rhoMat;
+real gravity = $gravity;
+real F       = -rhoMat * gravity;
+
+mesh th = gmshload("$meshFile");
+
+fespace Vh(th, [P1, P1]);
+Vh [uu, vv], [w, s];
+
+real sqrt2  = sqrt(2.);
+macro epsilon(u1,u2) [dx(u1), dy(u2), (dy(u1)+dx(u2))/sqrt2] // EOM
+macro div(u,v) (dx(u)+dy(v)) // EOM
+
+real mu     = E / (2*(1 + sigma));
+real lambda = E*sigma / ((1 + sigma)*(1 - 2*sigma));
+
+solve Elasticity([uu, vv], [w, s], solver=sparsesolver)
+    = int2d(th)(
+          lambda*div(w,s)*div(uu,vv)
+        + 2.*mu*(epsilon(w,s)' * epsilon(uu,vv))
+    )
+    - int2d(th)(F*s)
+    + on(1, uu=0, vv=0);
+
+exportArray(uu[]);
+exportArray(vv[]);
+
+real[int] xcoords(th.nv);
+real[int] 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);

examples/Freefem/2D-beam/gmsh_generate_beam2D.py

@@ -0,0 +1,75 @@
+import json
+import os
+import sys
+import gmsh
+
+RESULTS_DIR = "results"
+
+def generate_beam2D(length, height, nx, ny, filename):
+    """
+    Generate a 2D beam mesh (triangular elements)
+    """
+    gmsh.initialize()
+    gmsh.model.add("beam2d")
+
+    # Corner points (tag order matches border traversal in the .edp)
+    gmsh.model.geo.addPoint(0,      0,      0, tag=1)
+    gmsh.model.geo.addPoint(length, 0,      0, tag=2)
+    gmsh.model.geo.addPoint(length, height, 0, tag=3)
+    gmsh.model.geo.addPoint(0,      height, 0, tag=4)
+
+    # Boundary lines — tags chosen to match .edp labels directly
+    bottom = gmsh.model.geo.addLine(1, 2, tag=2)  # label 2, y=0
+    right  = gmsh.model.geo.addLine(2, 3, tag=3)  # label 3, x=length
+    top    = gmsh.model.geo.addLine(3, 4, tag=4)  # label 4, y=height
+    left   = gmsh.model.geo.addLine(4, 1, tag=1)  # label 1, x=0 (clamped)
+
+    loop = gmsh.model.geo.addCurveLoop([bottom, right, top, left], tag=1)
+    surf = gmsh.model.geo.addPlaneSurface([loop], tag=1)
+
+    gmsh.model.geo.synchronize()
+
+    # Physical groups — boundary tags match .edp labels 1–4
+    gmsh.model.addPhysicalGroup(1, [left],   tag=1, name="Fixed")
+    gmsh.model.addPhysicalGroup(1, [bottom], tag=2, name="Bottom")
+    gmsh.model.addPhysicalGroup(1, [right],  tag=3, name="Right")
+    gmsh.model.addPhysicalGroup(1, [top],    tag=4, name="Top")
+    gmsh.model.addPhysicalGroup(2, [surf],   tag=5, name="Beam")
+
+    # Structured transfinite layout; no recombine → triangular elements
+    gmsh.model.mesh.setTransfiniteCurve(bottom, nx)
+    gmsh.model.mesh.setTransfiniteCurve(top,    nx)
+    gmsh.model.mesh.setTransfiniteCurve(left,   ny)
+    gmsh.model.mesh.setTransfiniteCurve(right,  ny)
+    gmsh.model.mesh.setTransfiniteSurface(surf)
+
+    gmsh.model.mesh.generate(2)
+    gmsh.model.mesh.setOrder(1)  # ensure linear (first-order) elements only
+    _, node_coords, _ = gmsh.model.mesh.getNodes()
+    x_positions = node_coords[0::3]
+    y_positions = node_coords[1::3]
+
+    os.makedirs(RESULTS_DIR, exist_ok=True)
+    msh_path = os.path.join(RESULTS_DIR, filename)
+    gmsh.option.setNumber("Mesh.MshFileVersion", 2.2)  # gmshload requires MSH v2
+    gmsh.write(msh_path)
+    gmsh.finalize()
+
+    return msh_path, x_positions, y_positions
+
+
+if __name__ == "__main__":
+    config_file = sys.argv[1] if len(sys.argv) > 1 else "params.json"
+    with open(config_file) as f:
+        cfg = json.load(f)
+
+    msh_path, x_positions, y_positions = generate_beam2D(
+        length=float(cfg["length"]),
+        height=float(cfg["height"]),
+        nx=int(cfg["nx"]),
+        ny=int(cfg["ny"]),
+        filename=cfg.get("meshfile"),
+    )
+    print("Wrote:", msh_path)
+    print("Node x-positions (Gmsh order):", x_positions)
+    print("Node y-positions (Gmsh order):", y_positions)

examples/Freefem/2D-beam/params.json

@@ -0,0 +1,12 @@
+{
+    "length": 1.0,
+    "height" : 0.2,
+    "nx": 20,
+    "ny": 8,
+    "force": 1000.0,
+    "rho": 20,
+    "gravity": 9.8,
+    "youngModulus": 210000,
+    "poissonRatio": 0.3,
+    "meshfile": "beam2d.msh"
+}