Add 3D circular beam MMS validation: SOFA vs FreeFEM vs Saint-Venant analytical

examples/Freefem/validation/3D/3D_circulaire/beam3d_circular_tet.py

@@ -0,0 +1,70 @@
+import json
+import os
+import sys
+import gmsh
+
+RESULTS_DIR = "results"
+
+
+def generate_beam3D_circular_tet(length, radius, mesh_size, filename):
+
+    gmsh.initialize()
+    gmsh.model.add("beam3d_circular_tet")
+
+    disk_fixed = gmsh.model.occ.addDisk(0, 0, 0, radius, radius,
+                                         zAxis=[1, 0, 0], xAxis=[0, 1, 0])
+    gmsh.model.occ.synchronize()
+
+
+    out = gmsh.model.occ.extrude([(2, disk_fixed)], length, 0, 0)
+    gmsh.model.occ.synchronize()
+
+    vol_tag     = [e[1] for e in out if e[0] == 3][0]
+    surf_tags   = [e[1] for e in out if e[0] == 2]
+
+    disk_loaded = None
+    lateral = None
+    for s in surf_tags:
+        xmin, ymin, zmin, xmax, ymax, zmax = gmsh.model.occ.getBoundingBox(2, s)
+        if abs(xmax - xmin) < 1e-6:
+            disk_loaded = s
+        else:
+            lateral = s
+
+    assert disk_loaded is not None and lateral is not None, \
+        "Could not identify end cap"
+
+    gmsh.model.addPhysicalGroup(2, [disk_fixed],  tag=1, name="Fixed")
+    gmsh.model.addPhysicalGroup(2, [disk_loaded], tag=2, name="Loaded")
+    gmsh.model.addPhysicalGroup(2, [lateral],     tag=3, name="Lateral")
+    gmsh.model.addPhysicalGroup(3, [vol_tag],     tag=4, name="Beam")
+
+    gmsh.model.mesh.setSize(gmsh.model.getEntities(0), mesh_size)
+
+    gmsh.model.mesh.generate(3)
+    gmsh.model.mesh.setOrder(1) 
+    _, node_coords, _ = gmsh.model.mesh.getNodes()
+
+    os.makedirs(RESULTS_DIR, exist_ok=True)
+    msh_path = os.path.join(RESULTS_DIR, filename)
+    gmsh.option.setNumber("Mesh.MshFileVersion", 2.2) 
+    gmsh.write(msh_path)
+    gmsh.finalize()
+
+    return msh_path, len(node_coords) // 3
+
+
+if __name__ == "__main__":
+    config_file = sys.argv[1] if len(sys.argv) > 1 else "params.json"
+    with open(config_file) as f:
+        all_cfg = json.load(f)
+    cfg = all_cfg["beam3d_circle_tet"]
+
+    msh_path, n_nodes = generate_beam3D_circular_tet(
+        length=float(cfg["length"]),
+        radius=float(cfg["radius"]),
+        mesh_size=float(cfg["mesh_size"]),
+        filename=cfg.get("meshfile"),
+    )
+    print("Wrote:", msh_path)
+    print("Number of nodes:", n_nodes)

examples/Freefem/validation/3D/3D_circulaire/comparaison_script3d_circle.py

@@ -0,0 +1,154 @@
+# =============================================================================
+# comparaison_script3d_circle.py
+# 3D circular beam, axial traction: SOFA vs FreeFEM vs analytical (Saint-Venant)
+#
+# Same structure as comparaison_script3d.py (distributed-load case): generate
+# the mesh, run FreeFEM, run SOFA, match nodes by coordinates, compute RMS.
+# Here we additionally compare both solvers against the exact analytical
+# solution.
+# =============================================================================
+
+import os
+import json
+import numpy as np
+import matplotlib.pyplot as plt
+
+from pyfreefem import FreeFemRunner
+import sofa_beam3d_circle_traction as sofa_case
+
+RESULTS_DIR = "results"
+
+
+def _default_edp_path():
+    return os.path.join(os.path.dirname(os.path.abspath(__file__)),
+                         "freefem_beam3d_circle_traction.edp")
+
+
+def match_by_coordinates(coords_a, vals_a, coords_b, vals_b, tol=1e-6): 
+    matched_a, matched_b, matched_coords = [], [], []
+    for i, c in enumerate(coords_a):
+        d = np.linalg.norm(coords_b - c, axis=1)
+        j = np.argmin(d)
+        if d[j] > tol:
+            continue
+        matched_coords.append(c)
+        matched_a.append(vals_a[i])
+        matched_b.append(vals_b[j])
+    return np.array(matched_coords), np.array(matched_a), np.array(matched_b)
+
+
+def analytical_solution(coords, E, nu, q):
+    x, y, z = coords[:, 0], coords[:, 1], coords[:, 2]
+    ux = (q / E) * x
+    uy = -nu * (q / E) * y
+    uz = -nu * (q / E) * z
+    return ux, uy, uz
+
+
+def rms(a, b):
+    return float(np.sqrt(np.mean((a - b) ** 2)))
+
+
+def main():
+    with open("params.json") as f:
+        mesh_cfg = json.load(f)["beam3d_circle_tet"]
+    with open("params_beam3d_circle_traction.json") as f:
+        phys_cfg = json.load(f)
+
+    length = mesh_cfg["length"]
+    radius = mesh_cfg["radius"]
+    E = phys_cfg["youngModulus"]
+    nu = phys_cfg["poissonRatio"]
+    q = phys_cfg["q"]
+    exclusion_factor = phys_cfg.get("exclusionFactor", 2.0)
+
+    msh_path = os.path.join(RESULTS_DIR, mesh_cfg["meshfile"])
+    if not os.path.exists(msh_path):
+        raise FileNotFoundError(
+            f"{msh_path} not found -- run beam3d_circular_tet.py first "
+            f"(and check_mesh_conformity.py to validate it)."
+        )
+    print(f"Mesh: {msh_path}")
+
+    runner = FreeFemRunner(_default_edp_path())
+    exports = runner.execute({
+        "meshfile": os.path.abspath(msh_path),
+        "E": E, "nu": nu, "q": q, "radius": radius, "length": length,
+        "exclusionFactor": exclusion_factor,
+    })
+    ux_ff = exports["ux[]"] if "ux[]" in exports else exports["ux"]
+    uy_ff = exports["uy[]"] if "uy[]" in exports else exports["uy"]
+    uz_ff = exports["uz[]"] if "uz[]" in exports else exports["uz"]
+    xcoords = exports["xcoords"]
+    ycoords = exports["ycoords"]
+    zcoords = exports["zcoords"]
+    coords_ff = np.column_stack([xcoords, ycoords, zcoords])
+
+    coords_sofa, u_sofa = sofa_case.sofaRun(
+        mesh_file=msh_path, q=q, young_modulus=E, poisson_ratio=nu,
+    )
+    ux_sofa, uy_sofa, uz_sofa = u_sofa[:, 0], u_sofa[:, 1], u_sofa[:, 2] 
+
+    tol = 1e-6 * max(length, radius)
+    coords_m, ux_sofa_m, ux_ff_m = match_by_coordinates(
+        coords_sofa, ux_sofa, coords_ff, ux_ff, tol=tol)
+    _, uy_sofa_m, uy_ff_m = match_by_coordinates(
+        coords_sofa, uy_sofa, coords_ff, uy_ff, tol=tol)
+    _, uz_sofa_m, uz_ff_m = match_by_coordinates(
+        coords_sofa, uz_sofa, coords_ff, uz_ff, tol=tol)
+
+    rms_ux = rms(ux_sofa_m, ux_ff_m)
+    rms_uy = rms(uy_sofa_m, uy_ff_m)
+    rms_uz = rms(uz_sofa_m, uz_ff_m)
+
+    print(f"RMS_ux (SOFA vs FF) = {rms_ux:.6e}")
+    print(f"RMS_uy (SOFA vs FF) = {rms_uy:.6e}")
+    print(f"RMS_uz (SOFA vs FF) = {rms_uz:.6e}")
+
+    exclusion_x = exclusion_factor * radius
+    far_mask = coords_m[:, 0] >= exclusion_x
+    ux_an, uy_an, uz_an = analytical_solution(coords_m[far_mask], E, nu, q)
+
+    rms_ux_sofa_an = rms(ux_sofa_m[far_mask], ux_an)
+    rms_uy_sofa_an = rms(uy_sofa_m[far_mask], uy_an)
+    rms_uz_sofa_an = rms(uz_sofa_m[far_mask], uz_an)
+    rms_ux_ff_an   = rms(ux_ff_m[far_mask], ux_an)
+    rms_uy_ff_an   = rms(uy_ff_m[far_mask], uy_an)
+    rms_uz_ff_an   = rms(uz_ff_m[far_mask], uz_an)
+
+
+    os.makedirs(RESULTS_DIR, exist_ok=True)
+    out_path = os.path.join(RESULTS_DIR, "comparison_beam3d_circle_traction_results.txt")
+    with open(out_path, "w") as f:
+        f.write("x y z ux_sofa ux_ff uy_sofa uy_ff uz_sofa uz_ff\n")
+        f.write("-" * 100 + "\n")
+        for i in range(len(coords_m)):
+            f.write(f"{coords_m[i,0]:10.4f} {coords_m[i,1]:10.4f} {coords_m[i,2]:10.4f} "
+                     f"{ux_sofa_m[i]:14.6e} {ux_ff_m[i]:14.6e} "
+                     f"{uy_sofa_m[i]:14.6e} {uy_ff_m[i]:14.6e} "
+                     f"{uz_sofa_m[i]:14.6e} {uz_ff_m[i]:14.6e}\n")
+        f.write("\nRMS norms (SOFA vs FreeFEM)\n")
+        f.write(f"  RMS_ux = {rms_ux:.6e}\n  RMS_uy = {rms_uy:.6e}\n  RMS_uz = {rms_uz:.6e}\n")
+        f.write(f"\nRMS norms vs analytical (x >= {exclusion_x})\n")
+        f.write(f"  SOFA: ux={rms_ux_sofa_an:.6e} uy={rms_uy_sofa_an:.6e} uz={rms_uz_sofa_an:.6e}\n")
+        f.write(f"  FF:   ux={rms_ux_ff_an:.6e} uy={rms_uy_ff_an:.6e} uz={rms_uz_ff_an:.6e}\n")
+    print("\nWrote:", out_path)
+
+    fig, axes = plt.subplots(1, 3, figsize=(15, 5))
+    for ax, (a, b, name) in zip(axes, [(ux_sofa_m, ux_ff_m, "ux"),
+                                        (uy_sofa_m, uy_ff_m, "uy"),
+                                        (uz_sofa_m, uz_ff_m, "uz")]):
+        ax.scatter(a, b, s=15, alpha=0.6)
+        lims = [min(a.min(), b.min()), max(a.max(), b.max())]
+        ax.plot(lims, lims, "r--", linewidth=1)
+        ax.set_xlabel(f"{name}_sofa")
+        ax.set_ylabel(f"{name}_ff")
+        ax.set_title(name)
+    fig.suptitle("3D Circular Beam — Axial Traction — SOFA vs FreeFEM (parity)")
+    fig_path = os.path.join(RESULTS_DIR, "comparison_beam3d_circle_traction_fields.png")
+    fig.savefig(fig_path, dpi=120, bbox_inches="tight")
+    print("Wrote:", fig_path)
+
+
+if __name__ == "__main__":
+    main()

examples/Freefem/validation/3D/3D_circulaire/freefem_beam3d_circle_traction.edp

@@ -0,0 +1,75 @@
+IMPORT "io.edp"
+load "msh3"
+load "gmsh"
+
+DEFAULT (meshfile, "beam3d_circle_tet.msh")
+DEFAULT (E, 210000.0)
+DEFAULT (nu, 0.3)
+DEFAULT (F, 1000.0)
+DEFAULT (radius, 0.1)
+DEFAULT (length, 1.0)
+DEFAULT (exclusionFactor, 2.0)
+
+real E      = $E;
+real nu     = $nu;
+real F      = $F;
+real radius = $radius;
+real length = $length;
+real exclusionFactor = $exclusionFactor;
+
+real A = pi*radius^2;
+real q = F/A;             
+
+real lambda = E*nu/((1.+nu)*(1.-2.*nu));
+real mu     = E/(2.*(1.+nu));
+
+
+mesh3 Th = gmshload3("$meshfile");
+
+
+fespace Vh(Th, P1);
+Vh ux, uy, uz, vx, vy, vz;
+
+solve Elasticity([ux, uy, uz], [vx, vy, vz]) =
+    int3d(Th)(
+        lambda*(dx(ux)+dy(uy)+dz(uz))*(dx(vx)+dy(vy)+dz(vz))
+      + mu*( 2.*dx(ux)*dx(vx) + 2.*dy(uy)*dy(vy) + 2.*dz(uz)*dz(vz)
+           + (dy(ux)+dx(uy))*(dy(vx)+dx(vy))
+           + (dz(ux)+dx(uz))*(dz(vx)+dx(vz))
+           + (dz(uy)+dy(uz))*(dz(vy)+dy(vz)) )
+    )
+  - int2d(Th, 2)( q*vx )
+  + on(1, ux=0, uy=0, uz=0);
+
+exportArray(ux[]);
+exportArray(uy[]);
+exportArray(uz[]);
+
+real[int] xcoords(Th.nv);
+real[int] ycoords(Th.nv);
+real[int] zcoords(Th.nv);
+for (int i = 0; i < Th.nv; i++) {
+    xcoords[i] = Th(i).x;
+    ycoords[i] = Th(i).y;
+    zcoords[i] = Th(i).z;
+}
+exportArray(xcoords);
+exportArray(ycoords);
+exportArray(zcoords);
+
+real exclusionx = exclusionFactor*radius;
+
+real sumSqUx = 0., sumSqUy = 0., sumSqUz = 0.;
+int nCompared = 0;
+for (int i = 0; i < Th.nv; i++) {
+    real xi = Th(i).x, yi = Th(i).y, zi = Th(i).z;
+    if (xi < exclusionx) continue;
+    real uxAn = (F/(E*A))*xi;
+    real uyAn = -nu*(F/(E*A))*yi;
+    real uzAn = -nu*(F/(E*A))*zi;
+    sumSqUx += (ux[](i) - uxAn)^2;
+    sumSqUy += (uy[](i) - uyAn)^2;
+    sumSqUz += (uz[](i) - uzAn)^2;
+    nCompared++;
+}
+

examples/Freefem/validation/3D/3D_circulaire/params.json

@@ -0,0 +1,8 @@
+{
+  "beam3d_circle_tet": {
+  "length": 1.0,
+  "radius": 0.1,
+  "mesh_size": 0.05,
+  "meshfile": "beam3d_circle_tet.msh"
+}
+}

examples/Freefem/validation/3D/3D_circulaire/params_beam3d_circle_traction.json

@@ -0,0 +1,8 @@
+{
+  "length": 1.0,
+  "radius": 0.1,
+  "youngModulus": 210000.0,
+  "poissonRatio": 0.3,
+  "q": 31830.988618,
+  "exclusionFactor": 2.0
+}