diff --git a/examples/Freefem/validation/3D/3D_circulaire_downwards/beam3d_circular_tet.py b/examples/Freefem/validation/3D/3D_circulaire_downwards/beam3d_circular_tet.py new file mode 100644 index 00000000..0c0223d1 --- /dev/null +++ b/examples/Freefem/validation/3D/3D_circulaire_downwards/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) \ No newline at end of file diff --git a/examples/Freefem/validation/3D/3D_circulaire_downwards/comparaison_script3d_circle_bending.py b/examples/Freefem/validation/3D/3D_circulaire_downwards/comparaison_script3d_circle_bending.py new file mode 100644 index 00000000..b3d04f03 --- /dev/null +++ b/examples/Freefem/validation/3D/3D_circulaire_downwards/comparaison_script3d_circle_bending.py @@ -0,0 +1,130 @@ +# ============================================================================= +# comparaison_script3d_circle_bending.py +# 3D circular beam, transverse (bending) load: SOFA vs FreeFEM +# +# Same structure as comparaison_script3d_circle.py (axial traction case): +# generate the mesh, run FreeFEM, run SOFA, match nodes by coordinates, +# compute RMS. No analytical reference here -- the Saint-Venant traction +# formulas do not apply to bending; an Euler-Bernoulli/Timoshenko reference +# would need to be added separately if/when needed. +# ============================================================================= + +import os +import json +import numpy as np +import matplotlib.pyplot as plt + +from pyfreefem import FreeFemRunner +import sofa_beam3d_circle_bending as sofa_case + +RESULTS_DIR = "results" + + +def _default_edp_path(): + return os.path.join(os.path.dirname(os.path.abspath(__file__)), + "freefem_beam3d_circle_bending.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 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_bending.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"] + + 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, "F": q * (np.pi * radius**2), "radius": radius, "length": length, + }) + 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}") + + os.makedirs(RESULTS_DIR, exist_ok=True) + out_path = os.path.join(RESULTS_DIR, "comparison_beam3d_circle_bending_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") + 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 — Bending — SOFA vs FreeFEM (parity)") + fig_path = os.path.join(RESULTS_DIR, "comparison_beam3d_circle_bending_fields.png") + fig.savefig(fig_path, dpi=120, bbox_inches="tight") + print("Wrote:", fig_path) + + +if __name__ == "__main__": + main() \ No newline at end of file diff --git a/examples/Freefem/validation/3D/3D_circulaire_downwards/freefem_beam3d_circle_bending.edp b/examples/Freefem/validation/3D/3D_circulaire_downwards/freefem_beam3d_circle_bending.edp new file mode 100644 index 00000000..7c9a6e12 --- /dev/null +++ b/examples/Freefem/validation/3D/3D_circulaire_downwards/freefem_beam3d_circle_bending.edp @@ -0,0 +1,56 @@ +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) + +real E = $E; +real nu = $nu; +real F = $F; +real radius = $radius; +real length = $length; + +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*vy ) + + 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); diff --git a/examples/Freefem/validation/3D/3D_circulaire_downwards/params.json b/examples/Freefem/validation/3D/3D_circulaire_downwards/params.json new file mode 100644 index 00000000..2e2cae32 --- /dev/null +++ b/examples/Freefem/validation/3D/3D_circulaire_downwards/params.json @@ -0,0 +1,8 @@ +{ + "beam3d_circle_tet": { + "length": 1.0, + "radius": 0.1, + "mesh_size": 0.05, + "meshfile": "beam3d_circle_tet.msh" +} +} \ No newline at end of file diff --git a/examples/Freefem/validation/3D/3D_circulaire_downwards/params_beam3d_circle_bending.json b/examples/Freefem/validation/3D/3D_circulaire_downwards/params_beam3d_circle_bending.json new file mode 100644 index 00000000..fba09028 --- /dev/null +++ b/examples/Freefem/validation/3D/3D_circulaire_downwards/params_beam3d_circle_bending.json @@ -0,0 +1,8 @@ +{ + "length": 1.0, + "radius": 0.1, + "youngModulus": 210000.0, + "poissonRatio": 0.3, + "q": 150.9154943, + "exclusionFactor": 2.0 +} \ No newline at end of file diff --git a/examples/Freefem/validation/3D/3D_circulaire_downwards/sofa_beam3d_circle_bending.py b/examples/Freefem/validation/3D/3D_circulaire_downwards/sofa_beam3d_circle_bending.py new file mode 100644 index 00000000..73d81836 --- /dev/null +++ b/examples/Freefem/validation/3D/3D_circulaire_downwards/sofa_beam3d_circle_bending.py @@ -0,0 +1,168 @@ +import json +import os +import sys +import numpy as np +import Sofa +import Sofa.Core +import Sofa.Simulation + +RESULTS_DIR = "results" + + +def consistent_traction_forces(nodes, loaded_faces, q): + + N = len(nodes) + F = np.zeros((N, 3)) + for tri in loaded_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: + F[nid, 1] -= q * area / 3.0 + return F + + +def _default_mesh_path(): + return os.path.join(os.path.dirname(os.path.abspath(__file__)), + RESULTS_DIR, "beam3d_circle_tet.msh") + + +def _default_params_path(): + return os.path.join(os.path.dirname(os.path.abspath(__file__)), "params_beam3d_circle_bending.json") + + +def create_scene_args(rootNode, mesh_file, q, young_modulus, poisson_ratio, tol=1e-6): + 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") + + if not os.path.isfile(mesh_file): + raise FileNotFoundError( + f"Maillage introuvable : {mesh_file}\n" + f"-> generate it first : python beam3d_circular_tet.py params.json" + ) + + 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() + fixed_idx = np.where(np.isclose(nodes[:, 0], x_min, atol=tol))[0].tolist() + loaded_mask = np.all(np.isclose(nodes[tris, 0], x_max, atol=tol), axis=1) + loaded_faces = tris[loaded_mask] + + assert len(fixed_idx) > 0, "error " + assert len(loaded_faces) > 0, " (x_max) not found " + + F_nodal = consistent_traction_forces(nodes, loaded_faces, q) + 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="LoadedBending" + , indices=list(range(N)) + , forces=forces_list) + + 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() + , 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[:, :3] - pos0[:, :3] + x0 = pos0[:, :3] + + os.makedirs(RESULTS_DIR, exist_ok=True) + out_path = os.path.join(RESULTS_DIR, "sofa_beam3d_circle_bending_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") + + 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 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..f1dc8e0c --- /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.01, + "meshfile": "beam3d_circular_tet.msh" + } +} 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..14c4759a --- /dev/null +++ b/examples/Freefem/validation/3D_torsion/sofa_beam3d_torsion.py @@ -0,0 +1,274 @@ +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): + """ + Assemble les forces nodales consistantes pour la traction de torsion + t(y,z) = (-(T/J)*(z-zc), (T/J)*(y-yc)) sur la face chargee. + + Comme t est lineaire sur chaque triangle P1 et phi_i est lineaire, + F_i = int_T phi_i * t dA se calcule EXACTEMENT (pas de quadrature + approchee necessaire) via la formule de la matrice de masse P1 : + + F_i += (area / 12.0) * (2 * t(x_i) + t(x_j) + t(x_k)) + + Cela remplace l'ancienne version "lumpee" (F_i += t(x_i) * area / 3), + qui n'est exacte que si t est constante sur le triangle (ce qui etait + le cas pour la charge distribuee, mais pas pour la torsion). + """ + N = len(nodes) + F = np.zeros((N, 3)) + + def traction(nid): + y, z = nodes[nid, 1], nodes[nid, 2] + dy, dz = y - yc, z - zc + ty = -(T / J) * dz + tz = (T / J) * dy + return ty, tz + + 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)) + + t_local = [traction(nid) for nid in tri] # t_local[l] <-> tri[l] + + for local_i in range(3): + nid_i = tri[local_i] + local_j, local_k = [l for l in range(3) if l != local_i] + ty_i, tz_i = t_local[local_i] + ty_j, tz_j = t_local[local_j] + ty_k, tz_k = t_local[local_k] + + F[nid_i, 1] += (area / 12.0) * (2 * ty_i + ty_j + ty_k) + F[nid_i, 2] += (area / 12.0) * (2 * tz_i + tz_j + tz_k) + + 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