diff --git a/examples/Freefem/validation/3D/3D_circulaire/beam3d_circular_tet.py b/examples/Freefem/validation/3D/3D_circulaire/beam3d_circular_tet.py new file mode 100644 index 00000000..0c0223d1 --- /dev/null +++ b/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) \ No newline at end of file diff --git a/examples/Freefem/validation/3D/3D_circulaire/comparaison_script3d_circle.py b/examples/Freefem/validation/3D/3D_circulaire/comparaison_script3d_circle.py new file mode 100644 index 00000000..6abd6871 --- /dev/null +++ b/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() \ No newline at end of file diff --git a/examples/Freefem/validation/3D/3D_circulaire/freefem_beam3d_circle_traction.edp b/examples/Freefem/validation/3D/3D_circulaire/freefem_beam3d_circle_traction.edp new file mode 100644 index 00000000..3e4e584c --- /dev/null +++ b/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++; +} + diff --git a/examples/Freefem/validation/3D/3D_circulaire/params.json b/examples/Freefem/validation/3D/3D_circulaire/params.json new file mode 100644 index 00000000..2e2cae32 --- /dev/null +++ b/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" +} +} \ No newline at end of file diff --git a/examples/Freefem/validation/3D/3D_circulaire/params_beam3d_circle_traction.json b/examples/Freefem/validation/3D/3D_circulaire/params_beam3d_circle_traction.json new file mode 100644 index 00000000..8ba89712 --- /dev/null +++ b/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 +} \ No newline at end of file diff --git a/examples/Freefem/validation/3D/3D_circulaire/sofa_beam3d_circle_traction.py b/examples/Freefem/validation/3D/3D_circulaire/sofa_beam3d_circle_traction.py new file mode 100644 index 00000000..b2936da9 --- /dev/null +++ b/examples/Freefem/validation/3D/3D_circulaire/sofa_beam3d_circle_traction.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, 0] += 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_traction.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="LoadedTraction" + , 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_traction_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