diff --git a/examples/Freefem/2D-beam/compare_sofa_freefem.py b/examples/Freefem/2D-beam/compare_sofa_freefem.py index c7ed3ed..569cd0b 100644 --- a/examples/Freefem/2D-beam/compare_sofa_freefem.py +++ b/examples/Freefem/2D-beam/compare_sofa_freefem.py @@ -54,106 +54,111 @@ def _write_freefem_results(path, x, y, u_x, u_y): 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)]: + for formulation2d in ["planeStrain", "planeStress"]: + print(f"Running comparison simulations for {formulation2d}") + # --- 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('\\', '/'), + 'formulation' : formulation2d + }) + 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", f"freefem_{formulation2d}_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('\\', '/'), + formulation=formulation2d + ) + 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(f"results/comparison_{formulation2d}_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(f"results/plot_{formulation2d}_displacement_field.png", dpi=150) + plt.close(fig) + + # Fig 2 — deformed beams overlaid 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.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") - fig.savefig(f"results/plot_error_{comp}.png", dpi=150) + ax.legend() + fig.savefig(f"results/plot_{formulation2d}_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_{formulation2d}_{comp}.png", dpi=150) + plt.close(fig) + diff --git a/examples/Freefem/2D-beam/freefem_beam2d.edp b/examples/Freefem/2D-beam/freefem_beam2d.edp index 00d5076..3b89990 100644 --- a/examples/Freefem/2D-beam/freefem_beam2d.edp +++ b/examples/Freefem/2D-beam/freefem_beam2d.edp @@ -6,9 +6,10 @@ DEFAULT (rhoMat, "7.85e-6") DEFAULT (gravity, "9810") DEFAULT (poissonRatio, "0.3") DEFAULT (meshFile, "") +DEFAULT (formulation, "planeStrain") real E = $youngModulus; -real sigma = $poissonRatio; +real nu = $poissonRatio; real rhoMat = $rhoMat; real gravity = $gravity; real F = -rhoMat * gravity; @@ -22,8 +23,12 @@ 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)); +real mu = E / (2*(1 + nu)); +IFEQ (formulation, "planeStrain") +real lambda = E*nu / ((1 + nu)*(1 - 2*nu)); // 3D but z is constrained +ELSE +real lambda = E*nu / (1 - nu*nu); // Pure 2D +ENDIF solve Elasticity([uu, vv], [w, s], solver=sparsesolver) = int2d(th)( diff --git a/examples/Freefem/2D-beam/sofa_beam2d.py b/examples/Freefem/2D-beam/sofa_beam2d.py index 9ea8d1c..dbb8f8d 100644 --- a/examples/Freefem/2D-beam/sofa_beam2d.py +++ b/examples/Freefem/2D-beam/sofa_beam2d.py @@ -37,7 +37,16 @@ def onAnimateEndEvent(self, event): f.write(f"{x:12.6f} {y:12.6f} {ux:12.6f} {uy:12.6f}\n") -def create_scene_args(rootNode, mesh_file, height, young_modulus, poisson_ratio, rho, gravity): +def create_scene_args(rootNode + , mesh_file + , height + , young_modulus + , poisson_ratio + , rho + , gravity + , formulation): + dim_template = "Vec3d" if formulation == "planeStrain" else "Vec2d" + requiredPlugins = [ "Elasticity", "Sofa.Component.Constraint.Projective", @@ -85,7 +94,7 @@ def create_scene_args(rootNode, mesh_file, height, young_modulus, poisson_ratio, , filename=mesh_file) dofs = Beam.addObject('MechanicalObject' , name="dofs" - , template="Vec3d" + , template=dim_template , src="@loader") with Beam.addChild('triangles') as Triangles: @@ -94,15 +103,16 @@ def create_scene_args(rootNode, mesh_file, height, young_modulus, poisson_ratio, , src="@../loader") Triangles.addObject('TriangleSetTopologyModifier') Triangles.addObject('TriangleSetGeometryAlgorithms' - , template="Vec3d" + , template=dim_template , drawTriangles=True) Triangles.addObject('MeshMatrixMass' , name="mass" + , template=f"{dim_template},{dim_template}" , massDensity=rho , topology="@topology") Triangles.addObject('LinearSmallStrainFEMForceField' , name="FEM" - , template="Vec3d" + , template=dim_template , youngModulus=young_modulus , poissonRatio=poisson_ratio , topology="@topology") @@ -111,18 +121,18 @@ def create_scene_args(rootNode, mesh_file, height, young_modulus, poisson_ratio, eps = 1e-5 Beam.addObject('BoxROI' , name="fixed_roi" - , template="Vec3d" + , template=dim_template , box=[-eps, -eps, -1.0, eps, height + eps, 1.0] , drawBoxes=True) Beam.addObject('FixedProjectiveConstraint' - , template="Vec3d" + , template=dim_template , indices="@fixed_roi.indices") os.makedirs(RESULTS_DIR, exist_ok=True) exporter = rootNode.addObject( DisplacementExporter( dofs_node = dofs, - output_file = os.path.join(RESULTS_DIR, "sofa_results.txt"), + output_file = os.path.join(RESULTS_DIR, f"sofa_{formulation}_results.txt"), name = "exportCtrl" ) ) @@ -140,19 +150,21 @@ def createScene(rootNode): , rho = float(cfg["rho"]) , gravity = float(cfg["gravity"]) , mesh_file = os.path.join(RESULTS_DIR, cfg.get("meshfile", "beam2d.msh")) + , formulation = "planeStrain" ) return rootNode -def sofaRun(height, young_modulus, poisson_ratio, rho, gravity, mesh_file=None): +def sofaRun(height, young_modulus, poisson_ratio, rho, gravity, mesh_file, formulation): root = Sofa.Core.Node("root") _, exporter = create_scene_args(root , mesh_file = mesh_file - , height = height + , height = height , young_modulus= young_modulus , poisson_ratio= poisson_ratio , rho = rho - , gravity = gravity) + , gravity = gravity + , formulation = formulation) Sofa.Simulation.init(root) Sofa.Simulation.animate(root, root.dt.value) return exporter.pos_initial, exporter.u @@ -169,4 +181,5 @@ def sofaRun(height, young_modulus, poisson_ratio, rho, gravity, mesh_file=None): , rho = float(cfg["rho"]) , gravity = float(cfg["gravity"]) , mesh_file = os.path.join(RESULTS_DIR, cfg.get("meshfile", "beam2d.msh")) + , formulation = "planeStrain" ) diff --git a/examples/Freefem/3D-Elasticity/3d-beam.edp b/examples/Freefem/3D-Elasticity/3d-beam.edp new file mode 100644 index 0000000..4caee62 --- /dev/null +++ b/examples/Freefem/3D-Elasticity/3d-beam.edp @@ -0,0 +1,66 @@ +IMPORT "io.edp" +load "msh3" +load "gmsh" +DEFAULT (youngModulus, 10000) +DEFAULT (poissonRatio, 0.3) +DEFAULT (rhoMat, 20) +DEFAULT (grav, 10) +DEFAULT (meshFile, "results/beam3d.msh") + +real E = $youngModulus; +real sigma = $poissonRatio; +real rho = $rhoMat; +real g = $grav; + +real mu = E / (2.*(1. + sigma)); +real lambda = E * sigma / ((1. + sigma) * (1. - 2.*sigma)); + +mesh3 Th = gmshload3("$meshFile"); + +fespace Vh(Th, [P1, P1, P1]); +Vh [ux, uy, uz], [vx, vy, vz]; + +macro Epsilon(u1,u2,u3) [dx(u1),dy(u2),dz(u3), + (dy(u3)+dz(u2))/sqrt(2.), + (dz(u1)+dx(u3))/sqrt(2.), + (dx(u2)+dy(u1))/sqrt(2.)] // EOM +macro Div(u1,u2,u3) (dx(u1)+dy(u2)+dz(u3)) // EOM + +problem Elasticite3D([ux,uy,uz],[vx,vy,vz], solver=sparsesolver) = + int3d(Th)( + lambda * Div(vx,vy,vz) * Div(ux,uy,uz) + + 2.*mu * (Epsilon(vx,vy,vz)' * Epsilon(ux,uy,uz)) + ) + + int3d(Th)(rho * g * vz) + + on(1, ux=0, uy=0, uz=0); // FIXED: tag 1 = face gauche x=0 (clamped) + +Elasticite3D; + +cout << "uz simule = " << uz[].linfty << " m" << endl; + +fespace Sh(Th, P1); +Sh ux2, uy2, uz2; +ux2 = ux; uy2 = uy; uz2 = uz; + +real[int] dispX(Th.nv); +real[int] dispY(Th.nv); +real[int] dispZ(Th.nv); +real[int] coordX(Th.nv); +real[int] coordY(Th.nv); +real[int] coordZ(Th.nv); + +for(int i = 0; i < Th.nv; i++){ + dispX[i] = ux2[][i]; + dispY[i] = uy2[][i]; + dispZ[i] = uz2[][i]; + coordX[i] = Th(i).x; + coordY[i] = Th(i).y; + coordZ[i] = Th(i).z; +} + +exportArray(dispX); +exportArray(dispY); +exportArray(dispZ); +exportArray(coordX); +exportArray(coordY); +exportArray(coordZ); diff --git a/examples/Freefem/3D-Elasticity/compare_ff_sofa.py b/examples/Freefem/3D-Elasticity/compare_ff_sofa.py new file mode 100644 index 0000000..2d85d4c --- /dev/null +++ b/examples/Freefem/3D-Elasticity/compare_ff_sofa.py @@ -0,0 +1,162 @@ +""" +3D Beam Simulation - Comparison FreeFEM vs SOFA +""" +import json +import os +import sys +import numpy as np +import matplotlib.pyplot as plt +from gmsh_generate_beam3D import generate_beam3D +from sofa_beam3d 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, z, u_x, u_y, u_z): + with open(path, 'w') as f: + f.write(f"{'x':>12} {'y':>12} {'z':>12} {'ux':>12} {'uy':>12} {'uz':>12}\n") + f.write("-" * 80 + "\n") + for xi, yi, zi, uxi, uyi, uzi in zip(x, y, z, u_x, u_y, u_z): + f.write(f"{xi:12.6f} {yi:12.6f} {zi:12.6f} {uxi:12.6f} {uyi:12.6f} {uzi: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"]) + width = float(cfg["width"]) + nx = int(cfg["nx"]) + ny = int(cfg["ny"]) + nz = int(cfg["nz"]) + rho = float(cfg["rho"]) + gravity = float(cfg["gravity"]) + young_modulus = float(cfg["youngModulus"]) + poisson_ratio = float(cfg["poissonRatio"]) + mesh_filename = cfg.get("meshfile", "beam3d.msh") + + # --- Run Gmsh --- + path_before_gmsh = os.environ.get('PATH', '') + msh_path, x_positions, y_positions, z_positions = generate_beam3D( + length=length, + height=height, + width=width, + nx=nx, + ny=ny, + nz=nz, + filename=mesh_filename + ) + + # --- Run FreeFEM --- + os.environ['PATH'] = path_before_gmsh + runner = FreeFemRunner("3d-beam.edp") + results = runner.execute({ + 'youngModulus' : young_modulus, + 'poissonRatio' : poisson_ratio, + 'rhoMat' : rho, + 'grav' : gravity, + 'meshFile' : os.path.abspath(msh_path).replace('\\', '/'), + }) + x_ff = results['coordX'] + y_ff = results['coordY'] + z_ff = results['coordZ'] + u_x_ff = results['dispX'] + u_y_ff = results['dispY'] + u_z_ff = results['dispZ'] + + os.makedirs("results", exist_ok=True) + _write_freefem_results( + os.path.join("results", "freefem_3d_results.txt"), + x_ff, y_ff, z_ff, u_x_ff, u_y_ff, u_z_ff + ) + + # --- Run SOFA --- + pos_initial, u_sofa = sofaRun( + height = height, + width = width, + 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] + z_sofa = pos_initial[:, 2] + u_x_sofa = u_sofa[:, 0] + u_y_sofa = u_sofa[:, 1] + u_z_sofa = u_sofa[:, 2] + + ## --- Compare Results --- + err_ux = u_x_sofa - u_x_ff + err_uy = u_y_sofa - u_y_ff + err_uz = u_z_sofa - u_z_ff + + with open("results/comparison_3d_results.txt", 'w') as f: + header = (f"{'x':>12} {'y':>12} {'z':>12} " + f"{'ux_ff':>12} {'ux_sofa':>12} " + f"{'uy_ff':>12} {'uy_sofa':>12} " + f"{'uz_ff':>12} {'uz_sofa':>12} " + f"{'err_ux':>12} {'err_uy':>12} {'err_uz':>12}") + f.write(header + "\n") + f.write("-" * len(header) + "\n") + for vals in zip(x_ff, y_ff, z_ff, + u_x_ff, u_x_sofa, + u_y_ff, u_y_sofa, + u_z_ff, u_z_sofa, + err_ux, err_uy, err_uz): + 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") + f.write(f" ||uz_sofa - uz_ff||_2 = {_l2(u_z_sofa, u_z_ff):.6e}\n") + + ## --- Plots --- + + + tol = min(height, width) * 0.1 + mask_ff = (np.abs(y_ff - height/2) < tol) & (np.abs(z_ff - width/2) < tol) + mask_sofa = (np.abs(y_sofa - height/2) < tol) & (np.abs(z_sofa - width/2) < tol) + + fig, axes = plt.subplots(3, 1, figsize=(12, 10), constrained_layout=True) + fig.suptitle("Displacement along beam axis — FreeFEM vs SOFA") + + for ax, comp, uff, usofa, label in zip( + axes, + ['ux', 'uy', 'uz'], + [u_x_ff, u_y_ff, u_z_ff], + [u_x_sofa, u_y_sofa, u_z_sofa], + ['$u_x$', '$u_y$', '$u_z$'], + ): + ax.scatter(x_ff[mask_ff], uff[mask_ff], s=10, label="FreeFEM", alpha=0.8) + ax.scatter(x_sofa[mask_sofa], usofa[mask_sofa], s=10, label="SOFA", alpha=0.8, marker='x') + ax.set_xlabel("x (m)") + ax.set_ylabel(f"{label} (m)") + ax.set_title(f"{label} along beam axis") + ax.legend() + + fig.savefig("results/plot_3d_displacement_axis.png", dpi=150) + plt.close(fig) + + # Fig 2 — deformed shape (projection XZ — plan de flexion) + fig, ax = plt.subplots(figsize=(14, 5)) + ax.scatter(x_ff + u_x_ff, z_ff + u_z_ff, s=4, label="FreeFEM", alpha=0.7) + ax.scatter(x_sofa + u_x_sofa, z_sofa + u_z_sofa, s=4, label="SOFA", alpha=0.7, marker='x') + ax.set_aspect('equal') + ax.set_title("Deformed beam (XZ projection) — FreeFEM vs SOFA") + ax.set_xlabel("x (m)") + ax.set_ylabel("z (m)") + ax.legend() + fig.savefig("results/plot_3d_deformed_beam_XZ.png", dpi=150) + plt.close(fig) + + \ No newline at end of file diff --git a/examples/Freefem/3D-Elasticity/gmsh_generate_beam3D.py b/examples/Freefem/3D-Elasticity/gmsh_generate_beam3D.py new file mode 100644 index 0000000..c9b5f66 --- /dev/null +++ b/examples/Freefem/3D-Elasticity/gmsh_generate_beam3D.py @@ -0,0 +1,107 @@ +import json +import os +import sys +import gmsh + +RESULTS_DIR = "results" + +def generate_beam3D(length, height, width, nx, ny, nz, filename): + """ + Generate a 3D beam mesh (tetrahedral elements) + """ + gmsh.initialize() + gmsh.model.add("beam3d") + + # Corner points (bottom rectangle) + 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 + 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() + + # Structured transfinite layout on base surface + 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) + + # Recombine triangles → quads on the base face before extrusion + #gmsh.model.mesh.setRecombine(2, surf) + + # Extrude surface to create volume + extruded = gmsh.model.geo.extrude( + [(2, surf)], 0, 0, width, + numElements=[nz], + #recombine=True, + ) + + gmsh.model.geo.synchronize() + + # The volume created by extrusion + volume_tag = extruded[1][1] + + # Ensure the volume respects nx/ny/nz — required for the subdivision to + # produce the correct number of tets. + # gmsh.model.mesh.setTransfiniteVolume(volume_tag) + + # Physical groups — dim=2 surface groups so FreeFEM's on() finds boundary triangles + # (in 3D, boundary elements are triangles on surfaces, not edges on curves) + _eps = 1e-9 + def _surfs(xmin, ymin, zmin, xmax, ymax, zmax): + return [t for _, t in gmsh.model.getEntitiesInBoundingBox( + xmin, ymin, zmin, xmax, ymax, zmax, dim=2)] + + gmsh.model.addPhysicalGroup(2, _surfs(-_eps, -_eps, -_eps, _eps, height+_eps, width+_eps), tag=1, name="Fixed") + gmsh.model.addPhysicalGroup(2, _surfs(-_eps, -_eps, -_eps, length+_eps, _eps, width+_eps), tag=2, name="Bottom") + gmsh.model.addPhysicalGroup(2, _surfs(length-_eps,-_eps, -_eps, length+_eps, height+_eps, width+_eps), tag=3, name="Right") + gmsh.model.addPhysicalGroup(2, _surfs(-_eps, height-_eps, -_eps, length+_eps, height+_eps, width+_eps), tag=4, name="Top") + gmsh.model.addPhysicalGroup(2, _surfs(-_eps, -_eps, -_eps, length+_eps, height+_eps, _eps), tag=6, name="BaseZ") + gmsh.model.addPhysicalGroup(2, _surfs(-_eps, -_eps, width-_eps, length+_eps, height+_eps, width+_eps), tag=7, name="TopZ") + gmsh.model.addPhysicalGroup(3, [volume_tag], tag=5, name="Beam") + #gmsh.option.setNumber("Mesh.SubdivisionAlgorithm", 2) + + gmsh.model.mesh.generate(3) + gmsh.model.mesh.setOrder(1) + + _, node_coords, _ = gmsh.model.mesh.getNodes() + x_positions = node_coords[0::3] + y_positions = node_coords[1::3] + z_positions = node_coords[2::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, z_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, z_positions = generate_beam3D( + length=float(cfg["length"]), + height=float(cfg["height"]), + width=float(cfg["width"]), + nx=int(cfg["nx"]), + ny=int(cfg["ny"]), + nz=int(cfg["nz"]), + filename=cfg.get("meshfile"), + ) + print("Wrote:", msh_path) diff --git a/examples/Freefem/3D-Elasticity/params.json b/examples/Freefem/3D-Elasticity/params.json new file mode 100644 index 0000000..b172549 --- /dev/null +++ b/examples/Freefem/3D-Elasticity/params.json @@ -0,0 +1,13 @@ +{ + "length": 1.0, + "height" : 0.2, + "width" : 0.2, + "nx": 20, + "ny": 7, + "nz": 7, + "rho": 20, + "gravity": 9.8, + "youngModulus": 210000, + "poissonRatio": 0.3, + "meshfile": "beam3d.msh" +} diff --git a/examples/Freefem/3D-Elasticity/sofa_beam3d.py b/examples/Freefem/3D-Elasticity/sofa_beam3d.py new file mode 100644 index 0000000..08248c2 --- /dev/null +++ b/examples/Freefem/3D-Elasticity/sofa_beam3d.py @@ -0,0 +1,191 @@ +""" +3D Beam Simulation - SOFA Scene +""" +import json +import os +import sys +import Sofa +import Sofa.Core +import Sofa.Simulation +import SofaRuntime + +RESULTS_DIR = "results" + + +class DisplacementExporter(Sofa.Core.Controller): + """Stores 3D nodal displacements and writes them to file.""" + + def __init__(self, dofs_node, output_file, *args, **kwargs): + super().__init__(*args, **kwargs) + self.dofs_node = dofs_node + self.output_file = output_file + self.pos_initial = None + self.u = None + + def onSimulationInitDoneEvent(self, event): + self.pos_initial = self.dofs_node.position.array()[:, :3].copy() + + def onAnimateEndEvent(self, event): + pos_final = self.dofs_node.position.array()[:, :3] + self.u = pos_final - self.pos_initial + + with open(self.output_file, 'w') as f: + f.write(f"{'x':>12} {'y':>12} {'z':>12} {'ux':>12} {'uy':>12} {'uz':>12} \n") + f.write("-" * 54 + "\n") + for (x, y,z), (ux, uy, uz) in zip(self.pos_initial, self.u): + f.write(f"{x:12.6f} {y:12.6f} {z:12.6f} {ux:12.6f} {uy:12.6f} {uz:12.6f} \n") + + +# FIXED: added `width` parameter (was used in BoxROI but never defined) +# FIXED: removed `formulation` parameter (no longer needed in 3D-only code) +def create_scene_args(rootNode + , mesh_file + , height + , width # ADDED: needed for BoxROI y-extent + , young_modulus + , poisson_ratio + , rho + , gravity): + dim_template = "Vec3d" + + requiredPlugins = [ + "Elasticity", + "Sofa.Component.Constraint.Projective", + "Sofa.Component.Engine.Select", + "Sofa.Component.IO.Mesh", + "Sofa.Component.LinearSolver.Direct", + "Sofa.Component.Mass", + "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.gravity = [0, 0, -gravity] + rootNode.dt = 1.0 + + rootNode.addObject('DefaultAnimationLoop') + rootNode.addObject('VisualStyle', displayFlags=["showBehaviorModels", "showForceFields", "showWireframe"]) + + with rootNode.addChild('Beam') as Beam: + Beam.addObject('NewtonRaphsonSolver' + , name="newtonSolver" + , printLog=False + , warnWhenLineSearchFails=True + , maxNbIterationsNewton=1 + , maxNbIterationsLineSearch=1 + , lineSearchCoefficient=1 + , relativeSuccessiveStoppingThreshold=0 + , absoluteResidualStoppingThreshold=1e-7 + , absoluteEstimateDifferenceThreshold=1e-12 + , relativeInitialStoppingThreshold=1e-12 + , relativeEstimateDifferenceThreshold=0) + Beam.addObject('SparseLDLSolver' + , name="linearSolver" + , template="CompressedRowSparseMatrixd") + Beam.addObject('StaticSolver' + , name="staticSolver" + , newtonSolver="@newtonSolver" + , linearSolver="@linearSolver") + + Beam.addObject('MeshGmshLoader' + , name="loader" + , filename=mesh_file) + dofs = Beam.addObject('MechanicalObject' + , name="dofs" + , template=dim_template + , src="@loader") + + with Beam.addChild('tetrahedra') as Tetrahedra: + Tetrahedra.addObject('TetrahedronSetTopologyContainer' + , name="topology" + , src="@../loader") + Tetrahedra.addObject('TetrahedronSetTopologyModifier') + Tetrahedra.addObject('TetrahedronSetGeometryAlgorithms' + , template=dim_template + , drawTriangles=True) + Tetrahedra.addObject('MeshMatrixMass' + , name="mass" + , template=f"{dim_template},{dim_template}" + , massDensity=rho + , topology="@topology") + Tetrahedra.addObject('LinearSmallStrainFEMForceField' + , name="FEM" + , template=dim_template + , youngModulus=young_modulus + , poissonRatio=poisson_ratio + , topology="@topology") + + # Fix the left wall (x=0, label 4 in gmsh matching FreeFEM's on(4,...)) + eps = 1e-5 + Beam.addObject('BoxROI' + , name="fixed_roi" + , template=dim_template + , box=[-eps, -eps, -eps, eps, width + eps, height + eps] # FIXED: width now defined + , drawBoxes=True) + Beam.addObject('FixedProjectiveConstraint' + , template=dim_template + , indices="@fixed_roi.indices") + + os.makedirs(RESULTS_DIR, exist_ok=True) + exporter = rootNode.addObject( + DisplacementExporter( + dofs_node = dofs, + output_file = os.path.join(RESULTS_DIR, "sofa_3d_results.txt"), # FIXED: removed formulation (no longer a parameter) + name = "exportCtrl" + ) + ) + + return rootNode, exporter + + +def createScene(rootNode): + with open("params.json") as f: + cfg = json.load(f) + create_scene_args(rootNode + , mesh_file = os.path.join(RESULTS_DIR, cfg.get("meshfile", "beam3d.msh")) + , height = float(cfg["height"]) + , width = float(cfg["width"]) # ADDED: read width from params.json + , young_modulus= float(cfg["youngModulus"]) + , poisson_ratio= float(cfg["poissonRatio"]) + , rho = float(cfg["rho"]) + , gravity = float(cfg["gravity"]) + # FIXED: removed formulation= (no longer a parameter) + ) + return rootNode + + +# FIXED: removed `formulation` parameter, added `width` +def sofaRun(height, width, young_modulus, poisson_ratio, rho, gravity, mesh_file): + root = Sofa.Core.Node("root") + _, exporter = create_scene_args(root + , mesh_file = mesh_file + , height = height + , width = width # ADDED + , young_modulus= young_modulus + , poisson_ratio= poisson_ratio + , rho = rho + , gravity = gravity + # FIXED: removed formulation= + ) + Sofa.Simulation.init(root) + Sofa.Simulation.animate(root, root.dt.value) + return exporter.pos_initial, exporter.u + + +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) + + sofaRun(height = float(cfg["height"]) + , width = float(cfg["width"]) # ADDED + , young_modulus = float(cfg["youngModulus"]) + , poisson_ratio = float(cfg["poissonRatio"]) + , rho = float(cfg["rho"]) + , gravity = float(cfg["gravity"]) + , mesh_file = os.path.join(RESULTS_DIR, cfg.get("meshfile", "beam3d.msh")) + # FIXED: removed formulation= + ) \ No newline at end of file