From 181fe08f559464df62481630bdfc72ca04b8aabd Mon Sep 17 00:00:00 2001 From: fatima Date: Thu, 30 Jul 2026 12:25:52 +0200 Subject: [PATCH] Add 3D beam torsion validation (circular section): mesh generation, SOFA vs FreeFEM vs Saint-Venant analytical comparison --- .../3D_torsion/compare_beam3d_torsion.py | 204 +++++++++++++++ .../3D_torsion/freefem_beam3d_torsion.edp | 99 +++++++ .../generate_beam3d_circular_tet.py | 72 +++++ .../Freefem/validation/3D_torsion/params.json | 8 + .../3D_torsion/params_beam3d_torsion.json | 6 + .../3D_torsion/sofa_beam3d_torsion.py | 247 ++++++++++++++++++ 6 files changed, 636 insertions(+) create mode 100644 examples/Freefem/validation/3D_torsion/compare_beam3d_torsion.py create mode 100644 examples/Freefem/validation/3D_torsion/freefem_beam3d_torsion.edp create mode 100644 examples/Freefem/validation/3D_torsion/generate_beam3d_circular_tet.py create mode 100644 examples/Freefem/validation/3D_torsion/params.json create mode 100644 examples/Freefem/validation/3D_torsion/params_beam3d_torsion.json create mode 100644 examples/Freefem/validation/3D_torsion/sofa_beam3d_torsion.py diff --git a/examples/Freefem/validation/3D_torsion/compare_beam3d_torsion.py b/examples/Freefem/validation/3D_torsion/compare_beam3d_torsion.py new file mode 100644 index 00000000..670c894e --- /dev/null +++ b/examples/Freefem/validation/3D_torsion/compare_beam3d_torsion.py @@ -0,0 +1,204 @@ +import json +import os +import sys +import numpy as np +import matplotlib.pyplot as plt + +from sofa_beam3d_torsion import sofaRun, MESH_DIR, DEFAULT_MESH_FILENAME +from pyfreefem import FreeFemRunner + + +def _rms(a, b): + return np.linalg.norm(a - b) / np.sqrt(a.size) + + +def _rel_rms(u_ref, u_test): + denom = np.linalg.norm(u_ref) + return float(np.linalg.norm(u_test - u_ref) / denom) if denom > 0 else float("nan") + + +def _default_params_path(): + return os.path.join(os.path.dirname(os.path.abspath(__file__)), "params_beam3d_torsion.json") + + +def _default_mesh_path(): + return os.path.join(MESH_DIR, DEFAULT_MESH_FILENAME) + + +def _to_freefem_path(path): + return path.replace(os.sep, "/") + + +def _rewrite_freefem_output_in_place(raw, out_path): + header = (f"{'x':>12} {'y':>12} {'z':>12} " + f"{'ux':>14} {'uy':>14} {'uz':>14}") + with open(out_path, 'w') as f: + f.write(header + "\n") + f.write("-" * len(header) + "\n") + for x, y, z, ux, uy, uz in raw: + f.write(f"{x:12.6f} {y:12.6f} {z:12.6f} " + f"{ux:+14.6e} {uy:+14.6e} {uz:+14.6e}\n") + + +def _pair_by_coordinates(x_a, y_a, z_a, x_b, y_b, z_b, tol=1e-6, snap=1e-6): + + x_a, y_a, z_a = map(np.asarray, (x_a, y_a, z_a)) + x_b, y_b, z_b = map(np.asarray, (x_b, y_b, z_b)) + + print(f"[diag] n_sofa={x_a.size} n_freefem={x_b.size}") + + if x_a.size != x_b.size: + raise ValueError( + f"Node COUNT mismatch: SOFA has {x_a.size} nodes, " + f"FreeFEM has {x_b.size} nodes. " + ) + + def snap_(v): + return np.round(v / snap) * snap + + xs_a, ys_a, zs_a = snap_(x_a), snap_(y_a), snap_(z_a) + xs_b, ys_b, zs_b = snap_(x_b), snap_(y_b), snap_(z_b) + + order_a = np.lexsort((zs_a, ys_a, xs_a)) + order_b = np.lexsort((zs_b, ys_b, xs_b)) + + da = np.stack([x_a[order_a], y_a[order_a], z_a[order_a]], axis=1) + db = np.stack([x_b[order_b], y_b[order_b], z_b[order_b]], axis=1) + diff = np.linalg.norm(da - db, axis=1) + + print(f"[diag] max sorted-coordinate discrepancy = {diff.max():.6e} (tol={tol:.1e})") + if diff.max() >= 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 + + +def _analytical_displacement(x0, torque, radius, young_modulus, poisson_ratio, yc, zc): + + G = young_modulus / (2.0 * (1.0 + poisson_ratio)) + J = np.pi * radius ** 4 / 2.0 + theta_prime = torque / (G * J) + + x, y, z = x0[:, 0], x0[:, 1], x0[:, 2] + ux = np.zeros_like(x) + uy = -theta_prime * x * (z - zc) + uz = theta_prime * x * (y - yc) + return np.column_stack([ux, uy, uz]), theta_prime + + +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) + + T = float(cfg["T"]) + radius = float(cfg["radius"]) + young_modulus = float(cfg["youngModulus"]) + poisson_ratio = float(cfg["poissonRatio"]) + mesh_file = _default_mesh_path() + + os.makedirs("results", exist_ok=True) + ff_out_path = os.path.join(os.path.dirname(os.path.abspath(__file__)), + "results", "freefem_beam3d_torsion_raw.txt") + + runner = FreeFemRunner("freefem_beam3d_torsion.edp") + runner.execute({ + 'T': T, + 'radius': radius, + 'youngModulus': young_modulus, + 'poissonRatio': poisson_ratio, + 'meshFile': _to_freefem_path(mesh_file), + 'outFile': _to_freefem_path(ff_out_path), + }) + + if not os.path.isfile(ff_out_path): + raise RuntimeError( + ) + + raw = np.loadtxt(ff_out_path) + _rewrite_freefem_output_in_place(raw, ff_out_path) + + x_ff, y_ff, z_ff = raw[:, 0], raw[:, 1], raw[:, 2] + ux_ff, uy_ff, uz_ff = raw[:, 3], raw[:, 4], raw[:, 5] + + # ========== Run SOFA =========== + pos0_sofa, u_sofa = sofaRun(mesh_file=mesh_file, T=T, radius=radius, + young_modulus=young_modulus, + poisson_ratio=poisson_ratio) + x_sofa, y_sofa, z_sofa = pos0_sofa[:, 0], pos0_sofa[:, 1], pos0_sofa[:, 2] + ux_sofa, uy_sofa, uz_sofa = u_sofa[:, 0], u_sofa[:, 1], u_sofa[:, 2] + + perm = _pair_by_coordinates(x_sofa, y_sofa, z_sofa, x_ff, y_ff, z_ff) + ux_ff_p = ux_ff[perm] + uy_ff_p = uy_ff[perm] + uz_ff_p = uz_ff[perm] + u_ff_p = np.column_stack([ux_ff_p, uy_ff_p, uz_ff_p]) + u_sofa_full = np.column_stack([ux_sofa, uy_sofa, uz_sofa]) + + yc = 0.5 * (y_sofa.min() + y_sofa.max()) + zc = 0.5 * (z_sofa.min() + z_sofa.max()) + u_ana, theta_prime = _analytical_displacement( + pos0_sofa, T, radius, young_modulus, poisson_ratio, yc, zc + ) + + rms_ux = _rms(ux_sofa, ux_ff_p) + rms_uy = _rms(uy_sofa, uy_ff_p) + rms_uz = _rms(uz_sofa, uz_ff_p) + + rel_sofa_ff = _rel_rms(u_sofa_full, u_ff_p) + rel_sofa_ana = _rel_rms(u_ana, u_sofa_full) + rel_ff_ana = _rel_rms(u_ana, u_ff_p) + + with open("results/comparison_beam3d_torsion_results.txt", 'w') as f: + header = (f"{'x':>10} {'y':>10} {'z':>10} {'ux_sofa':>14} {'ux_ff':>14} " + f"{'uy_sofa':>14} {'uy_ff':>14} {'uz_sofa':>14} {'uz_ff':>14}") + f.write(header + "\n") + f.write("-" * len(header) + "\n") + for x, y, z, uxs, uxf, uys, uyf, uzs, uzf in zip( + x_sofa, y_sofa, z_sofa, ux_sofa, ux_ff_p, uy_sofa, uy_ff_p, uz_sofa, uz_ff_p): + f.write(f"{x:10.4f} {y:10.4f} {z:10.4f} {uxs:+14.6e} {uxf:+14.6e} " + f"{uys:+14.6e} {uyf:+14.6e} {uzs:+14.6e} {uzf:+14.6e}\n") + + f.write("\n") + f.write(f"theta' (analytical) = {theta_prime:.6g} rad/m\n") + f.write("RMS norms\n") + f.write("-" * 40 + "\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") + f.write(f" RMS_uz (SOFA vs FF) = {rms_uz:.6e}\n") + f.write(f" Relatif SOFA vs FF = {rel_sofa_ff:.3%}\n") + f.write(f" Relatif SOFA vs Analytique = {rel_sofa_ana:.3%}\n") + f.write(f" Relatif FF vs Analytique = {rel_ff_ana:.3%}\n") + + print("=" * 70) + print("Validation : poutre 3D section circulaire, torsion pure") + print("=" * 70) + print(f"theta' (analytical) = {theta_prime:.6g} rad/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}") + print("-" * 70) + print(f"Relatif SOFA vs FF = {rel_sofa_ff:.3%}") + print(f"Relatif SOFA vs Analytique = {rel_sofa_ana:.3%}") + print(f"Relatif FF vs Analytique = {rel_ff_ana:.3%}") + + fig, axes = plt.subplots(1, 3, figsize=(15, 5)) + for ax, (u_s, u_f, label) in zip(axes, [ + (ux_sofa, ux_ff_p, 'ux'), (uy_sofa, uy_ff_p, 'uy'), (uz_sofa, uz_ff_p, 'uz')]): + ax.scatter(u_s, u_f, s=15, alpha=0.8) + lo = min(u_s.min(), u_f.min()) + hi = max(u_s.max(), u_f.max()) + ax.plot([lo, hi], [lo, hi], 'r--', linewidth=1) + ax.set_xlabel(f'{label}_sofa') + ax.set_ylabel(f'{label}_ff') + ax.set_title(label) + + fig.suptitle("3D Beam - Torsion- SOFA vs FreeFEM ", fontsize=14) + plt.tight_layout() + fig.savefig("results/comparison_beam3d_torsion_fields.png", dpi=150) + plt.close(fig) \ No newline at end of file diff --git a/examples/Freefem/validation/3D_torsion/freefem_beam3d_torsion.edp b/examples/Freefem/validation/3D_torsion/freefem_beam3d_torsion.edp new file mode 100644 index 00000000..89feb017 --- /dev/null +++ b/examples/Freefem/validation/3D_torsion/freefem_beam3d_torsion.edp @@ -0,0 +1,99 @@ +load "gmsh" +load "msh3" + +DEFAULT (T, 0.08) +DEFAULT (radius, 0.1) +DEFAULT (youngModulus, 10000.0) +DEFAULT (poissonRatio, 0.3) +DEFAULT (meshFile, "beam3d_circular_tet.msh") +DEFAULT (outFile, "freefem_beam3d_torsion_out.txt") + +real Torque = $T; +real radius = $radius; +real E = $youngModulus; +real nu = $poissonRatio; + +real mu = E / (2.*(1.+nu)); +real lambda = E*nu / ((1.+nu)*(1.-2.*nu)); + +mesh3 Th = gmshload3("$meshFile"); + +// =============== Geometrie =============== +real xmin=1e30, xmax=-1e30, ymin=1e30, ymax=-1e30, zmin=1e30, zmax=-1e30; +for (int i = 0; i < Th.nv; i++) { + xmin = min(xmin, Th(i).x); xmax = max(xmax, Th(i).x); + ymin = min(ymin, Th(i).y); ymax = max(ymax, Th(i).y); + zmin = min(zmin, Th(i).z); zmax = max(zmax, Th(i).z); +} +real length = xmax - xmin; +real yc = 0.5*(ymin+ymax); +real zc = 0.5*(zmin+zmax); + +real J = pi*radius^4/2.; +real G = mu; +real thetaPrime = Torque/(G*J); +real thetaTotal = thetaPrime*length; +cout << "length=" << length << " yc=" << yc << " zc=" << zc << " J=" << J << endl; +cout << "theta' = " << thetaPrime << " rad/m, theta_total = " << thetaTotal << " rad" << endl; +if (abs(thetaTotal) > 0.1) + cout << " total torsion's angle : small-strain " << endl; + + +fespace Wh(Th, [P1, P1, P1]); +Wh [ux, uy, uz], [vx, vy, vz]; + + +varf vElasticity([ux, uy, uz], [vx, vy, vz]) = + int3d(Th)( + lambda*(dx(ux)+dy(uy)+dz(uz))*(dx(vx)+dy(vy)+dz(vz)) + + 2.*mu*( dx(ux)*dx(vx) + dy(uy)*dy(vy) + dz(uz)*dz(vz) ) + + mu*( dy(ux)+dx(uy) )*( dy(vx)+dx(vy) ) + + mu*( dz(ux)+dx(uz) )*( dz(vx)+dx(vz) ) + + mu*( dz(uy)+dy(uz) )*( dz(vy)+dy(vz) ) + ); + +matrix A = vElasticity(Wh, Wh); +real[int] b(Wh.ndof); +b = 0.0; + +func taux = 0.; +func tauy = -(Torque/J)*(z-zc); +func tauz = (Torque/J)*(y-yc); + +varf vTraction([ux, uy, uz], [vx, vy, vz]) = int2d(Th, 2)(taux*vx + tauy*vy + tauz*vz); +real[int] bTraction = vTraction(0, Wh); +b += bTraction; + +cout << " Aire face charged (attendu " << pi*radius^2 << ") = " + << int2d(Th, 2)(1.) << endl; +cout << "Applicated moment (attendu " << Torque << ") = " + << int2d(Th, 2)((y-yc)*tauz - (z-zc)*tauy) << endl; + +real tgv = 1e30; +int nFixed = 0; +for (int i = 0; i < Th.nv; i++) { + if (abs(Th(i).x - xmin) < 1e-6) { + A(3*i, 3*i) = tgv; b[3*i] = 0.0; + A(3*i+1, 3*i+1) = tgv; b[3*i+1] = 0.0; + A(3*i+2, 3*i+2) = tgv; b[3*i+2] = 0.0; + nFixed++; + } +} + +set(A, solver=sparsesolver); +real[int] sol = A^-1 * b; +ux[] = sol; + +cout << " amplitude uy sur maillage : min=" << uy[].min + << " max=" << uy[].max << endl; + + +{ + ofstream fout("$outFile"); + fout.precision(12); + for (int i = 0; i < Th.nv; i++) { + fout << Th(i).x << " " << Th(i).y << " " << Th(i).z << " " + << sol[3*i] << " " << sol[3*i+1] << " " << sol[3*i+2] << endl; + } +} + diff --git a/examples/Freefem/validation/3D_torsion/generate_beam3d_circular_tet.py b/examples/Freefem/validation/3D_torsion/generate_beam3d_circular_tet.py new file mode 100644 index 00000000..2396321e --- /dev/null +++ b/examples/Freefem/validation/3D_torsion/generate_beam3d_circular_tet.py @@ -0,0 +1,72 @@ +import json +import os +import sys +import gmsh + +MESH_DIR = os.path.join(os.path.dirname(os.path.abspath(__file__)), "mesh") +DEFAULT_MESH_FILENAME = "beam3d_circular_tet.msh" + + +def generate_beam3D_circular_tet(length, radius, mesh_size, filename=None): + + if filename is None: + filename = DEFAULT_MESH_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, \ + "Verify the faces " + 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(MESH_DIR, exist_ok=True) + msh_path = os.path.join(MESH_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", DEFAULT_MESH_FILENAME), + ) + print("Wrote:", msh_path) + print("Number of nodes:", n_nodes) \ No newline at end of file diff --git a/examples/Freefem/validation/3D_torsion/params.json b/examples/Freefem/validation/3D_torsion/params.json new file mode 100644 index 00000000..f4f14328 --- /dev/null +++ b/examples/Freefem/validation/3D_torsion/params.json @@ -0,0 +1,8 @@ +{ + "beam3d_circle_tet": { + "length": 1.0, + "radius": 0.1, + "mesh_size": 0.06, + "meshfile": "beam3d_circular_tet.msh" + } +} \ No newline at end of file diff --git a/examples/Freefem/validation/3D_torsion/params_beam3d_torsion.json b/examples/Freefem/validation/3D_torsion/params_beam3d_torsion.json new file mode 100644 index 00000000..f9947124 --- /dev/null +++ b/examples/Freefem/validation/3D_torsion/params_beam3d_torsion.json @@ -0,0 +1,6 @@ +{ + "T": 0.09, + "radius": 0.1, + "youngModulus": 10000.0, + "poissonRatio": 0.3 +} \ No newline at end of file diff --git a/examples/Freefem/validation/3D_torsion/sofa_beam3d_torsion.py b/examples/Freefem/validation/3D_torsion/sofa_beam3d_torsion.py new file mode 100644 index 00000000..3e6619d1 --- /dev/null +++ b/examples/Freefem/validation/3D_torsion/sofa_beam3d_torsion.py @@ -0,0 +1,247 @@ +import json +import os +import sys +import numpy as np +import Sofa +import Sofa.Core +import Sofa.Simulation + +RESULTS_DIR = "results" +MESH_DIR = os.path.join(os.path.dirname(os.path.abspath(__file__)), "mesh") +DEFAULT_MESH_FILENAME = "beam3d_circular_tet.msh" + + +def torsion_consistent_forces(nodes, end_faces, T, J, yc, zc): + N = len(nodes) + F = np.zeros((N, 3)) + for tri in end_faces: + pts = nodes[tri, :] + v1 = pts[1] - pts[0] + v2 = pts[2] - pts[0] + area = 0.5 * np.linalg.norm(np.cross(v1, v2)) + for nid in tri: + y, z = nodes[nid, 1], nodes[nid, 2] + dy, dz = y - yc, z - zc + ty = -(T / J) * dz + tz = (T / J) * dy + F[nid, 1] += ty * area / 3.0 + F[nid, 2] += tz * area / 3.0 + return F + + +def _check_small_strain(T, radius, young_modulus, poisson_ratio, length, + theta_length_limit=0.1): + J = np.pi * radius**4 / 2.0 + G = young_modulus / (2.0 * (1.0 + poisson_ratio)) + theta = T / (G * J) + theta_total = theta * length + if abs(theta_total) > theta_length_limit: + print( + f"estimated total Torsion's Angle = {theta_total:.3g} rad " + f"(> {theta_length_limit} rad). small-strain linear modal it's not validated", + file=sys.stderr, + ) + return theta, theta_total + + +def _verify_torsion(x0, u, radius, yc, zc, theta_total, tol_x=1e-6, + radius_rel_tol=0.05, angle_abs_tol=0.05): + x = x0[:, 0] + x_max = x.max() + end_mask = np.isclose(x, x_max, atol=tol_x) + + y0e, z0e = x0[end_mask, 1], x0[end_mask, 2] + uxe, uye, uze = u[end_mask, 0], u[end_mask, 1], u[end_mask, 2] + + y1e, z1e = y0e + uye, z0e + uze + + r0 = np.sqrt((y0e - yc) ** 2 + (z0e - zc) ** 2) + r1 = np.sqrt((y1e - yc) ** 2 + (z1e - zc) ** 2) + valid = r0 > 0.1 * radius + + r_rel_err = np.abs(r1[valid] - r0[valid]) / r0[valid] + + angle0 = np.arctan2(z0e[valid] - zc, y0e[valid] - yc) + angle1 = np.arctan2(z1e[valid] - zc, y1e[valid] - yc) + dangle = np.mod(angle1 - angle0 + np.pi, 2 * np.pi) - np.pi + + + ok_radius = r_rel_err.max() < radius_rel_tol + ok_angle = abs(dangle.mean() - theta_total) < angle_abs_tol + + if ok_radius and ok_angle: + print(" It's a real torsion ") + else: + print( + " The deformation it's not a torsion ===> verify the T & youngModulus ", + file=sys.stderr, + ) + + return { + "r0": r0, "r1": r1, "r_rel_err": r_rel_err, + "dangle_mean": dangle.mean(), "dangle_std": dangle.std(), + "ux_mean": uxe.mean(), "ux_max_abs": np.abs(uxe).max(), + "ok": ok_radius and ok_angle, + } + +def _default_mesh_path(): + return os.path.join(MESH_DIR, DEFAULT_MESH_FILENAME) + + +def _default_params_path(): + return os.path.join(os.path.dirname(os.path.abspath(__file__)), "params_beam3d_torsion.json") + + +def create_scene_args(rootNode, mesh_file, T, radius, young_modulus, poisson_ratio, tol=1e-6): + if not os.path.isfile(mesh_file): + raise FileNotFoundError() + + requiredPlugins = [ + "Elasticity", + "Sofa.Component.Constraint.Projective", + "Sofa.Component.IO.Mesh", + "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"]) + + template = "Vec3d" + + 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") + + loader = Beam.addObject('MeshGmshLoader', name="loader", filename=mesh_file) + + nodes = np.array(loader.position.value) + tets = np.array(loader.tetrahedra.value) + tris = np.array(loader.triangles.value) + N = len(nodes) + + x_min = nodes[:, 0].min() + x_max = nodes[:, 0].max() + length = x_max - x_min + fixed_idx = np.where(np.isclose(nodes[:, 0], x_min, atol=tol))[0].tolist() + + end_mask = np.all(np.isclose(nodes[tris, 0], x_max, atol=tol), axis=1) + end_faces = tris[end_mask] + + if len(fixed_idx) == 0: + raise RuntimeError() + if len(end_faces) == 0: + raise RuntimeError() + + _check_small_strain(T, radius, young_modulus, poisson_ratio, length) + + yc = 0.5 * (nodes[:, 1].min() + nodes[:, 1].max()) + zc = 0.5 * (nodes[:, 2].min() + nodes[:, 2].max()) + J = np.pi * radius**4 / 2.0 + + F_nodal = torsion_consistent_forces(nodes, end_faces, T, J, yc, zc) + forces_list = F_nodal.tolist() + + dofs = Beam.addObject('MechanicalObject' + , name="dofs" + , template=template + , position="@loader.position" + , showObject=True + , showObjectScale=0.01) + + Beam.addObject('TetrahedronSetTopologyContainer' + , name="topology" + , src="@loader") + Beam.addObject('TetrahedronSetTopologyModifier') + + 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="TorqueTraction" + , indices=list(range(N)) + , forces=forces_list + , showArrowSize=0.0 + , showColor=[1.0, 0.2, 0.0, 1.0]) + + + return rootNode, dofs, nodes.copy() + + +def createScene(rootNode): + with open(_default_params_path()) as f: + cfg = json.load(f) + create_scene_args(rootNode + , mesh_file=_default_mesh_path() + , T=float(cfg["T"]) + , radius=float(cfg["radius"]) + , young_modulus=float(cfg["youngModulus"]) + , poisson_ratio=float(cfg["poissonRatio"])) + return rootNode + + +def sofaRun(mesh_file, T, radius, young_modulus, poisson_ratio): + root = Sofa.Core.Node("root") + _, dofs, pos0 = create_scene_args(root + , mesh_file=mesh_file + , T=T + , radius=radius + , 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[:, :3] - pos0[:, :3] + x0 = pos0[:, :3] + + os.makedirs(RESULTS_DIR, exist_ok=True) + out_path = os.path.join(RESULTS_DIR, "sofa_beam3d_torsion_results.txt") + with open(out_path, 'w') as f: + f.write(f"{'x0':>12} {'y0':>12} {'z0':>12} {'ux':>12} {'uy':>12} {'uz':>12}\n") + f.write("-" * 78 + "\n") + for (xi, yi, zi), (uxi, uyi, uzi) in zip(x0, u): + f.write(f"{xi:12.6f} {yi:12.6f} {zi:12.6f} {uxi:12.6f} {uyi:12.6f} {uzi:12.6f}\n") + + yc = 0.5 * (x0[:, 1].min() + x0[:, 1].max()) + zc = 0.5 * (x0[:, 2].min() + x0[:, 2].max()) + length = x0[:, 0].max() - x0[:, 0].min() + _, theta_total = _check_small_strain(T, radius, young_modulus, poisson_ratio, length) + _verify_torsion(x0, u, radius, yc, zc, theta_total) + + 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() + , T=float(cfg["T"]) + , radius=float(cfg["radius"]) + , young_modulus=float(cfg["youngModulus"]) + , poisson_ratio=float(cfg["poissonRatio"])) \ No newline at end of file