From ce50c3a9f8caa26dd2c84364dd063cfa44792f29 Mon Sep 17 00:00:00 2001 From: fatima Date: Fri, 13 Mar 2026 16:01:10 +0100 Subject: [PATCH 1/5] Add FreeFem++ elasticity simulation files --- examples/Freefem/3d-cube.edp | 78 +++++++++++++++ examples/Freefem/barre-1d-solution-man.edp | 44 +++++++++ examples/Freefem/barre_exacte.edp | 42 +++++++++ examples/Freefem/beam-2d.edp | 104 ++++++++++++++++++++ examples/Freefem/elasticity_3D.edp | 105 +++++++++++++++++++++ 5 files changed, 373 insertions(+) create mode 100644 examples/Freefem/3d-cube.edp create mode 100644 examples/Freefem/barre-1d-solution-man.edp create mode 100644 examples/Freefem/barre_exacte.edp create mode 100644 examples/Freefem/beam-2d.edp create mode 100644 examples/Freefem/elasticity_3D.edp diff --git a/examples/Freefem/3d-cube.edp b/examples/Freefem/3d-cube.edp new file mode 100644 index 00000000..74916a82 --- /dev/null +++ b/examples/Freefem/3d-cube.edp @@ -0,0 +1,78 @@ +load "msh3" +load "medit" +load "iovtk" + + +// Pb elasticite dans un cube +mesh3 Th = cube(10, 10, 10); + +// Vérifier les labels du cube +cout << "Labels bords : " << endl; +int[int] labs = labels(Th); +for(int i = 0; i < labs.n; i++) + cout << " label " << labs[i] << endl; + +real E = 210e9; +real sigma = 0.3; +real rho = 7800.; +real gravity = -9.81 * rho; +real mu = E / (2.*(1.+sigma)); +real lambda = E*sigma / ((1.+sigma)*(1.-2.*sigma)); + +// Zmin/zmax +real zmin = 1e100, zmax = -1e100; +for (int i = 0; i < Th.nv; i++) { + zmin = min(zmin, Th(i).z); + zmax = max(zmax, Th(i).z); +} +cout << "zmin=" << zmin << " zmax=" << zmax << endl; + +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)(gravity * vz) + + on(1, ux=0, uy=0, uz=0); // encastrement + +Elasticite3D; + +cout << "Deplacement max uz = " << uz[].linfty << " m" << endl; +cout << "Theorique uz_max = " << rho*9.81*(zmax-zmin)^2/(2.*E) << " m" << endl; + +fespace Wh(Th, P1); +Wh normu = sqrt(ux^2 + uy^2 + uz^2); +medit("Norme", Th, normu, wait=1); +medit("uz", Th, uz, wait=1); + +real umax = max(ux[].linfty, max(uy[].linfty, uz[].linfty)); +real coef = 0.1*(zmax-zmin)/umax; +cout << "Coefficient amplification = " << coef << endl; +mesh3 Thd = movemesh3(Th, transfo=[x+coef*ux, y+coef*uy, z+coef*uz]); +//savemesh(Thd, "cube_deforme.mesh"); +medit("Cube deforme", Thd, wait=1); + + + + + + + +// Exporter la solution en format VTK +int[int] order = [1, 1, 1]; // Ordre P1 +savevtk("cube_deplacement.vtu", Th, [ux, uy, uz], dataname="Deplacement", order=order); +savevtk("cube_norme.vtu", Th, normu, dataname="Norme", order=order); +savevtk("cube_uz.vtu", Th, uz, dataname="Uz", order=order); + +// Exporter le maillage déformé +mesh3 Thd = movemesh3(Th, transfo=[x + 1e6*ux, y + 1e6*uy, z + 1e6*uz]); +savevtk("cube_deforme.vtu", Thd, [ux, uy, uz], dataname="Deplacement", order=order); + +cout << "Fichiers VTK exportés. Ouvrez-les avec ParaView." << endl; \ No newline at end of file diff --git a/examples/Freefem/barre-1d-solution-man.edp b/examples/Freefem/barre-1d-solution-man.edp new file mode 100644 index 00000000..406c39f2 --- /dev/null +++ b/examples/Freefem/barre-1d-solution-man.edp @@ -0,0 +1,44 @@ +// Solution manufacturée normalisée +// On résout : -u'' = pi^2*sin(pi*x) EA simplifie =1 +real L = 1.0; +real E = 1.0; +real A = 1.0; +real pi = 3.14159265358979; + +int[int] nvals = [2, 4, 8, 16, 32, 64]; + +ofstream file("convergence_grossier.txt"); +file << "# n\th\terrL2\terrH1" << endl; + +for (int idx = 0; idx < nvals.n; idx++) { + int n = nvals[idx]; + real h = L / n; + + meshL Th = segment(n, [x * L]); + fespace Vh(Th, P1); + Vh u, v; + + problem Traction(u, v) = + int1d(Th)( dx(u) * dx(v) ) + - int1d(Th)( pi * pi * sin(pi * x) * v ) + + on(1, u = 0) + + on(2, u = 0); + + Traction; + + Vh uexact = sin(pi * x); + + // Erreur sur les noeuds uniquement +real errL2 = 0, errH1 = 0; + +errL2 = sqrt(int1d(Th)( (u - sin(pi*x))^2 )); +errH1 = sqrt(int1d(Th)( (dx(u) - pi*cos(pi*x))^2 )); + + + cout << "n = " << n + << ", h = " << h + << ", errL2 = " << errL2 + << ", errH1 = " << errH1 << endl; + + file << n << "\t" << h << "\t" << errL2 << "\t" << errH1 << endl; +} \ No newline at end of file diff --git a/examples/Freefem/barre_exacte.edp b/examples/Freefem/barre_exacte.edp new file mode 100644 index 00000000..55adbef6 --- /dev/null +++ b/examples/Freefem/barre_exacte.edp @@ -0,0 +1,42 @@ +// === Élasticité linéaire - BARRE en traction === +// Version pour étude de convergence + +real L = 1.0; +real E = 1e6; +real A = 0.001; +real F = 1000; + +// Différentes valeurs de n pour l'étude de convergence + +int[int] nvals = [2, 4, 8, 16, 32, 64]; +ofstream file("convergence_data.txt"); +file << "# n\th\terrL2\terrH1" << endl; + +for (int idx = 0; idx < nvals.n; idx++) { + int n = nvals[idx]; + real h = L/n; + + meshL Th = segment(n, [x*L]); + fespace Vh(Th, P1); + Vh u, v; + + problem Traction(u, v) = + int1d(Th)( E * A * dx(u) * dx(v) ) + - int0d(Th, 2)( F * v ) + + on(1, u=0); + + Traction; + + Vh uexact = (F/(E*A)) * x; + + // Erreur L2 + real errL2 = sqrt(int1d(Th)((u - uexact)^2)); + + // Erreur H1 (semi-norme) + real errH1 = sqrt(int1d(Th)((dx(u) - dx(uexact))^2)); + + cout << "n = " << n << ", h = " << h << ", errL2 = " << errL2 << ", errH1 = " << errH1 << endl; + + // Sauvegarde + file << n << "\t " << h << "\t " << errL2 << "\t " << errH1 << endl; +} \ No newline at end of file diff --git a/examples/Freefem/beam-2d.edp b/examples/Freefem/beam-2d.edp new file mode 100644 index 00000000..2f9269ef --- /dev/null +++ b/examples/Freefem/beam-2d.edp @@ -0,0 +1,104 @@ +// Beam under its own weight + +// Parameters +real E = 21.5; +real sigma = 0.29; +real gravity = -0.05; +int n=140; +int m=35; + + +// Maillage non structure " avec courbes parametriques " +border a(t=2, 0){x=0; y=t ;label=1;} // gauche cote de la poutre à encastre +border b(t=0, 10){x=t; y=0 ;label=2;} // bas (libre) +border c(t=0, 2){x=10; y=t ;label=3;} // droite (libre) +border d(t=0, 10){x=10-t; y=2; label=4;} // haut (libre) +mesh th = buildmesh(b(n) + c(m) + d(n) + a(m)); // un maillage equilibre previlige celui là +// n=100 unite de long pour les cotes haut et bas, et n=25 unite de haut pour les cote gauche et droite +//mesh th = buildmesh(b(n) + c(n) + d(n) + a(n)); // evite celui là le maillage est raffine uniquement aux extrimite a et c aka gauche et droite +// Fespace element finis espace P1 ; +plot (th,cmm="maillage-poutre", wait=1); +fespace Vh(th, [P1, P1]); +Vh [uu, vv], [w, s]; + +// Macro +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 + +// Coefficients de Lamé +real mu = E/(2*(1 + sigma)); +real lambda = E*sigma/((1 + sigma)*(1 - 2*sigma)); + +// Problème +solve Elasticity ([uu, vv], [w, s], solver=sparsesolver) + = int2d(th)( + lambda*div(w, s)*div(uu, vv) + + 2.*mu*(epsilon(w, s)' * epsilon(uu, vv)) + ) + - int2d(th)(gravity*s) + + on(1, uu=0, vv=0); + +// Résultats +cout << "Max displacement = " << uu[].linfty << ", " << vv[].linfty << endl; + +// Visualisation +plot([uu, vv], wait=1, cmm="Deplacement"); + + +mesh th1 = movemesh(th, [x + uu, y + vv]); +plot( th1, wait=1, cmm="Maillage initial et deformation"); + +// Espace pour les contraintes +fespace Wh(th, P1); +Wh sigmaxx, sigmayy, sigmaxy, sigmavm; + +// Calcul des contraintes +sigmaxx = lambda * (dx(uu) + dy(vv)) + 2*mu*dx(uu); +sigmayy = lambda * (dx(uu) + dy(vv)) + 2*mu*dy(vv); +sigmaxy = mu*(dy(uu) + dx(vv)); + +// Contrainte de von Mises pour plasticité +sigmavm = sqrt(sigmaxx^2 + sigmayy^2 - sigmaxx*sigmayy + 3*sigmaxy^2); + +// Visualisation +plot(sigmaxx, wait=1, fill=1, value=1, cmm="Contrainte sigma_xx (MPa)"); +plot(sigmavm, wait=1, fill=1, value=1, cmm="Contrainte von Mises (MPa)"); + +// Valeurs extrêmes +cout << "sigma_xx min = " << sigmaxx[].min << ", max = " << sigmaxx[].max << endl; +cout << "sigma_vm max = " << sigmavm[].max << endl; + + + +// Calcul des réactions par intégration sur le bord encastré (label=1) +real reactionx = int1d(th, 1)( (sigmaxx * N.x + sigmaxy * N.y) ); +real reactiony = int1d(th, 1)( (sigmaxy * N.x + sigmayy * N.y) ); + +// Poids total +real poidstotal = abs(gravity) * int2d(th)(1); + +// Affichage des résultats + +cout << "Reaction horizontale a l'encastrement : " << reactionx << endl; +cout << "Reaction verticale a l'encastrement : " << reactiony << endl; +cout << "Poids total de la poutre : " << poidstotal << endl; +cout << "Rapport reaction verticale / poids : " << reactiony / poidstotal << endl; + +// Test d'équilibre +real erreurequilibre = abs(reactiony/poidstotal - 1) * 100; +cout << "Erreur d'equilibre : " << erreurequilibre << " %" << endl; + +if (erreurequilibre < 5) + cout << "equilibre verifie " << endl; +else + cout << " Attention : desequilibre detecte !" << endl; + +// Vérification du moment à l'encastrement +real moment = int1d(th, 1)( (y - 1) * (sigmaxy * N.x + sigmayy * N.y) - (x) * (sigmaxx * N.x + sigmaxy * N.y) ); +cout << "Moment a l'encastrement : " << moment << endl; + + + + + diff --git a/examples/Freefem/elasticity_3D.edp b/examples/Freefem/elasticity_3D.edp new file mode 100644 index 00000000..902cc751 --- /dev/null +++ b/examples/Freefem/elasticity_3D.edp @@ -0,0 +1,105 @@ +load "msh3" +load "medit" +load "iovtk" + // remarque : la comparaison theorique qu'on a obtenue n'est pas vraiment valide + // car on compare une geometrie tres difficle à une poutre, mais on reste tjrs dans le cadre du respect de la loi d la physique +mesh3 Th("C:/Program Files (x86)/FreeFem++/examples/3d/dodecaedre01.mesh"); + +// Calculer zmin +real zmin = 1e100, zmax = -1e100; +for (int i = 0; i < Th.nv; i++) { + zmin = min(zmin, Th(i).z); + zmax = max(zmax, Th(i).z); +} +cout << "zmin=" << zmin << " zmax=" << zmax << endl; + +// Changer le label des faces du bas (proches de zmin) +Th = change(Th, flabel=(abs(z - zmin) < 0.1) ? 1 : 0); + +// Vérifier les labels modifiés +cout << "Nouveaux labels :" << endl; +for (int i = 0; i < Th.nbe; i++) { + if (Th.be(i).label == 1) + cout << "Face " << i << " label = " << Th.be(i).label << endl; +} + +// Vérifier tous les labels présents +int[int] labs = labels(Th); +cout << "Labels presents dans le maillage :" << endl; +for(int i = 0; i < labs.n; i++) + cout << " label " << labs[i] << endl; + +// Paramètres physiques +real E = 210e9; +real sigma = 0.3; +real rho = 7800.; +real gravity = -9.81 * rho; +real mu = E / (2.*(1.+sigma)); +real lambda = E*sigma / ((1.+sigma)*(1.-2.*sigma)); + +// Espace éléments finis +fespace Vh(Th, [P1,P1,P1]); +Vh [ux,uy,uz],[vx,vy,vz]; + +// Macros (attention aux parenthèses superflues) +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 + +// Problème d'élasticité on fixe les faces label=1 +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)(gravity * vz) + + on(1, ux=0, uy=0, uz=0); + +// Résolution +Elasticite3D; + +// Résultats +cout << "Deplacement max |ux| = " << ux[].linfty << " m" << endl; +cout << "Deplacement max |uy| = " << uy[].linfty << " m" << endl; +cout << "Deplacement max |uz| = " << uz[].linfty << " m" << endl; + +// Création d'un espace scalaire P1 +fespace Wh(Th, P1); +Wh normu; + + +normu = sqrt(ux*ux + uy*uy + uz*uz); + +cout << "Deplacement max total = " << normu[].linfty << " m" << endl; + + + +cout << "Theorique (poutre) = " << rho*9.81*(zmax-zmin)^2/(2.*E) << " m" << endl; + +// Visualisation +medit("Norme du deplacement", Th, normu, wait=1); +medit("Deplacement vertical uz", Th, uz, wait=1); + +// Maillage déformé +real umax = max(ux[].linfty, max(uy[].linfty, uz[].linfty)); +real coef = 0.1*(zmax-zmin)/umax; +mesh3 Thd = movemesh3(Th, transfo=[x+coef*ux, y+coef*uy, z+coef*uz]); +medit("Dodecaedre deforme", Thd, wait=1); + + + + + + +// Exporter la solution en format VTK +//int[int] order = [1, 1, 1]; +//savevtk("cube_deplacement.vtu", Th, [ux, uy, uz], dataname="Deplacement", order=order); +//savevtk("cube_norme.vtu", Th, normu, dataname="Norme", order=order); +//savevtk("cube_uz.vtu", Th, uz, dataname="Uz", order=order); + +// Exporter le maillage déformé +//mesh3 Thd = movemesh3(Th, transfo=[x + 1e6*ux, y + 1e6*uy, z + 1e6*uz]); +//savevtk("cube_deforme.vtu", Thd, [ux, uy, uz], dataname="Deplacement", order=order); + From c5a92f94362d75bf487fe8b23d5e0830eb4e3817 Mon Sep 17 00:00:00 2001 From: fatima Date: Fri, 27 Mar 2026 17:01:33 +0100 Subject: [PATCH 2/5] Ajout des scripts Python pour simulations barre 1D et beam 2D avec FreeFem et SOFA --- barre_1d_pyfreefem.py | 91 +++++ barre_1d_sofa.py | 120 +++++++ beam_2d_freefem.py | 114 ++++++ beam_2d_sofa.py | 213 +++++++++++ comparaison.py | 796 ++++++++++++++++++++++++++++++++++++++++++ 5 files changed, 1334 insertions(+) create mode 100644 barre_1d_pyfreefem.py create mode 100644 barre_1d_sofa.py create mode 100644 beam_2d_freefem.py create mode 100644 beam_2d_sofa.py create mode 100644 comparaison.py diff --git a/barre_1d_pyfreefem.py b/barre_1d_pyfreefem.py new file mode 100644 index 00000000..7c5d4d8d --- /dev/null +++ b/barre_1d_pyfreefem.py @@ -0,0 +1,91 @@ +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt +import subprocess +import tempfile +import os + +# PARAMÈTRES +L = 1.0 +E = 1e6 +A = 0.001 +N = 64 +OUTPUT_FF = "freefem_deplacement.txt" + + +FREEFEM_CODE = """ +real L = {L}; +real E = {E}; +real A = {A}; +int n = {N}; + +meshL Th = segment(n, [x*L]); +fespace Vh(Th, P1); +Vh u, v; + +problem Traction(u, v) = + int1d(Th)( E * A * dx(u) * dx(v) ) + - int1d(Th)( 1.0 * v ) + + on(1, u=0); + +Traction; + +ofstream file("{OUTPUT_FF}"); +file << "x u" << endl; +for (int i = 0; i < Th.nv; i++) {{ + real xi = Th(i).x; + real ui = u(xi, 0); + file << xi << " " << ui << endl; +}} + +cout << "Export OK : " << Th.nv << " noeuds ecrits dans {OUTPUT_FF}" << endl; +""".format(L=L, E=E, A=A, N=N, OUTPUT_FF=OUTPUT_FF) + +# EXÉCUTION FREEFEM +print("Lancement FreeFEM++...") +with tempfile.NamedTemporaryFile(mode='w', suffix='.edp', + delete=False, encoding='utf-8') as f: + f.write(FREEFEM_CODE) + temp_file = f.name + +try: + result = subprocess.run( + ['FreeFem++', temp_file], + capture_output=True, text=True, encoding='utf-8' + ) + print(result.stdout) + if result.returncode != 0: + print("ERREUR FreeFEM :\n", result.stderr) + exit(1) +finally: + os.unlink(temp_file) + +# LECTURE RÉSULTATS +if not os.path.exists(OUTPUT_FF): + print(f"Fichier {OUTPUT_FF} non trouvé — vérifier FreeFEM.") + exit(1) + +df_ff = pd.read_csv(OUTPUT_FF, sep=r'\s+') +df_ff = df_ff.sort_values('x').reset_index(drop=True) +x_ff = df_ff['x'].values +u_ff = df_ff['u'].values +#print(f"\nFreeFEM : {len(x_ff)} nœuds chargés depuis '{OUTPUT_FF}'") +#print(df_ff.to_string(index=False)) + +# SOLUTION EXACTE +# -EA u'' = 1, u(0)=0, u'(L)=0 → u(x) = x*(L - x/2) / (EA) +x_exact = np.linspace(0, L, 500) +u_exact = x_exact * (L - x_exact / 2.0) / (E * A) + +# ─── FIGURE +fig, ax = plt.subplots(figsize=(7, 5)) + +ax.plot(x_exact, u_exact, '-', color='gray', lw=1.5, label='Solution exacte') +ax.plot(x_ff, u_ff, 'o-', color='#1a6faf', lw=2, markersize=4, + label=f'FreeFEM (N={N})') +ax.set_xlabel("x (m)") +ax.set_ylabel("u(x) (m)") +ax.set_title("Déplacement axial — Barre 1D ") +ax.legend() +ax.grid(True, alpha=0.3) + diff --git a/barre_1d_sofa.py b/barre_1d_sofa.py new file mode 100644 index 00000000..7e8ebd47 --- /dev/null +++ b/barre_1d_sofa.py @@ -0,0 +1,120 @@ +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt +import subprocess +import tempfile +import os + +# ─── PARAMÈTRES +L = 1.0 +E = 1e6 +A = 0.001 +N = 64 # Nombre d'éléments +OUTPUT_FF = "freefem_deplacement.txt" # Sortie FreeFEM +OUTPUT_SOFA = "sofa_deplacement.txt" # Fichier SOFA à exporter + + +FREEFEM_CODE = """ +real L = {L}; +real E = {E}; +real A = {A}; +int n = {N}; + +meshL Th = segment(n, [x*L]); +fespace Vh(Th, P1); +Vh u, v; + +problem Traction(u, v) = + int1d(Th)( E * A * dx(u) * dx(v) ) + - int1d(Th)( 1.0 * v ) + + on(1, u=0); + +Traction; + +// Export noeuds : on parcourt les points du maillage +ofstream file("{OUTPUT_FF}"); +file << "x u" << endl; +for (int i = 0; i < Th.nv; i++) {{ + real xi = Th(i).x; + real ui = u(xi, 0); + file << xi << " " << ui << endl; +}} + +cout << "Export OK : " << Th.nv << " noeuds ecrits dans {OUTPUT_FF}" << endl; +""".format(L=L, E=E, A=A, N=N, OUTPUT_FF=OUTPUT_FF) + +# EXÉCUTION FREEFEM +print("Lancement FreeFEM++...") +with tempfile.NamedTemporaryFile(mode='w', suffix='.edp', + delete=False, encoding='utf-8') as f: + f.write(FREEFEM_CODE) + temp_file = f.name + +try: + result = subprocess.run( + ['FreeFem++', temp_file], + capture_output=True, text=True, encoding='utf-8' + ) + print(result.stdout) + if result.returncode != 0: + print("ERREUR FreeFEM :\n", result.stderr) + exit(1) +finally: + os.unlink(temp_file) + +# ─── LECTURE FREEFEM +if not os.path.exists(OUTPUT_FF): + print(f"Fichier {OUTPUT_FF} non trouvé — vérifier FreeFEM.") + exit(1) + +df_ff = pd.read_csv(OUTPUT_FF, sep=r'\s+') +x_ff = df_ff['x'].values +u_ff = df_ff['u'].values +print(f"\nFreeFEM : {len(x_ff)} nœuds chargés depuis '{OUTPUT_FF}'") + +# ─── SOLUTION EXACTE +# -EA u'' = 1, u(0)=0, u'(L)=0 → u(x) = (1/EA) * x*(L - x/2) +x_exact = np.linspace(0, L, 500) +u_exact = (1.0 / (E * A)) * x_exact * (L - x_exact / 2.0) + +# ─── LECTURE SOFA +sofa_disponible = os.path.exists(OUTPUT_SOFA) +if sofa_disponible: + df_sofa = pd.read_csv(OUTPUT_SOFA, sep=r'\s+') + # Colonnes attendues : x u (même convention que FreeFEM) + x_sofa = df_sofa.iloc[:, 0].values + u_sofa = df_sofa.iloc[:, 1].values + print(f"SOFA : {len(x_sofa)} nœuds chargés depuis '{OUTPUT_SOFA}'") + + # Interpolation FreeFEM sur les points SOFA pour comparaison ponctuelle + u_ff_interp = np.interp(x_sofa, x_ff, u_ff) + diff = u_sofa - u_ff_interp + err_max = np.max(np.abs(diff)) + err_rms = np.sqrt(np.mean(diff**2)) + print(f"\n--- Comparaison SOFA vs FreeFEM ---") + print(f" Erreur max : {err_max:.4e}") + print(f" Erreur RMS : {err_rms:.4e}") +else: + print(f"\n(Fichier SOFA '{OUTPUT_SOFA}' absent — comparaison désactivée)") + +# ─── FIGURE +fig, axes = plt.subplots(1, 2 if sofa_disponible else 1, + figsize=(13 if sofa_disponible else 7, 5)) +if not sofa_disponible: + axes = [axes] + +# Panneau 1 — Champ de déplacement +ax = axes[0] +ax.plot(x_exact, u_exact, '-', color='gray', lw=1.5, label='Solution exacte') +ax.plot(x_ff, u_ff, 'o-', color='#1a6faf', lw=2, markersize=4, + label=f'FreeFEM (N={N})') +if sofa_disponible: + ax.plot(x_sofa, u_sofa, 's--', color='#d62728', lw=2, markersize=4, + label='SOFA') +ax.set_xlabel("x (m)") +ax.set_ylabel("u(x) (m)") +ax.set_title("Déplacement axial — Barre 1D (charge uniforme)") +ax.legend() +ax.grid(True, alpha=0.3) + + diff --git a/beam_2d_freefem.py b/beam_2d_freefem.py new file mode 100644 index 00000000..7d4d36a5 --- /dev/null +++ b/beam_2d_freefem.py @@ -0,0 +1,114 @@ +import numpy as np +import pandas as pd +import subprocess +import tempfile +import os +import time + +# ========== PARAMÈTRES À DÉFINIR ========== +E = 21.5 # Module d'Young (GPa) +NU = 0.29 # Coefficient de Poisson +GRAVITY = -0.05 # Poids propre (N/m³) +N_LONG = 140 # Nombre d'éléments sur la longueur +N_HAUT = 35 # Nombre d'éléments sur la hauteur +OUTPUT_FILE = "freefem_results.txt" + +FREEFEM_CODE = """ +real E = {E}; +real sigma = {NU}; +real gravity = {GRAVITY}; +int n = {N_LONG}; +int m = {N_HAUT}; + +border a(t=2, 0){{x=0; y=t; label=1;}} +border b(t=0,10){{x=t; y=0; label=2;}} +border c(t=0, 2){{x=10; y=t; label=3;}} +border d(t=0,10){{x=10-t; y=2; label=4;}} + +mesh th = buildmesh(b(n) + c(m) + d(n) + a(m)); + +fespace Vh(th, [P1, P1]); +Vh [uu, vv], [w, s]; + +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)); + +solve Elasticity([uu,vv],[w,s], solver=sparsesolver) + = int2d(th)( + lambda * div(w,s) * div(uu,vv) + + 2.*mu * (epsilon(w,s)' * epsilon(uu,vv)) + ) + - int2d(th)(gravity*s) + + on(1, uu=0, vv=0); + +fespace Wh(th, P1); +Wh sigmaxx, sigmayy, sigmaxy, sigmavm; +sigmaxx = lambda*(dx(uu)+dy(vv)) + 2*mu*dx(uu); +sigmayy = lambda*(dx(uu)+dy(vv)) + 2*mu*dy(vv); +sigmaxy = mu*(dy(uu)+dx(vv)); +sigmavm = sqrt(sigmaxx^2 + sigmayy^2 - sigmaxx*sigmayy + 3*sigmaxy^2); + +cout << "Max uy = " << vv[].linfty << endl; +cout << "Max ux = " << uu[].linfty << endl; +cout << "sigmaxx max = " << sigmaxx[].max << endl; +cout << "sigmavm max = " << sigmavm[].max << endl; + + +{{ + ofstream fout("{OUTPUT_FILE}"); + fout << "x y ux uy sigmaxx sigmayy sigmaxy sigmavm" << endl; + for(int i=0; i= -0.01) & (ff["x"] <= 10.01) & + (ff["y"] >= -0.01) & (ff["y"] <= 2.01)].copy() + +# --- SOFA --- +# Tentative lecture espace-séparé puis virgule +try: + sf = pd.read_csv(FILE_SOFA, sep=r'\s+', comment='#') + if sf.shape[1] < 4: + raise ValueError +except Exception: + sf = pd.read_csv(FILE_SOFA, sep=',', comment='#') +sf.columns = [c.strip() for c in sf.columns] + +# Normalisation des noms de colonnes SOFA si besoin +col_map = {} +for c in sf.columns: + cl = c.lower().strip() + if cl in ("x",): col_map[c] = "x" + if cl in ("y",): col_map[c] = "y" + if cl in ("ux","u_x","disp_x"): col_map[c] = "ux" + if cl in ("uy","u_y","disp_y"): col_map[c] = "uy" +sf = sf.rename(columns=col_map) + +# Filtrage domaine SOFA (certains exports incluent z≠0) +sf = sf[(sf["x"] >= -0.01) & (sf["x"] <= 10.01) & + (sf["y"] >= -0.01) & (sf["y"] <= 2.01)].copy() + +print(f"\n FreeFEM : {len(ff):>6} nœuds | SOFA : {len(sf):>6} nœuds") +print(f" Colonnes FreeFEM : {list(ff.columns)}") +print(f" Colonnes SOFA : {list(sf.columns)}") + + +# GRILLE COMMUNE & INTERPOLATION + +gx, gy = np.mgrid[0:10:NX_GRID*1j, 0:2:NY_GRID*1j] +pts_grid = np.column_stack([gx.ravel(), gy.ravel()]) + +def interp_field(df, col, method="linear"): + """Interpolation griddata sur la grille commune.""" + pts = df[["x", "y"]].values + vals = df[col].values + return griddata(pts, vals, (gx, gy), method=method).reshape(gx.shape) + +# FreeFEM → linéaire (maillage dense non-structuré) +ff_g = {f: interp_field(ff, f, method="linear") for f in FIELDS} + +# SOFA → cubique si grille sparse, linéaire sinon +sofa_method = "linear" if len(sf) > 200 else "cubic" +sf_g = {f: interp_field(sf, f, method=sofa_method) for f in FIELDS} + +# Différences absolues et relatives +diff_g = {f: sf_g[f] - ff_g[f] for f in FIELDS} +rel_g = {} +for f in FIELDS: + ref = np.abs(ff_g[f]) + with np.errstate(invalid="ignore", divide="ignore"): + rel_g[f] = np.where(ref > 1e-12, + np.abs(diff_g[f]) / ref * 100.0, + np.nan) + +print(f"\n Interpolation SOFA : méthode='{sofa_method}'") +print(f" Grille commune : {NX_GRID}×{NY_GRID} = {NX_GRID*NY_GRID} points") + +# ══════════════════════════════════════════════════════════════════════════════ +# 3. FIGURE 1 — VISUALISATION DES MAILLAGES +# ══════════════════════════════════════════════════════════════════════════════ +fig1, axes1 = plt.subplots(1, 2, figsize=(15, 5)) +fig1.suptitle("Visualisation des maillages — FreeFEM vs SOFA", + fontsize=14, fontweight="bold") + +def plot_mesh(ax, df, title, color_node, color_tri): + x = df["x"].values + y = df["y"].values + # Triangulation de Delaunay pour visualiser la connectivité + if len(x) > 3: + tri = Delaunay(np.column_stack([x, y])) + # Filtrer les triangles trop grands (bord de domaine) + simp = tri.simplices + cx = x[simp].mean(axis=1) + cy = y[simp].mean(axis=1) + mask = (cx >= 0) & (cx <= 10) & (cy >= 0) & (cy <= 2) + # Taille max des arêtes + def edge_max(s): + pts = np.column_stack([x[s], y[s]]) + d = [np.linalg.norm(pts[i]-pts[j]) + for i,j in [(0,1),(1,2),(2,0)]] + return max(d) + edge_ok = np.array([edge_max(s) < 1.5 for s in simp]) + simp_ok = simp[mask & edge_ok] + triang = mtri.Triangulation(x, y, simp_ok) + ax.triplot(triang, color=color_tri, lw=0.4, alpha=0.6) + ax.scatter(x, y, s=3, c=color_node, zorder=3, alpha=0.8) + ax.set_xlim(-0.2, 10.2) + ax.set_ylim(-0.1, 2.1) + ax.set_aspect("equal") + ax.set_xlabel("x (m)") + ax.set_ylabel("y (m)") + ax.set_title(f"{title}\n({len(x)} nœuds)", fontsize=11) + ax.grid(True, alpha=0.2) + # Annotation densité + ax.text(0.98, 0.04, + f"Δx̄ ≈ {10/np.sqrt(len(x)):.3f} m", + transform=ax.transAxes, ha="right", fontsize=9, + bbox=dict(boxstyle="round,pad=0.3", fc="white", alpha=0.7)) + +plot_mesh(axes1[0], ff, "FreeFEM (maillage triangulaire non-structuré)", + "#1a6faf", "#6ab0de") +plot_mesh(axes1[1], sf, "SOFA (grille régulière + triangulation Q1→T3)", + "#d62728", "#f4a09a") + +plt.tight_layout() +plt.savefig("fig1_maillages.png", dpi=150, bbox_inches="tight") +print("\n ✔ fig1_maillages.png") + +# ══════════════════════════════════════════════════════════════════════════════ +# 4. FIGURE 2 — CARTES DE DÉPLACEMENTS (3 × 2) +# ══════════════════════════════════════════════════════════════════════════════ +fig2, axes2 = plt.subplots(3, 2, figsize=(15, 11)) +fig2.suptitle("Cartes de déplacements — FreeFEM | SOFA | Différence", + fontsize=14, fontweight="bold") + +rows = [ + ("FreeFEM", ff_g, "#1a6faf"), + ("SOFA", sf_g, "#d62728"), + ("Différence SOFA − FreeFEM", diff_g, None), +] + +for row_i, (rname, rdata, _) in enumerate(rows): + for col_i, field in enumerate(FIELDS): + ax = axes2[row_i, col_i] + dat = rdata[field] + if row_i < 2: + cmap = "RdBu_r" + vmax = np.nanpercentile(np.abs(dat), 99) + im = ax.contourf(gx, gy, dat, levels=50, + cmap=cmap, vmin=-vmax, vmax=vmax) + else: + cmap = "coolwarm" + vmax = np.nanpercentile(np.abs(dat), 97) + im = ax.contourf(gx, gy, dat, levels=50, + cmap=cmap, vmin=-vmax, vmax=vmax) + cb = fig2.colorbar(im, ax=ax, shrink=0.85, pad=0.02) + cb.ax.tick_params(labelsize=8) + ax.set_title(f"{rname} — {LABELS[field]}", fontsize=10) + ax.set_xlabel("x (m)") + ax.set_ylabel("y (m)") + ax.set_aspect("equal") + # Valeur extrême + vext = np.nanmax(np.abs(dat)) + ax.text(0.01, 0.04, f"max|·| = {vext:.4e}", + transform=ax.transAxes, fontsize=8, + bbox=dict(boxstyle="round,pad=0.2", fc="white", alpha=0.75)) + +plt.tight_layout() +plt.savefig("fig2_deplacements.png", dpi=150, bbox_inches="tight") +print(" ✔ fig2_deplacements.png") + +# ══════════════════════════════════════════════════════════════════════════════ +# 5. FIGURE 3 — PROFILS LIGNE MÉDIANE y = 1 +# ══════════════════════════════════════════════════════════════════════════════ +# Stratégie : on lit le profil depuis la GRILLE COMMUNE (déjà interpolée) +# pour les deux codes → pas de problème de nœuds manquants. +# On superpose aussi les nœuds bruts quand ils existent près de y=1. +# ── Tolérance adaptée à la densité de chaque maillage ────────────────────── +# FreeFEM : maillage dense (~4800 nœuds) → beaucoup de pts près de y=1 +# SOFA : grille 11×4 → les 4 lignes y sont à 0, 0.67, 1.33, 2 +# → aucun nœud à y=1 ! On prend la tolérance suffisante +# pour capturer y=0.67 ou y=1.33 (écart ≈ 0.33) +def tol_for(df, y_target=1.0, n_min=3): + """Trouve la tolérance minimale qui capture au moins n_min nœuds.""" + y_unique = np.sort(np.unique(np.round(df["y"].values, 6))) + nearest = y_unique[np.argmin(np.abs(y_unique - y_target))] + gap = np.abs(nearest - y_target) + 0.05 # petite marge + return max(gap, 0.12) + +xs_line = np.linspace(0, 10, 400) + +fig3, axes3 = plt.subplots(1, 2, figsize=(14, 5)) +fig3.suptitle("Profils de déplacement — ligne médiane y ≈ 1\n" + "(courbes depuis grille commune ; points = nœuds bruts)", + fontsize=12, fontweight="bold") + +for ax, field in zip(axes3, FIELDS): + + # ── A. Courbes depuis la grille commune (interpolée) ───────────────── + iy_common = np.argmin(np.abs(gy[0, :] - 1.0)) + xs_g = gx[:, iy_common] + + ax.plot(xs_g, ff_g[field][:, iy_common], + color="#1a6faf", lw=2.2, label="FreeFEM (grille commune)") + ax.plot(xs_g, sf_g[field][:, iy_common], + color="#d62728", lw=2.2, ls="--", label="SOFA (grille commune)") + + # ── B. Nœuds bruts FreeFEM (nombreux → scatter léger) ──────────────── + tol_ff = tol_for(ff) + sub_ff = ff[np.abs(ff["y"] - 1.0) < tol_ff].sort_values("x") + ax.scatter(sub_ff["x"], sub_ff[field], + s=6, c="#1a6faf", alpha=0.25, marker="o", + label=f"FreeFEM nœuds (|y−1|<{tol_ff:.2f})") + + # ── C. Nœuds bruts SOFA (peu nombreux → gros marqueurs) ────────────── + tol_sf = tol_for(sf) + sub_sf = sf[np.abs(sf["y"] - 1.0) < tol_sf].sort_values("x") + if len(sub_sf) > 0: + ax.scatter(sub_sf["x"], sub_sf[field], + s=60, c="#d62728", alpha=0.8, marker="s", zorder=5, + label=f"SOFA nœuds (|y−1|<{tol_sf:.2f}, n={len(sub_sf)})") + else: + ax.text(0.5, 0.5, "Aucun nœud SOFA proche de y=1\n(grille trop sparse)", + transform=ax.transAxes, ha="center", fontsize=9, + color="#d62728", alpha=0.7) + + # ── D. Différence sur axe secondaire ───────────────────────────────── + diff_line = diff_g[field][:, iy_common] + ax2 = ax.twinx() + ax2.fill_between(xs_g, 0, diff_line, + color="#2ca02c", alpha=0.15, label="Diff (aire)") + ax2.plot(xs_g, diff_line, color="#2ca02c", lw=1.3, + ls=":", label="Diff SOFA−FF") + ax2.axhline(0, color="#2ca02c", lw=0.5, ls="-", alpha=0.4) + ax2.set_ylabel("Δ = SOFA − FreeFEM", color="#2ca02c", fontsize=9) + ax2.tick_params(axis="y", labelcolor="#2ca02c") + + ax.set_xlabel("x (m)") + ax.set_ylabel(field) + ax.set_title(LABELS[field], fontsize=11) + ax.axvline(0, color="gray", lw=0.8, ls=":") + ax.grid(True, alpha=0.25) + lines_a, labs_a = ax.get_legend_handles_labels() + lines_b, labs_b = ax2.get_legend_handles_labels() + ax.legend(lines_a + lines_b, labs_a + labs_b, + fontsize=8, loc="lower left", ncol=2) + +plt.tight_layout() +plt.savefig("fig3_profil_mediane.png", dpi=150, bbox_inches="tight") +print(" ✔ fig3_profil_mediane.png") + +# ══════════════════════════════════════════════════════════════════════════════ +# 6. FIGURE 4 — CARTES D'ERREUR RELATIVE (%) +# ══════════════════════════════════════════════════════════════════════════════ +fig4, axes4 = plt.subplots(1, 2, figsize=(14, 4.5)) +fig4.suptitle("Erreur relative |SOFA − FreeFEM| / |FreeFEM| × 100 (%)", + fontsize=13, fontweight="bold") + +for ax, field in zip(axes4, FIELDS): + dat = rel_g[field] + vmax = np.nanpercentile(dat, 95) + im = ax.contourf(gx, gy, dat, levels=50, + cmap="hot_r", vmin=0, vmax=vmax) + cb = fig4.colorbar(im, ax=ax, shrink=0.9) + cb.set_label("%") + ax.set_title(f"Erreur relative — {LABELS[field]}", fontsize=10) + ax.set_xlabel("x (m)") + ax.set_ylabel("y (m)") + ax.set_aspect("equal") + med = np.nanmedian(dat) + mn = np.nanmean(dat) + ax.text(0.01, 0.04, + f"Médiane : {med:.2f}% Moyenne : {mn:.2f}%", + transform=ax.transAxes, fontsize=8, + bbox=dict(boxstyle="round,pad=0.2", fc="white", alpha=0.75)) + +plt.tight_layout() +plt.savefig("fig4_erreur_relative.png", dpi=150, bbox_inches="tight") +print(" ✔ fig4_erreur_relative.png") + +# ══════════════════════════════════════════════════════════════════════════════ +# 7. FIGURE 5 — HISTOGRAMMES DES ERREURS +# ══════════════════════════════════════════════════════════════════════════════ +fig5, axes5 = plt.subplots(1, 2, figsize=(12, 4)) +fig5.suptitle("Distribution des erreurs absolues SOFA − FreeFEM", + fontsize=13, fontweight="bold") + +for ax, field in zip(axes5, FIELDS): + d = diff_g[field].ravel() + d = d[~np.isnan(d)] + ax.hist(d, bins=80, color="#5a9fd4", edgecolor="white", + linewidth=0.3, alpha=0.85) + ax.axvline(np.mean(d), color="#d62728", lw=1.8, ls="-", + label=f"Moyenne : {np.mean(d):.2e}") + ax.axvline(np.median(d), color="#2ca02c", lw=1.8, ls="--", + label=f"Médiane : {np.median(d):.2e}") + ax.set_xlabel(f"Δ{field} = SOFA − FreeFEM") + ax.set_ylabel("Nombre de points") + ax.set_title(LABELS[field], fontsize=10) + ax.legend(fontsize=9) + ax.grid(True, alpha=0.3) + # Annotations stats + txt = (f"σ = {np.std(d):.2e}\n" + f"max|·| = {np.abs(d).max():.2e}\n" + f"RMSE = {np.sqrt(np.mean(d**2)):.2e}") + ax.text(0.97, 0.97, txt, transform=ax.transAxes, + va="top", ha="right", fontsize=8, + bbox=dict(boxstyle="round,pad=0.3", fc="lightyellow", alpha=0.8)) + +plt.tight_layout() +plt.savefig("fig5_histogrammes.png", dpi=150, bbox_inches="tight") +print(" ✔ fig5_histogrammes.png") + +# ══════════════════════════════════════════════════════════════════════════════ +# 8. TABLEAU RÉCAPITULATIF COMPLET +# ══════════════════════════════════════════════════════════════════════════════ +print("\n" + "═" * 72) +print(f" {'Quantité':<30} {'FreeFEM':>12} {'SOFA':>12} " + f"{'Δ absolu':>12} {'Δ relatif%':>11}") +print("═" * 72) + +metrics = { + # Déplacements verticaux + "max |u_y|" : lambda d: np.abs(d["uy"]).max(), + "min u_y (flèche)": lambda d: d["uy"].min(), + "mean u_y" : lambda d: d["uy"].mean(), + "RMSE u_y (grille)" : None, # traité séparément + # Déplacements horizontaux + "max |u_x|" : lambda d: np.abs(d["ux"]).max(), + "min u_x" : lambda d: d["ux"].min(), + "mean u_x" : lambda d: d["ux"].mean(), + "RMSE u_x (grille)" : None, + # Erreurs relatives médianes + "Erreur rel. méd. u_y (%)": None, + "Erreur rel. méd. u_x (%)": None, + "Erreur rel. moy. u_y (%)": None, + "Erreur rel. moy. u_x (%)": None, +} + +# Calcul RMSE sur grille commune +rmse = {} +for f in FIELDS: + d = diff_g[f].ravel() + d = d[~np.isnan(d)] + rmse[f] = np.sqrt(np.mean(d**2)) + +# Erreurs relatives médianes / moyennes sur grille +rel_med = {f: np.nanmedian(rel_g[f]) for f in FIELDS} +rel_moy = {f: np.nanmean(rel_g[f]) for f in FIELDS} + +def pct(vff, vsf): + return abs(vsf - vff) / (abs(vff) + 1e-15) * 100 + +rows_table = [ + ("── DÉPLACEMENT VERTICAL u_y ──────────────────────────────", None, None), + ("max |u_y|", np.abs(ff["uy"]).max(), np.abs(sf["uy"]).max()), + ("min u_y (flèche)", ff["uy"].min(), sf["uy"].min()), + ("mean u_y", ff["uy"].mean(), sf["uy"].mean()), + ("RMSE u_y (grille)", None, rmse["uy"]), + ("── DÉPLACEMENT HORIZONTAL u_x ────────────────────────────", None, None), + ("max |u_x|", np.abs(ff["ux"]).max(), np.abs(sf["ux"]).max()), + ("min u_x", ff["ux"].min(), sf["ux"].min()), + ("mean u_x", ff["ux"].mean(), sf["ux"].mean()), + ("RMSE u_x (grille)", None, rmse["ux"]), + ("── ERREURS RELATIVES SUR GRILLE COMMUNE ──────────────────", None, None), + ("Erreur rel. médiane u_y (%)", rel_med["uy"], None), + ("Erreur rel. médiane u_x (%)", rel_med["ux"], None), + ("Erreur rel. moyenne u_y (%)", rel_moy["uy"], None), + ("Erreur rel. moyenne u_x (%)", rel_moy["ux"], None), + ("── INFORMATIONS MAILLAGE ──────────────────────────────────", None, None), + ("Nœuds FreeFEM", len(ff), None), + ("Nœuds SOFA", len(sf), None), + ("Rapport nœuds FF/SOFA", len(ff)/max(len(sf),1), None), +] + +for row in rows_table: + label, vff, vsf = row + if vff is None and vsf is None: + print(f"\n {label}") + continue + if vsf is None: + # Valeur unique (erreur relative, info maillage) + print(f" {label:<45} {vff:>12.4f}") + continue + if vff is None: + # RMSE (pas de valeur FF individuelle) + print(f" {label:<45} {'—':>12} {vsf:>12.4e}") + continue + delta = vsf - vff + dp = pct(vff, vsf) + print(f" {label:<45} {vff:>12.6f} {vsf:>12.6f} " + f"{delta:>12.2e} {dp:>10.2f}%") + +print("═" * 72) + +# ══════════════════════════════════════════════════════════════════════════════ +# 9. EXPORT CSV DU TABLEAU +# ══════════════════════════════════════════════════════════════════════════════ +summary_rows = [] +for f in FIELDS: + ff_max = np.abs(ff[f]).max() + sf_max = np.abs(sf[f]).max() + ff_min = ff[f].min() + sf_min = sf[f].min() + summary_rows.append({ + "champ": f, + "FF_max_abs": ff_max, + "SF_max_abs": sf_max, + "delta_max_abs": sf_max - ff_max, + "pct_max_abs": pct(ff_max, sf_max), + "FF_min": ff_min, + "SF_min": sf_min, + "delta_min": sf_min - ff_min, + "pct_min": pct(ff_min, sf_min), + "RMSE_grille": rmse[f], + "err_rel_med_pct": rel_med[f], + "err_rel_moy_pct": rel_moy[f], + "noeuds_FF": len(ff), + "noeuds_SF": len(sf), + }) + +pd.DataFrame(summary_rows).to_csv("comparaison_summary.csv", index=False) +print("\n ✔ comparaison_summary.csv") + +# ══════════════════════════════════════════════════════════════════════════════ +# 10. AFFICHAGE +# ══════════════════════════════════════════════════════════════════════════════ +import os, subprocess, sys + +output_dir = os.path.dirname(os.path.abspath(__file__)) +figures = ["fig1_maillages.png", "fig2_deplacements.png", + "fig3_profil_mediane.png", "fig4_erreur_relative.png", + "fig5_histogrammes.png"] + +print("\n Fichiers générés :") +for fname in figures + ["comparaison_summary.csv"]: + fpath = os.path.join(output_dir, fname) + size = os.path.getsize(fpath) // 1024 if os.path.exists(fpath) else 0 + print(f" • {fname} ({size} Ko)") + +# Ouvrir automatiquement les images avec le viewer Windows/Mac/Linux +print("\n Ouverture des figures...") +for fname in figures: + fpath = os.path.join(output_dir, fname) + if not os.path.exists(fpath): + continue + try: + if sys.platform.startswith("win"): + os.startfile(fpath) # Windows : ouvre avec la visionneuse par défaut + elif sys.platform == "darwin": + subprocess.Popen(["open", fpath]) + else: + subprocess.Popen(["xdg-open", fpath]) + except Exception as e: + print(f" Impossible d'ouvrir {fname} : {e}") + + +try: + plt.show() +except Exception: + pass + + + + + + +# ══════════════════════════════════════════════════════════════════════════════ +# 11. FIGURE 6 — VISUALISATION DES NŒUDS AVEC LEURS NUMÉROS +# ══════════════════════════════════════════════════════════════════════════════ + +def plot_nodes_with_numbers(ax, df, title, color_node='blue', show_all=False, max_nodes_to_show=100): + """ + Affiche les nœuds avec leurs numéros. + + Parameters: + ----------- + ax : axes matplotlib + df : DataFrame avec colonnes x, y + title : str + color_node : couleur des nœuds + show_all : bool - si True affiche tous les nœuds, sinon limite à max_nodes_to_show + max_nodes_to_show : int - nombre maximum de nœuds à étiqueter + """ + x = df["x"].values + y = df["y"].values + n_nodes = len(x) + + # Tracer tous les nœuds + ax.scatter(x, y, s=20, c=color_node, zorder=3, alpha=0.7, edgecolors='black', linewidth=0.5) + + # Décider quels nœuds étiqueter + if show_all or n_nodes <= max_nodes_to_show: + # Afficher tous les nœuds + indices = range(n_nodes) + else: + # Afficher seulement quelques nœuds (bords + quelques internes) + # Identifier les nœuds uniques par ligne y + y_unique = np.unique(np.round(y, 6)) + n_y = len(y_unique) + n_x = n_nodes // n_y + + indices_to_label = [] + + # Ajouter les coins + corners = [0, n_x-1, n_nodes - n_x, n_nodes - 1] + indices_to_label.extend(corners) + + # Ajouter quelques nœuds sur les bords + step = max(1, n_x // 10) + for i in range(0, n_x, step): + indices_to_label.append(i) # bord bas + indices_to_label.append(n_nodes - n_x + i) # bord haut + + # Ajouter quelques nœuds internes + step_y = max(1, n_y // 5) + step_x = max(1, n_x // 5) + for iy in range(0, n_y, step_y): + for ix in range(0, n_x, step_x): + idx = iy * n_x + ix + if idx not in indices_to_label: + indices_to_label.append(idx) + + indices = list(set(indices_to_label)) # Éliminer les doublons + + # Afficher les numéros + for idx in indices: + if idx < n_nodes: + ax.annotate(str(idx), + (x[idx], y[idx]), + xytext=(5, 5), + textcoords='offset points', + fontsize=8, + bbox=dict(boxstyle="round,pad=0.2", fc="white", alpha=0.8, ec="gray"), + alpha=0.8) + + ax.set_xlim(-0.2, 10.2) + ax.set_ylim(-0.1, 2.1) + ax.set_aspect("equal") + ax.set_xlabel("x (m)") + ax.set_ylabel("y (m)") + ax.set_title(f"{title}\n({n_nodes} nœuds", fontsize=11) + ax.grid(True, alpha=0.3) + + # Ajouter une légende sur le nombre de nœuds affichés + if n_nodes > max_nodes_to_show: + ax.text(0.02, 0.98, f"Affichage de {len(indices)}/{n_nodes} nœuds", + transform=ax.transAxes, fontsize=8, + verticalalignment='top', + bbox=dict(boxstyle="round,pad=0.2", fc="yellow", alpha=0.7)) + + +# Créer la figure 6 +fig6, axes6 = plt.subplots(2, 2, figsize=(16, 12)) +fig6.suptitle("Visualisation des nœuds avec leurs numéros", + fontsize=14, fontweight="bold") + +# Sous-figure pour FreeFEM (tous les nœuds) +plot_nodes_with_numbers(axes6[0, 0], ff, "FreeFEM - Tous les nœuds", + color_node='#1a6faf', show_all=True) + +# Sous-figure pour FreeFEM (limité pour lisibilité) +plot_nodes_with_numbers(axes6[0, 1], ff, "FreeFEM - Échantillon de nœuds", + color_node='#1a6faf', show_all=False, max_nodes_to_show=100) + +# Sous-figure pour SOFA (tous les nœuds) +plot_nodes_with_numbers(axes6[1, 0], sf, "SOFA - Tous les nœuds", + color_node='#d62728', show_all=True) + +# Sous-figure pour SOFA (avec détails supplémentaires) +ax = axes6[1, 1] +x_sf = sf["x"].values +y_sf = sf["y"].values +ax.scatter(x_sf, y_sf, s=50, c='#d62728', zorder=3, alpha=0.8, edgecolors='black', linewidth=0.5) + +# Afficher tous les numéros pour SOFA (peu nombreux) +for idx, (xi, yi) in enumerate(zip(x_sf, y_sf)): + ax.annotate(str(idx), + (xi, yi), + xytext=(8, 8), + textcoords='offset points', + fontsize=10, + fontweight='bold', + bbox=dict(boxstyle="round,pad=0.3", fc="white", alpha=0.9, ec="red", lw=1.5), + alpha=0.9) + +# Ajouter les coordonnées pour quelques nœuds clés +for idx, (xi, yi) in enumerate(zip(x_sf, y_sf)): + if idx % 5 == 0 or idx < 4 or idx >= len(x_sf)-4: # Afficher les coordonnées pour certains nœuds + ax.annotate(f"({xi:.1f}, {yi:.1f})", + (xi, yi), + xytext=(8, -15), + textcoords='offset points', + fontsize=7, + alpha=0.7) + +ax.set_xlim(-0.2, 10.2) +ax.set_ylim(-0.1, 2.1) +ax.set_aspect("equal") +ax.set_xlabel("x (m)") +ax.set_ylabel("y (m)") +ax.set_title(f"SOFA - Détail des nœuds\n({len(sf)} nœuds, grille régulière)", fontsize=11) +ax.grid(True, alpha=0.3) + +plt.tight_layout() +plt.savefig("fig6_noeuds_avec_numeros.png", dpi=150, bbox_inches="tight") +print(" ✔ fig6_noeuds_avec_numeros.png") + + +# ══════════════════════════════════════════════════════════════════════════════ +# 12. FIGURE 7 — TOPOLOGIE DU MAILLAGE SOFA (grille structurée) +# ══════════════════════════════════════════════════════════════════════════════ + +def plot_sofa_grid_structure(ax, df, nx, ny): + """Visualise la structure de grille de SOFA avec les connexions""" + x = df["x"].values + y = df["y"].values + + # Déterminer nx, ny à partir des données si non fournis + if nx is None or ny is None: + y_unique = np.unique(np.round(y, 6)) + ny = len(y_unique) + nx = len(x) // ny + + # Tracer les lignes de la grille + for i in range(ny): + row_indices = slice(i*nx, (i+1)*nx) + ax.plot(x[row_indices], y[row_indices], 'b-', linewidth=1, alpha=0.5, color='gray') + + for j in range(nx): + col_indices = list(range(j, nx*ny, nx)) + ax.plot(x[col_indices], y[col_indices], 'b-', linewidth=1, alpha=0.5, color='gray') + + # Tracer les nœuds + ax.scatter(x, y, s=30, c='#d62728', zorder=3, alpha=0.8, edgecolors='black') + + # Numéroter les nœuds + for idx, (xi, yi) in enumerate(zip(x, y)): + ax.annotate(str(idx), + (xi, yi), + xytext=(5, 5), + textcoords='offset points', + fontsize=9, + bbox=dict(boxstyle="round,pad=0.2", fc="white", alpha=0.8)) + + ax.set_xlim(-0.2, 10.2) + ax.set_ylim(-0.1, 2.1) + ax.set_aspect("equal") + ax.set_xlabel("x (m)") + ax.set_ylabel("y (m)") + ax.set_title(f"Structure de la grille SOFA\n{nx} × {ny} = {nx*ny} nœuds", fontsize=11) + ax.grid(True, alpha=0.3) + +# Créer la figure 7 +fig7, ax7 = plt.subplots(1, 1, figsize=(12, 6)) + +# Déterminer nx, ny à partir des données SOFA +y_unique_sf = np.unique(np.round(sf["y"].values, 6)) +ny_sf = len(y_unique_sf) +nx_sf = len(sf) // ny_sf + +plot_sofa_grid_structure(ax7, sf, nx_sf, ny_sf) + +plt.tight_layout() +plt.savefig("fig7_grille_sofa.png", dpi=150, bbox_inches="tight") +print(" ✔ fig7_grille_sofa.png") + + +# ══════════════════════════════════════════════════════════════════════════════ +# 13. FIGURE 8 — CARTE DE CHALEUR DE LA DENSITÉ DE NŒUDS (FreeFEM) +# ══════════════════════════════════════════════════════════════════════════════ + +def plot_node_density_heatmap(ax, df, title, bins=50): + """Affiche une carte de chaleur de la densité de nœuds""" + x = df["x"].values + y = df["y"].values + + # Créer l'histogramme 2D + H, xedges, yedges = np.histogram2d(x, y, bins=bins, range=[[0, 10], [0, 2]]) + + # Afficher la carte de chaleur + im = ax.imshow(H.T, origin='lower', extent=[0, 10, 0, 2], + cmap='hot', aspect='auto', alpha=0.8) + + # Ajouter la colorbar + plt.colorbar(im, ax=ax, label='Nombre de nœuds par cellule') + + ax.set_xlabel("x (m)") + ax.set_ylabel("y (m)") + ax.set_title(title) + ax.grid(True, alpha=0.2) + +fig8, axes8 = plt.subplots(1, 2, figsize=(14, 5)) +fig8.suptitle("Densité de nœuds dans le maillage", fontsize=13, fontweight="bold") + +plot_node_density_heatmap(axes8[0], ff, "FreeFEM - Distribution des nœuds", bins=40) +plot_node_density_heatmap(axes8[1], sf, "SOFA - Distribution des nœuds (grille régulière)", bins=20) + +plt.tight_layout() +plt.savefig("fig8_densite_noeuds.png", dpi=150, bbox_inches="tight") +print(" ✔ fig8_densite_noeuds.png") + + +# ══════════════════════════════════════════════════════════════════════════════ +# Mettre à jour la liste des figures +figures += ["fig6_noeuds_avec_numeros.png", "fig7_grille_sofa.png", "fig8_densite_noeuds.png"] + + + + + + +# prise de note + +""" +POURQUOI ON PEUT COMPARER LES DEUX CODES MALGRÉ DES MAILLAGES DIFFÉRENTS ? +=========================================================================== + +PROBLÈME DE DÉPART : + • FreeFEM génère un maillage triangulaire NON-STRUCTURÉ avec ~4800 nœuds. + Les nœuds sont placés de façon irrégulière, denses là où le gradient est + fort (encastrement gauche), plus rares au centre. + + • SOFA utilise une RegularGridTopology 11×4 = 44 nœuds seulement, + placés sur une grille parfaitement régulière. + Les lignes y sont : 0.0, 0.67, 1.33, 2.0 + → aucun nœud à y = 1.0 exactement ! + +SOLUTION — LA GRILLE COMMUNE : + 1. On crée une grille uniforme 300×80 = 24 000 points couvrant [0,10]×[0,2]. + 2. On interpole chaque champ (ux, uy) sur cette grille depuis les nœuds bruts : + - FreeFEM : scipy.griddata(..., method='linear') + → triangulation Delaunay + interpolation linéaire par morceaux (P1) + → très fidèle car maillage dense + - SOFA : scipy.griddata(..., method='cubic') + → interpolation cubique car 44 points sont trop peu pour 'linear' + → lisse le champ entre les nœuds réguliers + + 3. Une fois les deux champs sur la même grille, on peut : + diff[i,j] = SOFA_interpolé[i,j] - FreeFEM_interpolé[i,j] + et calculer RMSE, erreur relative, histogrammes, etc. + +POURQUOI LA COMPARAISON A DU SENS ? + Les deux codes résolvent la MÊME équation (élasticité linéaire plane, + même E, nu, gravité, encastrement gauche). La différence vient de : + a) La formulation : plane_strain (Vec3, λ=Eν/((1+ν)(1-2ν))) vs + les coefficients Lamé 3D de FreeFEM (identiques à plane_strain) + b) La discrétisation : triangles P1 vs Q1 subdivisé en triangles + c) Le nombre de DDL : 44 nœuds SOFA << 4814 nœuds FreeFEM + → SOFA est moins précis, mais donne la bonne tendance physique. + + Le tableau récapitulatif vous donne l'écart % entre les deux pour + chaque quantité d'intérêt (flèche maximale, RMSE, erreur médiane...). +""" \ No newline at end of file From 46509007d3b69dd085442744ba056c408b0f7c94 Mon Sep 17 00:00:00 2001 From: Fimache Date: Fri, 27 Mar 2026 20:08:46 +0100 Subject: [PATCH 3/5] Convert XML scenes to Python scripts in examples/Freefem/sofa/ - Convert beam-2d.scn to beam-2d.py - Convert barre-traction.scn to barre-traction.py - Place scripts in examples/Freefem/sofa/ - Remove comparison scripts and duplicate .edp files (to be added in next PR) --- barre_1d_pyfreefem.py | 91 -- beam_2d_freefem.py | 114 --- comparaison.py | 796 ------------------ examples/Freefem/3d-cube.edp | 78 -- examples/Freefem/barre-1d-solution-man.edp | 44 - examples/Freefem/barre_exacte.edp | 42 - examples/Freefem/beam-2d.edp | 104 --- examples/Freefem/elasticity_3D.edp | 105 --- .../Freefem/sofa/barre-traction.py | 0 .../Freefem/sofa/beam-2d.py | 0 10 files changed, 1374 deletions(-) delete mode 100644 barre_1d_pyfreefem.py delete mode 100644 beam_2d_freefem.py delete mode 100644 comparaison.py delete mode 100644 examples/Freefem/3d-cube.edp delete mode 100644 examples/Freefem/barre-1d-solution-man.edp delete mode 100644 examples/Freefem/barre_exacte.edp delete mode 100644 examples/Freefem/beam-2d.edp delete mode 100644 examples/Freefem/elasticity_3D.edp rename barre_1d_sofa.py => examples/Freefem/sofa/barre-traction.py (100%) rename beam_2d_sofa.py => examples/Freefem/sofa/beam-2d.py (100%) diff --git a/barre_1d_pyfreefem.py b/barre_1d_pyfreefem.py deleted file mode 100644 index 7c5d4d8d..00000000 --- a/barre_1d_pyfreefem.py +++ /dev/null @@ -1,91 +0,0 @@ -import numpy as np -import pandas as pd -import matplotlib.pyplot as plt -import subprocess -import tempfile -import os - -# PARAMÈTRES -L = 1.0 -E = 1e6 -A = 0.001 -N = 64 -OUTPUT_FF = "freefem_deplacement.txt" - - -FREEFEM_CODE = """ -real L = {L}; -real E = {E}; -real A = {A}; -int n = {N}; - -meshL Th = segment(n, [x*L]); -fespace Vh(Th, P1); -Vh u, v; - -problem Traction(u, v) = - int1d(Th)( E * A * dx(u) * dx(v) ) - - int1d(Th)( 1.0 * v ) - + on(1, u=0); - -Traction; - -ofstream file("{OUTPUT_FF}"); -file << "x u" << endl; -for (int i = 0; i < Th.nv; i++) {{ - real xi = Th(i).x; - real ui = u(xi, 0); - file << xi << " " << ui << endl; -}} - -cout << "Export OK : " << Th.nv << " noeuds ecrits dans {OUTPUT_FF}" << endl; -""".format(L=L, E=E, A=A, N=N, OUTPUT_FF=OUTPUT_FF) - -# EXÉCUTION FREEFEM -print("Lancement FreeFEM++...") -with tempfile.NamedTemporaryFile(mode='w', suffix='.edp', - delete=False, encoding='utf-8') as f: - f.write(FREEFEM_CODE) - temp_file = f.name - -try: - result = subprocess.run( - ['FreeFem++', temp_file], - capture_output=True, text=True, encoding='utf-8' - ) - print(result.stdout) - if result.returncode != 0: - print("ERREUR FreeFEM :\n", result.stderr) - exit(1) -finally: - os.unlink(temp_file) - -# LECTURE RÉSULTATS -if not os.path.exists(OUTPUT_FF): - print(f"Fichier {OUTPUT_FF} non trouvé — vérifier FreeFEM.") - exit(1) - -df_ff = pd.read_csv(OUTPUT_FF, sep=r'\s+') -df_ff = df_ff.sort_values('x').reset_index(drop=True) -x_ff = df_ff['x'].values -u_ff = df_ff['u'].values -#print(f"\nFreeFEM : {len(x_ff)} nœuds chargés depuis '{OUTPUT_FF}'") -#print(df_ff.to_string(index=False)) - -# SOLUTION EXACTE -# -EA u'' = 1, u(0)=0, u'(L)=0 → u(x) = x*(L - x/2) / (EA) -x_exact = np.linspace(0, L, 500) -u_exact = x_exact * (L - x_exact / 2.0) / (E * A) - -# ─── FIGURE -fig, ax = plt.subplots(figsize=(7, 5)) - -ax.plot(x_exact, u_exact, '-', color='gray', lw=1.5, label='Solution exacte') -ax.plot(x_ff, u_ff, 'o-', color='#1a6faf', lw=2, markersize=4, - label=f'FreeFEM (N={N})') -ax.set_xlabel("x (m)") -ax.set_ylabel("u(x) (m)") -ax.set_title("Déplacement axial — Barre 1D ") -ax.legend() -ax.grid(True, alpha=0.3) - diff --git a/beam_2d_freefem.py b/beam_2d_freefem.py deleted file mode 100644 index 7d4d36a5..00000000 --- a/beam_2d_freefem.py +++ /dev/null @@ -1,114 +0,0 @@ -import numpy as np -import pandas as pd -import subprocess -import tempfile -import os -import time - -# ========== PARAMÈTRES À DÉFINIR ========== -E = 21.5 # Module d'Young (GPa) -NU = 0.29 # Coefficient de Poisson -GRAVITY = -0.05 # Poids propre (N/m³) -N_LONG = 140 # Nombre d'éléments sur la longueur -N_HAUT = 35 # Nombre d'éléments sur la hauteur -OUTPUT_FILE = "freefem_results.txt" - -FREEFEM_CODE = """ -real E = {E}; -real sigma = {NU}; -real gravity = {GRAVITY}; -int n = {N_LONG}; -int m = {N_HAUT}; - -border a(t=2, 0){{x=0; y=t; label=1;}} -border b(t=0,10){{x=t; y=0; label=2;}} -border c(t=0, 2){{x=10; y=t; label=3;}} -border d(t=0,10){{x=10-t; y=2; label=4;}} - -mesh th = buildmesh(b(n) + c(m) + d(n) + a(m)); - -fespace Vh(th, [P1, P1]); -Vh [uu, vv], [w, s]; - -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)); - -solve Elasticity([uu,vv],[w,s], solver=sparsesolver) - = int2d(th)( - lambda * div(w,s) * div(uu,vv) - + 2.*mu * (epsilon(w,s)' * epsilon(uu,vv)) - ) - - int2d(th)(gravity*s) - + on(1, uu=0, vv=0); - -fespace Wh(th, P1); -Wh sigmaxx, sigmayy, sigmaxy, sigmavm; -sigmaxx = lambda*(dx(uu)+dy(vv)) + 2*mu*dx(uu); -sigmayy = lambda*(dx(uu)+dy(vv)) + 2*mu*dy(vv); -sigmaxy = mu*(dy(uu)+dx(vv)); -sigmavm = sqrt(sigmaxx^2 + sigmayy^2 - sigmaxx*sigmayy + 3*sigmaxy^2); - -cout << "Max uy = " << vv[].linfty << endl; -cout << "Max ux = " << uu[].linfty << endl; -cout << "sigmaxx max = " << sigmaxx[].max << endl; -cout << "sigmavm max = " << sigmavm[].max << endl; - - -{{ - ofstream fout("{OUTPUT_FILE}"); - fout << "x y ux uy sigmaxx sigmayy sigmaxy sigmavm" << endl; - for(int i=0; i= -0.01) & (ff["x"] <= 10.01) & - (ff["y"] >= -0.01) & (ff["y"] <= 2.01)].copy() - -# --- SOFA --- -# Tentative lecture espace-séparé puis virgule -try: - sf = pd.read_csv(FILE_SOFA, sep=r'\s+', comment='#') - if sf.shape[1] < 4: - raise ValueError -except Exception: - sf = pd.read_csv(FILE_SOFA, sep=',', comment='#') -sf.columns = [c.strip() for c in sf.columns] - -# Normalisation des noms de colonnes SOFA si besoin -col_map = {} -for c in sf.columns: - cl = c.lower().strip() - if cl in ("x",): col_map[c] = "x" - if cl in ("y",): col_map[c] = "y" - if cl in ("ux","u_x","disp_x"): col_map[c] = "ux" - if cl in ("uy","u_y","disp_y"): col_map[c] = "uy" -sf = sf.rename(columns=col_map) - -# Filtrage domaine SOFA (certains exports incluent z≠0) -sf = sf[(sf["x"] >= -0.01) & (sf["x"] <= 10.01) & - (sf["y"] >= -0.01) & (sf["y"] <= 2.01)].copy() - -print(f"\n FreeFEM : {len(ff):>6} nœuds | SOFA : {len(sf):>6} nœuds") -print(f" Colonnes FreeFEM : {list(ff.columns)}") -print(f" Colonnes SOFA : {list(sf.columns)}") - - -# GRILLE COMMUNE & INTERPOLATION - -gx, gy = np.mgrid[0:10:NX_GRID*1j, 0:2:NY_GRID*1j] -pts_grid = np.column_stack([gx.ravel(), gy.ravel()]) - -def interp_field(df, col, method="linear"): - """Interpolation griddata sur la grille commune.""" - pts = df[["x", "y"]].values - vals = df[col].values - return griddata(pts, vals, (gx, gy), method=method).reshape(gx.shape) - -# FreeFEM → linéaire (maillage dense non-structuré) -ff_g = {f: interp_field(ff, f, method="linear") for f in FIELDS} - -# SOFA → cubique si grille sparse, linéaire sinon -sofa_method = "linear" if len(sf) > 200 else "cubic" -sf_g = {f: interp_field(sf, f, method=sofa_method) for f in FIELDS} - -# Différences absolues et relatives -diff_g = {f: sf_g[f] - ff_g[f] for f in FIELDS} -rel_g = {} -for f in FIELDS: - ref = np.abs(ff_g[f]) - with np.errstate(invalid="ignore", divide="ignore"): - rel_g[f] = np.where(ref > 1e-12, - np.abs(diff_g[f]) / ref * 100.0, - np.nan) - -print(f"\n Interpolation SOFA : méthode='{sofa_method}'") -print(f" Grille commune : {NX_GRID}×{NY_GRID} = {NX_GRID*NY_GRID} points") - -# ══════════════════════════════════════════════════════════════════════════════ -# 3. FIGURE 1 — VISUALISATION DES MAILLAGES -# ══════════════════════════════════════════════════════════════════════════════ -fig1, axes1 = plt.subplots(1, 2, figsize=(15, 5)) -fig1.suptitle("Visualisation des maillages — FreeFEM vs SOFA", - fontsize=14, fontweight="bold") - -def plot_mesh(ax, df, title, color_node, color_tri): - x = df["x"].values - y = df["y"].values - # Triangulation de Delaunay pour visualiser la connectivité - if len(x) > 3: - tri = Delaunay(np.column_stack([x, y])) - # Filtrer les triangles trop grands (bord de domaine) - simp = tri.simplices - cx = x[simp].mean(axis=1) - cy = y[simp].mean(axis=1) - mask = (cx >= 0) & (cx <= 10) & (cy >= 0) & (cy <= 2) - # Taille max des arêtes - def edge_max(s): - pts = np.column_stack([x[s], y[s]]) - d = [np.linalg.norm(pts[i]-pts[j]) - for i,j in [(0,1),(1,2),(2,0)]] - return max(d) - edge_ok = np.array([edge_max(s) < 1.5 for s in simp]) - simp_ok = simp[mask & edge_ok] - triang = mtri.Triangulation(x, y, simp_ok) - ax.triplot(triang, color=color_tri, lw=0.4, alpha=0.6) - ax.scatter(x, y, s=3, c=color_node, zorder=3, alpha=0.8) - ax.set_xlim(-0.2, 10.2) - ax.set_ylim(-0.1, 2.1) - ax.set_aspect("equal") - ax.set_xlabel("x (m)") - ax.set_ylabel("y (m)") - ax.set_title(f"{title}\n({len(x)} nœuds)", fontsize=11) - ax.grid(True, alpha=0.2) - # Annotation densité - ax.text(0.98, 0.04, - f"Δx̄ ≈ {10/np.sqrt(len(x)):.3f} m", - transform=ax.transAxes, ha="right", fontsize=9, - bbox=dict(boxstyle="round,pad=0.3", fc="white", alpha=0.7)) - -plot_mesh(axes1[0], ff, "FreeFEM (maillage triangulaire non-structuré)", - "#1a6faf", "#6ab0de") -plot_mesh(axes1[1], sf, "SOFA (grille régulière + triangulation Q1→T3)", - "#d62728", "#f4a09a") - -plt.tight_layout() -plt.savefig("fig1_maillages.png", dpi=150, bbox_inches="tight") -print("\n ✔ fig1_maillages.png") - -# ══════════════════════════════════════════════════════════════════════════════ -# 4. FIGURE 2 — CARTES DE DÉPLACEMENTS (3 × 2) -# ══════════════════════════════════════════════════════════════════════════════ -fig2, axes2 = plt.subplots(3, 2, figsize=(15, 11)) -fig2.suptitle("Cartes de déplacements — FreeFEM | SOFA | Différence", - fontsize=14, fontweight="bold") - -rows = [ - ("FreeFEM", ff_g, "#1a6faf"), - ("SOFA", sf_g, "#d62728"), - ("Différence SOFA − FreeFEM", diff_g, None), -] - -for row_i, (rname, rdata, _) in enumerate(rows): - for col_i, field in enumerate(FIELDS): - ax = axes2[row_i, col_i] - dat = rdata[field] - if row_i < 2: - cmap = "RdBu_r" - vmax = np.nanpercentile(np.abs(dat), 99) - im = ax.contourf(gx, gy, dat, levels=50, - cmap=cmap, vmin=-vmax, vmax=vmax) - else: - cmap = "coolwarm" - vmax = np.nanpercentile(np.abs(dat), 97) - im = ax.contourf(gx, gy, dat, levels=50, - cmap=cmap, vmin=-vmax, vmax=vmax) - cb = fig2.colorbar(im, ax=ax, shrink=0.85, pad=0.02) - cb.ax.tick_params(labelsize=8) - ax.set_title(f"{rname} — {LABELS[field]}", fontsize=10) - ax.set_xlabel("x (m)") - ax.set_ylabel("y (m)") - ax.set_aspect("equal") - # Valeur extrême - vext = np.nanmax(np.abs(dat)) - ax.text(0.01, 0.04, f"max|·| = {vext:.4e}", - transform=ax.transAxes, fontsize=8, - bbox=dict(boxstyle="round,pad=0.2", fc="white", alpha=0.75)) - -plt.tight_layout() -plt.savefig("fig2_deplacements.png", dpi=150, bbox_inches="tight") -print(" ✔ fig2_deplacements.png") - -# ══════════════════════════════════════════════════════════════════════════════ -# 5. FIGURE 3 — PROFILS LIGNE MÉDIANE y = 1 -# ══════════════════════════════════════════════════════════════════════════════ -# Stratégie : on lit le profil depuis la GRILLE COMMUNE (déjà interpolée) -# pour les deux codes → pas de problème de nœuds manquants. -# On superpose aussi les nœuds bruts quand ils existent près de y=1. -# ── Tolérance adaptée à la densité de chaque maillage ────────────────────── -# FreeFEM : maillage dense (~4800 nœuds) → beaucoup de pts près de y=1 -# SOFA : grille 11×4 → les 4 lignes y sont à 0, 0.67, 1.33, 2 -# → aucun nœud à y=1 ! On prend la tolérance suffisante -# pour capturer y=0.67 ou y=1.33 (écart ≈ 0.33) -def tol_for(df, y_target=1.0, n_min=3): - """Trouve la tolérance minimale qui capture au moins n_min nœuds.""" - y_unique = np.sort(np.unique(np.round(df["y"].values, 6))) - nearest = y_unique[np.argmin(np.abs(y_unique - y_target))] - gap = np.abs(nearest - y_target) + 0.05 # petite marge - return max(gap, 0.12) - -xs_line = np.linspace(0, 10, 400) - -fig3, axes3 = plt.subplots(1, 2, figsize=(14, 5)) -fig3.suptitle("Profils de déplacement — ligne médiane y ≈ 1\n" - "(courbes depuis grille commune ; points = nœuds bruts)", - fontsize=12, fontweight="bold") - -for ax, field in zip(axes3, FIELDS): - - # ── A. Courbes depuis la grille commune (interpolée) ───────────────── - iy_common = np.argmin(np.abs(gy[0, :] - 1.0)) - xs_g = gx[:, iy_common] - - ax.plot(xs_g, ff_g[field][:, iy_common], - color="#1a6faf", lw=2.2, label="FreeFEM (grille commune)") - ax.plot(xs_g, sf_g[field][:, iy_common], - color="#d62728", lw=2.2, ls="--", label="SOFA (grille commune)") - - # ── B. Nœuds bruts FreeFEM (nombreux → scatter léger) ──────────────── - tol_ff = tol_for(ff) - sub_ff = ff[np.abs(ff["y"] - 1.0) < tol_ff].sort_values("x") - ax.scatter(sub_ff["x"], sub_ff[field], - s=6, c="#1a6faf", alpha=0.25, marker="o", - label=f"FreeFEM nœuds (|y−1|<{tol_ff:.2f})") - - # ── C. Nœuds bruts SOFA (peu nombreux → gros marqueurs) ────────────── - tol_sf = tol_for(sf) - sub_sf = sf[np.abs(sf["y"] - 1.0) < tol_sf].sort_values("x") - if len(sub_sf) > 0: - ax.scatter(sub_sf["x"], sub_sf[field], - s=60, c="#d62728", alpha=0.8, marker="s", zorder=5, - label=f"SOFA nœuds (|y−1|<{tol_sf:.2f}, n={len(sub_sf)})") - else: - ax.text(0.5, 0.5, "Aucun nœud SOFA proche de y=1\n(grille trop sparse)", - transform=ax.transAxes, ha="center", fontsize=9, - color="#d62728", alpha=0.7) - - # ── D. Différence sur axe secondaire ───────────────────────────────── - diff_line = diff_g[field][:, iy_common] - ax2 = ax.twinx() - ax2.fill_between(xs_g, 0, diff_line, - color="#2ca02c", alpha=0.15, label="Diff (aire)") - ax2.plot(xs_g, diff_line, color="#2ca02c", lw=1.3, - ls=":", label="Diff SOFA−FF") - ax2.axhline(0, color="#2ca02c", lw=0.5, ls="-", alpha=0.4) - ax2.set_ylabel("Δ = SOFA − FreeFEM", color="#2ca02c", fontsize=9) - ax2.tick_params(axis="y", labelcolor="#2ca02c") - - ax.set_xlabel("x (m)") - ax.set_ylabel(field) - ax.set_title(LABELS[field], fontsize=11) - ax.axvline(0, color="gray", lw=0.8, ls=":") - ax.grid(True, alpha=0.25) - lines_a, labs_a = ax.get_legend_handles_labels() - lines_b, labs_b = ax2.get_legend_handles_labels() - ax.legend(lines_a + lines_b, labs_a + labs_b, - fontsize=8, loc="lower left", ncol=2) - -plt.tight_layout() -plt.savefig("fig3_profil_mediane.png", dpi=150, bbox_inches="tight") -print(" ✔ fig3_profil_mediane.png") - -# ══════════════════════════════════════════════════════════════════════════════ -# 6. FIGURE 4 — CARTES D'ERREUR RELATIVE (%) -# ══════════════════════════════════════════════════════════════════════════════ -fig4, axes4 = plt.subplots(1, 2, figsize=(14, 4.5)) -fig4.suptitle("Erreur relative |SOFA − FreeFEM| / |FreeFEM| × 100 (%)", - fontsize=13, fontweight="bold") - -for ax, field in zip(axes4, FIELDS): - dat = rel_g[field] - vmax = np.nanpercentile(dat, 95) - im = ax.contourf(gx, gy, dat, levels=50, - cmap="hot_r", vmin=0, vmax=vmax) - cb = fig4.colorbar(im, ax=ax, shrink=0.9) - cb.set_label("%") - ax.set_title(f"Erreur relative — {LABELS[field]}", fontsize=10) - ax.set_xlabel("x (m)") - ax.set_ylabel("y (m)") - ax.set_aspect("equal") - med = np.nanmedian(dat) - mn = np.nanmean(dat) - ax.text(0.01, 0.04, - f"Médiane : {med:.2f}% Moyenne : {mn:.2f}%", - transform=ax.transAxes, fontsize=8, - bbox=dict(boxstyle="round,pad=0.2", fc="white", alpha=0.75)) - -plt.tight_layout() -plt.savefig("fig4_erreur_relative.png", dpi=150, bbox_inches="tight") -print(" ✔ fig4_erreur_relative.png") - -# ══════════════════════════════════════════════════════════════════════════════ -# 7. FIGURE 5 — HISTOGRAMMES DES ERREURS -# ══════════════════════════════════════════════════════════════════════════════ -fig5, axes5 = plt.subplots(1, 2, figsize=(12, 4)) -fig5.suptitle("Distribution des erreurs absolues SOFA − FreeFEM", - fontsize=13, fontweight="bold") - -for ax, field in zip(axes5, FIELDS): - d = diff_g[field].ravel() - d = d[~np.isnan(d)] - ax.hist(d, bins=80, color="#5a9fd4", edgecolor="white", - linewidth=0.3, alpha=0.85) - ax.axvline(np.mean(d), color="#d62728", lw=1.8, ls="-", - label=f"Moyenne : {np.mean(d):.2e}") - ax.axvline(np.median(d), color="#2ca02c", lw=1.8, ls="--", - label=f"Médiane : {np.median(d):.2e}") - ax.set_xlabel(f"Δ{field} = SOFA − FreeFEM") - ax.set_ylabel("Nombre de points") - ax.set_title(LABELS[field], fontsize=10) - ax.legend(fontsize=9) - ax.grid(True, alpha=0.3) - # Annotations stats - txt = (f"σ = {np.std(d):.2e}\n" - f"max|·| = {np.abs(d).max():.2e}\n" - f"RMSE = {np.sqrt(np.mean(d**2)):.2e}") - ax.text(0.97, 0.97, txt, transform=ax.transAxes, - va="top", ha="right", fontsize=8, - bbox=dict(boxstyle="round,pad=0.3", fc="lightyellow", alpha=0.8)) - -plt.tight_layout() -plt.savefig("fig5_histogrammes.png", dpi=150, bbox_inches="tight") -print(" ✔ fig5_histogrammes.png") - -# ══════════════════════════════════════════════════════════════════════════════ -# 8. TABLEAU RÉCAPITULATIF COMPLET -# ══════════════════════════════════════════════════════════════════════════════ -print("\n" + "═" * 72) -print(f" {'Quantité':<30} {'FreeFEM':>12} {'SOFA':>12} " - f"{'Δ absolu':>12} {'Δ relatif%':>11}") -print("═" * 72) - -metrics = { - # Déplacements verticaux - "max |u_y|" : lambda d: np.abs(d["uy"]).max(), - "min u_y (flèche)": lambda d: d["uy"].min(), - "mean u_y" : lambda d: d["uy"].mean(), - "RMSE u_y (grille)" : None, # traité séparément - # Déplacements horizontaux - "max |u_x|" : lambda d: np.abs(d["ux"]).max(), - "min u_x" : lambda d: d["ux"].min(), - "mean u_x" : lambda d: d["ux"].mean(), - "RMSE u_x (grille)" : None, - # Erreurs relatives médianes - "Erreur rel. méd. u_y (%)": None, - "Erreur rel. méd. u_x (%)": None, - "Erreur rel. moy. u_y (%)": None, - "Erreur rel. moy. u_x (%)": None, -} - -# Calcul RMSE sur grille commune -rmse = {} -for f in FIELDS: - d = diff_g[f].ravel() - d = d[~np.isnan(d)] - rmse[f] = np.sqrt(np.mean(d**2)) - -# Erreurs relatives médianes / moyennes sur grille -rel_med = {f: np.nanmedian(rel_g[f]) for f in FIELDS} -rel_moy = {f: np.nanmean(rel_g[f]) for f in FIELDS} - -def pct(vff, vsf): - return abs(vsf - vff) / (abs(vff) + 1e-15) * 100 - -rows_table = [ - ("── DÉPLACEMENT VERTICAL u_y ──────────────────────────────", None, None), - ("max |u_y|", np.abs(ff["uy"]).max(), np.abs(sf["uy"]).max()), - ("min u_y (flèche)", ff["uy"].min(), sf["uy"].min()), - ("mean u_y", ff["uy"].mean(), sf["uy"].mean()), - ("RMSE u_y (grille)", None, rmse["uy"]), - ("── DÉPLACEMENT HORIZONTAL u_x ────────────────────────────", None, None), - ("max |u_x|", np.abs(ff["ux"]).max(), np.abs(sf["ux"]).max()), - ("min u_x", ff["ux"].min(), sf["ux"].min()), - ("mean u_x", ff["ux"].mean(), sf["ux"].mean()), - ("RMSE u_x (grille)", None, rmse["ux"]), - ("── ERREURS RELATIVES SUR GRILLE COMMUNE ──────────────────", None, None), - ("Erreur rel. médiane u_y (%)", rel_med["uy"], None), - ("Erreur rel. médiane u_x (%)", rel_med["ux"], None), - ("Erreur rel. moyenne u_y (%)", rel_moy["uy"], None), - ("Erreur rel. moyenne u_x (%)", rel_moy["ux"], None), - ("── INFORMATIONS MAILLAGE ──────────────────────────────────", None, None), - ("Nœuds FreeFEM", len(ff), None), - ("Nœuds SOFA", len(sf), None), - ("Rapport nœuds FF/SOFA", len(ff)/max(len(sf),1), None), -] - -for row in rows_table: - label, vff, vsf = row - if vff is None and vsf is None: - print(f"\n {label}") - continue - if vsf is None: - # Valeur unique (erreur relative, info maillage) - print(f" {label:<45} {vff:>12.4f}") - continue - if vff is None: - # RMSE (pas de valeur FF individuelle) - print(f" {label:<45} {'—':>12} {vsf:>12.4e}") - continue - delta = vsf - vff - dp = pct(vff, vsf) - print(f" {label:<45} {vff:>12.6f} {vsf:>12.6f} " - f"{delta:>12.2e} {dp:>10.2f}%") - -print("═" * 72) - -# ══════════════════════════════════════════════════════════════════════════════ -# 9. EXPORT CSV DU TABLEAU -# ══════════════════════════════════════════════════════════════════════════════ -summary_rows = [] -for f in FIELDS: - ff_max = np.abs(ff[f]).max() - sf_max = np.abs(sf[f]).max() - ff_min = ff[f].min() - sf_min = sf[f].min() - summary_rows.append({ - "champ": f, - "FF_max_abs": ff_max, - "SF_max_abs": sf_max, - "delta_max_abs": sf_max - ff_max, - "pct_max_abs": pct(ff_max, sf_max), - "FF_min": ff_min, - "SF_min": sf_min, - "delta_min": sf_min - ff_min, - "pct_min": pct(ff_min, sf_min), - "RMSE_grille": rmse[f], - "err_rel_med_pct": rel_med[f], - "err_rel_moy_pct": rel_moy[f], - "noeuds_FF": len(ff), - "noeuds_SF": len(sf), - }) - -pd.DataFrame(summary_rows).to_csv("comparaison_summary.csv", index=False) -print("\n ✔ comparaison_summary.csv") - -# ══════════════════════════════════════════════════════════════════════════════ -# 10. AFFICHAGE -# ══════════════════════════════════════════════════════════════════════════════ -import os, subprocess, sys - -output_dir = os.path.dirname(os.path.abspath(__file__)) -figures = ["fig1_maillages.png", "fig2_deplacements.png", - "fig3_profil_mediane.png", "fig4_erreur_relative.png", - "fig5_histogrammes.png"] - -print("\n Fichiers générés :") -for fname in figures + ["comparaison_summary.csv"]: - fpath = os.path.join(output_dir, fname) - size = os.path.getsize(fpath) // 1024 if os.path.exists(fpath) else 0 - print(f" • {fname} ({size} Ko)") - -# Ouvrir automatiquement les images avec le viewer Windows/Mac/Linux -print("\n Ouverture des figures...") -for fname in figures: - fpath = os.path.join(output_dir, fname) - if not os.path.exists(fpath): - continue - try: - if sys.platform.startswith("win"): - os.startfile(fpath) # Windows : ouvre avec la visionneuse par défaut - elif sys.platform == "darwin": - subprocess.Popen(["open", fpath]) - else: - subprocess.Popen(["xdg-open", fpath]) - except Exception as e: - print(f" Impossible d'ouvrir {fname} : {e}") - - -try: - plt.show() -except Exception: - pass - - - - - - -# ══════════════════════════════════════════════════════════════════════════════ -# 11. FIGURE 6 — VISUALISATION DES NŒUDS AVEC LEURS NUMÉROS -# ══════════════════════════════════════════════════════════════════════════════ - -def plot_nodes_with_numbers(ax, df, title, color_node='blue', show_all=False, max_nodes_to_show=100): - """ - Affiche les nœuds avec leurs numéros. - - Parameters: - ----------- - ax : axes matplotlib - df : DataFrame avec colonnes x, y - title : str - color_node : couleur des nœuds - show_all : bool - si True affiche tous les nœuds, sinon limite à max_nodes_to_show - max_nodes_to_show : int - nombre maximum de nœuds à étiqueter - """ - x = df["x"].values - y = df["y"].values - n_nodes = len(x) - - # Tracer tous les nœuds - ax.scatter(x, y, s=20, c=color_node, zorder=3, alpha=0.7, edgecolors='black', linewidth=0.5) - - # Décider quels nœuds étiqueter - if show_all or n_nodes <= max_nodes_to_show: - # Afficher tous les nœuds - indices = range(n_nodes) - else: - # Afficher seulement quelques nœuds (bords + quelques internes) - # Identifier les nœuds uniques par ligne y - y_unique = np.unique(np.round(y, 6)) - n_y = len(y_unique) - n_x = n_nodes // n_y - - indices_to_label = [] - - # Ajouter les coins - corners = [0, n_x-1, n_nodes - n_x, n_nodes - 1] - indices_to_label.extend(corners) - - # Ajouter quelques nœuds sur les bords - step = max(1, n_x // 10) - for i in range(0, n_x, step): - indices_to_label.append(i) # bord bas - indices_to_label.append(n_nodes - n_x + i) # bord haut - - # Ajouter quelques nœuds internes - step_y = max(1, n_y // 5) - step_x = max(1, n_x // 5) - for iy in range(0, n_y, step_y): - for ix in range(0, n_x, step_x): - idx = iy * n_x + ix - if idx not in indices_to_label: - indices_to_label.append(idx) - - indices = list(set(indices_to_label)) # Éliminer les doublons - - # Afficher les numéros - for idx in indices: - if idx < n_nodes: - ax.annotate(str(idx), - (x[idx], y[idx]), - xytext=(5, 5), - textcoords='offset points', - fontsize=8, - bbox=dict(boxstyle="round,pad=0.2", fc="white", alpha=0.8, ec="gray"), - alpha=0.8) - - ax.set_xlim(-0.2, 10.2) - ax.set_ylim(-0.1, 2.1) - ax.set_aspect("equal") - ax.set_xlabel("x (m)") - ax.set_ylabel("y (m)") - ax.set_title(f"{title}\n({n_nodes} nœuds", fontsize=11) - ax.grid(True, alpha=0.3) - - # Ajouter une légende sur le nombre de nœuds affichés - if n_nodes > max_nodes_to_show: - ax.text(0.02, 0.98, f"Affichage de {len(indices)}/{n_nodes} nœuds", - transform=ax.transAxes, fontsize=8, - verticalalignment='top', - bbox=dict(boxstyle="round,pad=0.2", fc="yellow", alpha=0.7)) - - -# Créer la figure 6 -fig6, axes6 = plt.subplots(2, 2, figsize=(16, 12)) -fig6.suptitle("Visualisation des nœuds avec leurs numéros", - fontsize=14, fontweight="bold") - -# Sous-figure pour FreeFEM (tous les nœuds) -plot_nodes_with_numbers(axes6[0, 0], ff, "FreeFEM - Tous les nœuds", - color_node='#1a6faf', show_all=True) - -# Sous-figure pour FreeFEM (limité pour lisibilité) -plot_nodes_with_numbers(axes6[0, 1], ff, "FreeFEM - Échantillon de nœuds", - color_node='#1a6faf', show_all=False, max_nodes_to_show=100) - -# Sous-figure pour SOFA (tous les nœuds) -plot_nodes_with_numbers(axes6[1, 0], sf, "SOFA - Tous les nœuds", - color_node='#d62728', show_all=True) - -# Sous-figure pour SOFA (avec détails supplémentaires) -ax = axes6[1, 1] -x_sf = sf["x"].values -y_sf = sf["y"].values -ax.scatter(x_sf, y_sf, s=50, c='#d62728', zorder=3, alpha=0.8, edgecolors='black', linewidth=0.5) - -# Afficher tous les numéros pour SOFA (peu nombreux) -for idx, (xi, yi) in enumerate(zip(x_sf, y_sf)): - ax.annotate(str(idx), - (xi, yi), - xytext=(8, 8), - textcoords='offset points', - fontsize=10, - fontweight='bold', - bbox=dict(boxstyle="round,pad=0.3", fc="white", alpha=0.9, ec="red", lw=1.5), - alpha=0.9) - -# Ajouter les coordonnées pour quelques nœuds clés -for idx, (xi, yi) in enumerate(zip(x_sf, y_sf)): - if idx % 5 == 0 or idx < 4 or idx >= len(x_sf)-4: # Afficher les coordonnées pour certains nœuds - ax.annotate(f"({xi:.1f}, {yi:.1f})", - (xi, yi), - xytext=(8, -15), - textcoords='offset points', - fontsize=7, - alpha=0.7) - -ax.set_xlim(-0.2, 10.2) -ax.set_ylim(-0.1, 2.1) -ax.set_aspect("equal") -ax.set_xlabel("x (m)") -ax.set_ylabel("y (m)") -ax.set_title(f"SOFA - Détail des nœuds\n({len(sf)} nœuds, grille régulière)", fontsize=11) -ax.grid(True, alpha=0.3) - -plt.tight_layout() -plt.savefig("fig6_noeuds_avec_numeros.png", dpi=150, bbox_inches="tight") -print(" ✔ fig6_noeuds_avec_numeros.png") - - -# ══════════════════════════════════════════════════════════════════════════════ -# 12. FIGURE 7 — TOPOLOGIE DU MAILLAGE SOFA (grille structurée) -# ══════════════════════════════════════════════════════════════════════════════ - -def plot_sofa_grid_structure(ax, df, nx, ny): - """Visualise la structure de grille de SOFA avec les connexions""" - x = df["x"].values - y = df["y"].values - - # Déterminer nx, ny à partir des données si non fournis - if nx is None or ny is None: - y_unique = np.unique(np.round(y, 6)) - ny = len(y_unique) - nx = len(x) // ny - - # Tracer les lignes de la grille - for i in range(ny): - row_indices = slice(i*nx, (i+1)*nx) - ax.plot(x[row_indices], y[row_indices], 'b-', linewidth=1, alpha=0.5, color='gray') - - for j in range(nx): - col_indices = list(range(j, nx*ny, nx)) - ax.plot(x[col_indices], y[col_indices], 'b-', linewidth=1, alpha=0.5, color='gray') - - # Tracer les nœuds - ax.scatter(x, y, s=30, c='#d62728', zorder=3, alpha=0.8, edgecolors='black') - - # Numéroter les nœuds - for idx, (xi, yi) in enumerate(zip(x, y)): - ax.annotate(str(idx), - (xi, yi), - xytext=(5, 5), - textcoords='offset points', - fontsize=9, - bbox=dict(boxstyle="round,pad=0.2", fc="white", alpha=0.8)) - - ax.set_xlim(-0.2, 10.2) - ax.set_ylim(-0.1, 2.1) - ax.set_aspect("equal") - ax.set_xlabel("x (m)") - ax.set_ylabel("y (m)") - ax.set_title(f"Structure de la grille SOFA\n{nx} × {ny} = {nx*ny} nœuds", fontsize=11) - ax.grid(True, alpha=0.3) - -# Créer la figure 7 -fig7, ax7 = plt.subplots(1, 1, figsize=(12, 6)) - -# Déterminer nx, ny à partir des données SOFA -y_unique_sf = np.unique(np.round(sf["y"].values, 6)) -ny_sf = len(y_unique_sf) -nx_sf = len(sf) // ny_sf - -plot_sofa_grid_structure(ax7, sf, nx_sf, ny_sf) - -plt.tight_layout() -plt.savefig("fig7_grille_sofa.png", dpi=150, bbox_inches="tight") -print(" ✔ fig7_grille_sofa.png") - - -# ══════════════════════════════════════════════════════════════════════════════ -# 13. FIGURE 8 — CARTE DE CHALEUR DE LA DENSITÉ DE NŒUDS (FreeFEM) -# ══════════════════════════════════════════════════════════════════════════════ - -def plot_node_density_heatmap(ax, df, title, bins=50): - """Affiche une carte de chaleur de la densité de nœuds""" - x = df["x"].values - y = df["y"].values - - # Créer l'histogramme 2D - H, xedges, yedges = np.histogram2d(x, y, bins=bins, range=[[0, 10], [0, 2]]) - - # Afficher la carte de chaleur - im = ax.imshow(H.T, origin='lower', extent=[0, 10, 0, 2], - cmap='hot', aspect='auto', alpha=0.8) - - # Ajouter la colorbar - plt.colorbar(im, ax=ax, label='Nombre de nœuds par cellule') - - ax.set_xlabel("x (m)") - ax.set_ylabel("y (m)") - ax.set_title(title) - ax.grid(True, alpha=0.2) - -fig8, axes8 = plt.subplots(1, 2, figsize=(14, 5)) -fig8.suptitle("Densité de nœuds dans le maillage", fontsize=13, fontweight="bold") - -plot_node_density_heatmap(axes8[0], ff, "FreeFEM - Distribution des nœuds", bins=40) -plot_node_density_heatmap(axes8[1], sf, "SOFA - Distribution des nœuds (grille régulière)", bins=20) - -plt.tight_layout() -plt.savefig("fig8_densite_noeuds.png", dpi=150, bbox_inches="tight") -print(" ✔ fig8_densite_noeuds.png") - - -# ══════════════════════════════════════════════════════════════════════════════ -# Mettre à jour la liste des figures -figures += ["fig6_noeuds_avec_numeros.png", "fig7_grille_sofa.png", "fig8_densite_noeuds.png"] - - - - - - -# prise de note - -""" -POURQUOI ON PEUT COMPARER LES DEUX CODES MALGRÉ DES MAILLAGES DIFFÉRENTS ? -=========================================================================== - -PROBLÈME DE DÉPART : - • FreeFEM génère un maillage triangulaire NON-STRUCTURÉ avec ~4800 nœuds. - Les nœuds sont placés de façon irrégulière, denses là où le gradient est - fort (encastrement gauche), plus rares au centre. - - • SOFA utilise une RegularGridTopology 11×4 = 44 nœuds seulement, - placés sur une grille parfaitement régulière. - Les lignes y sont : 0.0, 0.67, 1.33, 2.0 - → aucun nœud à y = 1.0 exactement ! - -SOLUTION — LA GRILLE COMMUNE : - 1. On crée une grille uniforme 300×80 = 24 000 points couvrant [0,10]×[0,2]. - 2. On interpole chaque champ (ux, uy) sur cette grille depuis les nœuds bruts : - - FreeFEM : scipy.griddata(..., method='linear') - → triangulation Delaunay + interpolation linéaire par morceaux (P1) - → très fidèle car maillage dense - - SOFA : scipy.griddata(..., method='cubic') - → interpolation cubique car 44 points sont trop peu pour 'linear' - → lisse le champ entre les nœuds réguliers - - 3. Une fois les deux champs sur la même grille, on peut : - diff[i,j] = SOFA_interpolé[i,j] - FreeFEM_interpolé[i,j] - et calculer RMSE, erreur relative, histogrammes, etc. - -POURQUOI LA COMPARAISON A DU SENS ? - Les deux codes résolvent la MÊME équation (élasticité linéaire plane, - même E, nu, gravité, encastrement gauche). La différence vient de : - a) La formulation : plane_strain (Vec3, λ=Eν/((1+ν)(1-2ν))) vs - les coefficients Lamé 3D de FreeFEM (identiques à plane_strain) - b) La discrétisation : triangles P1 vs Q1 subdivisé en triangles - c) Le nombre de DDL : 44 nœuds SOFA << 4814 nœuds FreeFEM - → SOFA est moins précis, mais donne la bonne tendance physique. - - Le tableau récapitulatif vous donne l'écart % entre les deux pour - chaque quantité d'intérêt (flèche maximale, RMSE, erreur médiane...). -""" \ No newline at end of file diff --git a/examples/Freefem/3d-cube.edp b/examples/Freefem/3d-cube.edp deleted file mode 100644 index 74916a82..00000000 --- a/examples/Freefem/3d-cube.edp +++ /dev/null @@ -1,78 +0,0 @@ -load "msh3" -load "medit" -load "iovtk" - - -// Pb elasticite dans un cube -mesh3 Th = cube(10, 10, 10); - -// Vérifier les labels du cube -cout << "Labels bords : " << endl; -int[int] labs = labels(Th); -for(int i = 0; i < labs.n; i++) - cout << " label " << labs[i] << endl; - -real E = 210e9; -real sigma = 0.3; -real rho = 7800.; -real gravity = -9.81 * rho; -real mu = E / (2.*(1.+sigma)); -real lambda = E*sigma / ((1.+sigma)*(1.-2.*sigma)); - -// Zmin/zmax -real zmin = 1e100, zmax = -1e100; -for (int i = 0; i < Th.nv; i++) { - zmin = min(zmin, Th(i).z); - zmax = max(zmax, Th(i).z); -} -cout << "zmin=" << zmin << " zmax=" << zmax << endl; - -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)(gravity * vz) - + on(1, ux=0, uy=0, uz=0); // encastrement - -Elasticite3D; - -cout << "Deplacement max uz = " << uz[].linfty << " m" << endl; -cout << "Theorique uz_max = " << rho*9.81*(zmax-zmin)^2/(2.*E) << " m" << endl; - -fespace Wh(Th, P1); -Wh normu = sqrt(ux^2 + uy^2 + uz^2); -medit("Norme", Th, normu, wait=1); -medit("uz", Th, uz, wait=1); - -real umax = max(ux[].linfty, max(uy[].linfty, uz[].linfty)); -real coef = 0.1*(zmax-zmin)/umax; -cout << "Coefficient amplification = " << coef << endl; -mesh3 Thd = movemesh3(Th, transfo=[x+coef*ux, y+coef*uy, z+coef*uz]); -//savemesh(Thd, "cube_deforme.mesh"); -medit("Cube deforme", Thd, wait=1); - - - - - - - -// Exporter la solution en format VTK -int[int] order = [1, 1, 1]; // Ordre P1 -savevtk("cube_deplacement.vtu", Th, [ux, uy, uz], dataname="Deplacement", order=order); -savevtk("cube_norme.vtu", Th, normu, dataname="Norme", order=order); -savevtk("cube_uz.vtu", Th, uz, dataname="Uz", order=order); - -// Exporter le maillage déformé -mesh3 Thd = movemesh3(Th, transfo=[x + 1e6*ux, y + 1e6*uy, z + 1e6*uz]); -savevtk("cube_deforme.vtu", Thd, [ux, uy, uz], dataname="Deplacement", order=order); - -cout << "Fichiers VTK exportés. Ouvrez-les avec ParaView." << endl; \ No newline at end of file diff --git a/examples/Freefem/barre-1d-solution-man.edp b/examples/Freefem/barre-1d-solution-man.edp deleted file mode 100644 index 406c39f2..00000000 --- a/examples/Freefem/barre-1d-solution-man.edp +++ /dev/null @@ -1,44 +0,0 @@ -// Solution manufacturée normalisée -// On résout : -u'' = pi^2*sin(pi*x) EA simplifie =1 -real L = 1.0; -real E = 1.0; -real A = 1.0; -real pi = 3.14159265358979; - -int[int] nvals = [2, 4, 8, 16, 32, 64]; - -ofstream file("convergence_grossier.txt"); -file << "# n\th\terrL2\terrH1" << endl; - -for (int idx = 0; idx < nvals.n; idx++) { - int n = nvals[idx]; - real h = L / n; - - meshL Th = segment(n, [x * L]); - fespace Vh(Th, P1); - Vh u, v; - - problem Traction(u, v) = - int1d(Th)( dx(u) * dx(v) ) - - int1d(Th)( pi * pi * sin(pi * x) * v ) - + on(1, u = 0) - + on(2, u = 0); - - Traction; - - Vh uexact = sin(pi * x); - - // Erreur sur les noeuds uniquement -real errL2 = 0, errH1 = 0; - -errL2 = sqrt(int1d(Th)( (u - sin(pi*x))^2 )); -errH1 = sqrt(int1d(Th)( (dx(u) - pi*cos(pi*x))^2 )); - - - cout << "n = " << n - << ", h = " << h - << ", errL2 = " << errL2 - << ", errH1 = " << errH1 << endl; - - file << n << "\t" << h << "\t" << errL2 << "\t" << errH1 << endl; -} \ No newline at end of file diff --git a/examples/Freefem/barre_exacte.edp b/examples/Freefem/barre_exacte.edp deleted file mode 100644 index 55adbef6..00000000 --- a/examples/Freefem/barre_exacte.edp +++ /dev/null @@ -1,42 +0,0 @@ -// === Élasticité linéaire - BARRE en traction === -// Version pour étude de convergence - -real L = 1.0; -real E = 1e6; -real A = 0.001; -real F = 1000; - -// Différentes valeurs de n pour l'étude de convergence - -int[int] nvals = [2, 4, 8, 16, 32, 64]; -ofstream file("convergence_data.txt"); -file << "# n\th\terrL2\terrH1" << endl; - -for (int idx = 0; idx < nvals.n; idx++) { - int n = nvals[idx]; - real h = L/n; - - meshL Th = segment(n, [x*L]); - fespace Vh(Th, P1); - Vh u, v; - - problem Traction(u, v) = - int1d(Th)( E * A * dx(u) * dx(v) ) - - int0d(Th, 2)( F * v ) - + on(1, u=0); - - Traction; - - Vh uexact = (F/(E*A)) * x; - - // Erreur L2 - real errL2 = sqrt(int1d(Th)((u - uexact)^2)); - - // Erreur H1 (semi-norme) - real errH1 = sqrt(int1d(Th)((dx(u) - dx(uexact))^2)); - - cout << "n = " << n << ", h = " << h << ", errL2 = " << errL2 << ", errH1 = " << errH1 << endl; - - // Sauvegarde - file << n << "\t " << h << "\t " << errL2 << "\t " << errH1 << endl; -} \ No newline at end of file diff --git a/examples/Freefem/beam-2d.edp b/examples/Freefem/beam-2d.edp deleted file mode 100644 index 2f9269ef..00000000 --- a/examples/Freefem/beam-2d.edp +++ /dev/null @@ -1,104 +0,0 @@ -// Beam under its own weight - -// Parameters -real E = 21.5; -real sigma = 0.29; -real gravity = -0.05; -int n=140; -int m=35; - - -// Maillage non structure " avec courbes parametriques " -border a(t=2, 0){x=0; y=t ;label=1;} // gauche cote de la poutre à encastre -border b(t=0, 10){x=t; y=0 ;label=2;} // bas (libre) -border c(t=0, 2){x=10; y=t ;label=3;} // droite (libre) -border d(t=0, 10){x=10-t; y=2; label=4;} // haut (libre) -mesh th = buildmesh(b(n) + c(m) + d(n) + a(m)); // un maillage equilibre previlige celui là -// n=100 unite de long pour les cotes haut et bas, et n=25 unite de haut pour les cote gauche et droite -//mesh th = buildmesh(b(n) + c(n) + d(n) + a(n)); // evite celui là le maillage est raffine uniquement aux extrimite a et c aka gauche et droite -// Fespace element finis espace P1 ; -plot (th,cmm="maillage-poutre", wait=1); -fespace Vh(th, [P1, P1]); -Vh [uu, vv], [w, s]; - -// Macro -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 - -// Coefficients de Lamé -real mu = E/(2*(1 + sigma)); -real lambda = E*sigma/((1 + sigma)*(1 - 2*sigma)); - -// Problème -solve Elasticity ([uu, vv], [w, s], solver=sparsesolver) - = int2d(th)( - lambda*div(w, s)*div(uu, vv) - + 2.*mu*(epsilon(w, s)' * epsilon(uu, vv)) - ) - - int2d(th)(gravity*s) - + on(1, uu=0, vv=0); - -// Résultats -cout << "Max displacement = " << uu[].linfty << ", " << vv[].linfty << endl; - -// Visualisation -plot([uu, vv], wait=1, cmm="Deplacement"); - - -mesh th1 = movemesh(th, [x + uu, y + vv]); -plot( th1, wait=1, cmm="Maillage initial et deformation"); - -// Espace pour les contraintes -fespace Wh(th, P1); -Wh sigmaxx, sigmayy, sigmaxy, sigmavm; - -// Calcul des contraintes -sigmaxx = lambda * (dx(uu) + dy(vv)) + 2*mu*dx(uu); -sigmayy = lambda * (dx(uu) + dy(vv)) + 2*mu*dy(vv); -sigmaxy = mu*(dy(uu) + dx(vv)); - -// Contrainte de von Mises pour plasticité -sigmavm = sqrt(sigmaxx^2 + sigmayy^2 - sigmaxx*sigmayy + 3*sigmaxy^2); - -// Visualisation -plot(sigmaxx, wait=1, fill=1, value=1, cmm="Contrainte sigma_xx (MPa)"); -plot(sigmavm, wait=1, fill=1, value=1, cmm="Contrainte von Mises (MPa)"); - -// Valeurs extrêmes -cout << "sigma_xx min = " << sigmaxx[].min << ", max = " << sigmaxx[].max << endl; -cout << "sigma_vm max = " << sigmavm[].max << endl; - - - -// Calcul des réactions par intégration sur le bord encastré (label=1) -real reactionx = int1d(th, 1)( (sigmaxx * N.x + sigmaxy * N.y) ); -real reactiony = int1d(th, 1)( (sigmaxy * N.x + sigmayy * N.y) ); - -// Poids total -real poidstotal = abs(gravity) * int2d(th)(1); - -// Affichage des résultats - -cout << "Reaction horizontale a l'encastrement : " << reactionx << endl; -cout << "Reaction verticale a l'encastrement : " << reactiony << endl; -cout << "Poids total de la poutre : " << poidstotal << endl; -cout << "Rapport reaction verticale / poids : " << reactiony / poidstotal << endl; - -// Test d'équilibre -real erreurequilibre = abs(reactiony/poidstotal - 1) * 100; -cout << "Erreur d'equilibre : " << erreurequilibre << " %" << endl; - -if (erreurequilibre < 5) - cout << "equilibre verifie " << endl; -else - cout << " Attention : desequilibre detecte !" << endl; - -// Vérification du moment à l'encastrement -real moment = int1d(th, 1)( (y - 1) * (sigmaxy * N.x + sigmayy * N.y) - (x) * (sigmaxx * N.x + sigmaxy * N.y) ); -cout << "Moment a l'encastrement : " << moment << endl; - - - - - diff --git a/examples/Freefem/elasticity_3D.edp b/examples/Freefem/elasticity_3D.edp deleted file mode 100644 index 902cc751..00000000 --- a/examples/Freefem/elasticity_3D.edp +++ /dev/null @@ -1,105 +0,0 @@ -load "msh3" -load "medit" -load "iovtk" - // remarque : la comparaison theorique qu'on a obtenue n'est pas vraiment valide - // car on compare une geometrie tres difficle à une poutre, mais on reste tjrs dans le cadre du respect de la loi d la physique -mesh3 Th("C:/Program Files (x86)/FreeFem++/examples/3d/dodecaedre01.mesh"); - -// Calculer zmin -real zmin = 1e100, zmax = -1e100; -for (int i = 0; i < Th.nv; i++) { - zmin = min(zmin, Th(i).z); - zmax = max(zmax, Th(i).z); -} -cout << "zmin=" << zmin << " zmax=" << zmax << endl; - -// Changer le label des faces du bas (proches de zmin) -Th = change(Th, flabel=(abs(z - zmin) < 0.1) ? 1 : 0); - -// Vérifier les labels modifiés -cout << "Nouveaux labels :" << endl; -for (int i = 0; i < Th.nbe; i++) { - if (Th.be(i).label == 1) - cout << "Face " << i << " label = " << Th.be(i).label << endl; -} - -// Vérifier tous les labels présents -int[int] labs = labels(Th); -cout << "Labels presents dans le maillage :" << endl; -for(int i = 0; i < labs.n; i++) - cout << " label " << labs[i] << endl; - -// Paramètres physiques -real E = 210e9; -real sigma = 0.3; -real rho = 7800.; -real gravity = -9.81 * rho; -real mu = E / (2.*(1.+sigma)); -real lambda = E*sigma / ((1.+sigma)*(1.-2.*sigma)); - -// Espace éléments finis -fespace Vh(Th, [P1,P1,P1]); -Vh [ux,uy,uz],[vx,vy,vz]; - -// Macros (attention aux parenthèses superflues) -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 - -// Problème d'élasticité on fixe les faces label=1 -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)(gravity * vz) - + on(1, ux=0, uy=0, uz=0); - -// Résolution -Elasticite3D; - -// Résultats -cout << "Deplacement max |ux| = " << ux[].linfty << " m" << endl; -cout << "Deplacement max |uy| = " << uy[].linfty << " m" << endl; -cout << "Deplacement max |uz| = " << uz[].linfty << " m" << endl; - -// Création d'un espace scalaire P1 -fespace Wh(Th, P1); -Wh normu; - - -normu = sqrt(ux*ux + uy*uy + uz*uz); - -cout << "Deplacement max total = " << normu[].linfty << " m" << endl; - - - -cout << "Theorique (poutre) = " << rho*9.81*(zmax-zmin)^2/(2.*E) << " m" << endl; - -// Visualisation -medit("Norme du deplacement", Th, normu, wait=1); -medit("Deplacement vertical uz", Th, uz, wait=1); - -// Maillage déformé -real umax = max(ux[].linfty, max(uy[].linfty, uz[].linfty)); -real coef = 0.1*(zmax-zmin)/umax; -mesh3 Thd = movemesh3(Th, transfo=[x+coef*ux, y+coef*uy, z+coef*uz]); -medit("Dodecaedre deforme", Thd, wait=1); - - - - - - -// Exporter la solution en format VTK -//int[int] order = [1, 1, 1]; -//savevtk("cube_deplacement.vtu", Th, [ux, uy, uz], dataname="Deplacement", order=order); -//savevtk("cube_norme.vtu", Th, normu, dataname="Norme", order=order); -//savevtk("cube_uz.vtu", Th, uz, dataname="Uz", order=order); - -// Exporter le maillage déformé -//mesh3 Thd = movemesh3(Th, transfo=[x + 1e6*ux, y + 1e6*uy, z + 1e6*uz]); -//savevtk("cube_deforme.vtu", Thd, [ux, uy, uz], dataname="Deplacement", order=order); - diff --git a/barre_1d_sofa.py b/examples/Freefem/sofa/barre-traction.py similarity index 100% rename from barre_1d_sofa.py rename to examples/Freefem/sofa/barre-traction.py diff --git a/beam_2d_sofa.py b/examples/Freefem/sofa/beam-2d.py similarity index 100% rename from beam_2d_sofa.py rename to examples/Freefem/sofa/beam-2d.py From 91b4b401c9a3114dea96a79628d31b4896dc3bec Mon Sep 17 00:00:00 2001 From: fatima Date: Mon, 30 Mar 2026 14:06:08 +0200 Subject: [PATCH 4/5] Add correct barre-traction.py for 1D traction simulation --- examples/Freefem/sofa/barre-traction.py | 148 ++++++++++++++++++++++++ 1 file changed, 148 insertions(+) create mode 100644 examples/Freefem/sofa/barre-traction.py diff --git a/examples/Freefem/sofa/barre-traction.py b/examples/Freefem/sofa/barre-traction.py new file mode 100644 index 00000000..c4af9891 --- /dev/null +++ b/examples/Freefem/sofa/barre-traction.py @@ -0,0 +1,148 @@ +import Sofa +import Sofa.Core +import Sofa.Simulation +import SofaRuntime +import numpy as np + +# Exporter les resultats +class ExportController(Sofa.Core.Controller): + """Récupère les déplacements après convergence et les écrit dans un .txt""" + + def __init__(self, dofs_node, output_file="sofa_deplacement.txt", *args, **kwargs): + super().__init__(*args, **kwargs) + self.dofs_node = dofs_node + self.output_file = output_file + self.exported = False + + def onAnimateEndEvent(self, event): + """Appelé automatiquement par SOFA à la fin de chaque pas de temps""" + if self.exported: + return + + + pos = self.dofs_node.position.array() + + x_deform = pos[:, 0] + u_x = x_deform - self.x_initial + + + order = np.argsort(self.x_initial) + x0_sorted = self.x_initial[order] + ux_sorted = u_x[order] + + with open(self.output_file, 'w') as f: + f.write("x_initial u_x\n") + for xi, ui in zip(x0_sorted, ux_sorted): + f.write(f"{xi:.6f} {ui:.6f}\n") + + print(f"[ExportController] Déplacements exportés → {self.output_file}") + print(f" u_x(L) = {ux_sorted[-1]:.6f} (attendu : 1.0)") + self.exported = True + + def onSimulationInitDoneEvent(self, event): + """Stocke les positions initiales au démarrage""" + pos = self.dofs_node.position.array() + self.x_initial = pos[:, 0].copy() + print(f"[ExportController] {len(self.x_initial)} noeuds initialisés.") + + +# creation de la scene +def createScene(rootNode): + L = 1.0 + E_eff = 1000.0 + F = 1000.0 + N = 9 + + rootNode.addObject('DefaultAnimationLoop') + rootNode.addObject('VisualStyle', + displayFlags="showBehaviorModels showForceFields showWireframe") + + plugins = rootNode.addChild('plugins') + for plugin in [ + "Elasticity", + "Sofa.Component.Constraint.Projective", + "Sofa.Component.Engine.Select", + "Sofa.Component.LinearSolver.Direct", + "Sofa.Component.MechanicalLoad", + "Sofa.Component.ODESolver.Backward", + "Sofa.Component.StateContainer", + "Sofa.Component.Topology.Container.Dynamic", + "Sofa.Component.Topology.Container.Grid", + "Sofa.Component.Visual", + "Sofa.GL.Component.Rendering3D", + ]: + plugins.addObject('RequiredPlugin', name=plugin) + + bar = rootNode.addChild('bar') + + bar.addObject('NewtonRaphsonSolver', + name="newtonSolver", + printLog=True, + warnWhenLineSearchFails=True, + maxNbIterationsNewton=1, + maxNbIterationsLineSearch=1, + lineSearchCoefficient=1, + relativeSuccessiveStoppingThreshold=0, + absoluteResidualStoppingThreshold=1e-7, + absoluteEstimateDifferenceThreshold=1e-12, + relativeInitialStoppingThreshold=1e-12, + relativeEstimateDifferenceThreshold=0) + + bar.addObject('SparseLDLSolver', + name="linearSolver", + template="CompressedRowSparseMatrixd") + + bar.addObject('StaticSolver', + name="staticSolver", + newtonSolver="@newtonSolver", + linearSolver="@linearSolver") + + bar.addObject('RegularGridTopology', + name="grid", + min="0 0 0", + max=f"{L} 0 0", + n=f"{N} 1 1") + + dofs = bar.addObject('MechanicalObject', + template="Vec3d", + name="dofs", + showObject=True, + showObjectScale=0.02) + + edges = bar.addChild('edges') + edges.addObject('EdgeSetTopologyContainer', + name="topology", + position="@../dofs.position", + edges="@../grid.edges") + edges.addObject('LinearSmallStrainFEMForceField', + name="FEM", + youngModulus=E_eff, + poissonRatio=0.0, + topology="@topology") + + bar.addObject('BoxROI', + name="fixed_roi", + box="-0.01 -1 -1 0.01 1 1", + drawBoxes=False) + bar.addObject('FixedProjectiveConstraint', + indices="@fixed_roi.indices") + + bar.addObject('BoxROI', + name="load_roi", + box=f"{L-0.01} -1 -1 {L+0.01} 1 1", + drawBoxes=False) + bar.addObject('ConstantForceField', + indices="@load_roi.indices", + forces=f"{F} 0 0", + showArrowSize=1e-4) + + + rootNode.addObject( + ExportController( + dofs_node = dofs, + output_file = "sofa_deplacement.txt", + name = "exportCtrl" + ) + ) + + return rootNode \ No newline at end of file From b234087534747d3c16e754d123861bc51544b0b0 Mon Sep 17 00:00:00 2001 From: Themis Skamagkis Date: Mon, 30 Mar 2026 16:26:43 +0200 Subject: [PATCH 5/5] apply corrections --- examples/Freefem/sofa/barre-traction.py | 17 ++++++++------ examples/Freefem/sofa/beam-2d.py | 31 ++++++++----------------- 2 files changed, 20 insertions(+), 28 deletions(-) diff --git a/examples/Freefem/sofa/barre-traction.py b/examples/Freefem/sofa/barre-traction.py index c4af9891..9b609899 100644 --- a/examples/Freefem/sofa/barre-traction.py +++ b/examples/Freefem/sofa/barre-traction.py @@ -1,3 +1,6 @@ +""" +1D Bar Simulation - Traction Load +""" import Sofa import Sofa.Core import Sofa.Simulation @@ -5,10 +8,10 @@ import numpy as np # Exporter les resultats -class ExportController(Sofa.Core.Controller): +class DisplacementExporter(Sofa.Core.Controller): """Récupère les déplacements après convergence et les écrit dans un .txt""" - def __init__(self, dofs_node, output_file="sofa_deplacement.txt", *args, **kwargs): + def __init__(self, dofs_node, output_file="sofa_displacement.txt", *args, **kwargs): super().__init__(*args, **kwargs) self.dofs_node = dofs_node self.output_file = output_file @@ -35,7 +38,7 @@ def onAnimateEndEvent(self, event): for xi, ui in zip(x0_sorted, ux_sorted): f.write(f"{xi:.6f} {ui:.6f}\n") - print(f"[ExportController] Déplacements exportés → {self.output_file}") + print(f"[DisplacementExporter] Déplacements exportés → {self.output_file}") print(f" u_x(L) = {ux_sorted[-1]:.6f} (attendu : 1.0)") self.exported = True @@ -43,7 +46,7 @@ def onSimulationInitDoneEvent(self, event): """Stocke les positions initiales au démarrage""" pos = self.dofs_node.position.array() self.x_initial = pos[:, 0].copy() - print(f"[ExportController] {len(self.x_initial)} noeuds initialisés.") + print(f"[DisplacementExporter] {len(self.x_initial)} noeuds initialisés.") # creation de la scene @@ -138,11 +141,11 @@ def createScene(rootNode): rootNode.addObject( - ExportController( + DisplacementExporter( dofs_node = dofs, - output_file = "sofa_deplacement.txt", + output_file = "sofa_displacement.txt", name = "exportCtrl" ) ) - return rootNode \ No newline at end of file + return rootNode diff --git a/examples/Freefem/sofa/beam-2d.py b/examples/Freefem/sofa/beam-2d.py index c5b8b3f3..85b48287 100644 --- a/examples/Freefem/sofa/beam-2d.py +++ b/examples/Freefem/sofa/beam-2d.py @@ -1,13 +1,5 @@ """ -beam_2d_sofa.py -=============== -Equivalent Python de beam-2d.scn (SOFA) -Poutre 2D sous gravité — formulation déformations planes (Vec3d) -Equivalent exact de beam-2d.edp (FreeFEM) - - -Export : - Génère sofa_results.txt avec x, y, ux, uy pour chaque nœud +2D Beam Simulation - Plane Deformation Under Gravity """ import Sofa @@ -24,14 +16,13 @@ GRAVITY = -0.05 NX, NY = 141, 36 TOTAL_MASS = 20.0 -OUTPUT_FILE = "sofa_results.txt" +OUTPUT_FILE = "sofa_displacement.txt" -class ExportController(Sofa.Core.Controller): +class DisplacementExporter(Sofa.Core.Controller): """ - Contrôleur SOFA : exporte les positions déformées après convergence. - Répond à onAnimateEndEvent (fin de chaque pas de temps). + Recovers nodal displacements after convergence and writes them in a .txt file """ def __init__(self, *args, **kwargs): @@ -46,20 +37,18 @@ def onAnimateEndEvent(self, event): if self.exported: return - # Laisser le temps à la simulation de converger - # 50 itérations suffisent généralement if self.iteration < 50: return root = self.getContext().getRootContext() strain = root.getChild("beam_plane_strain") if strain is None: - print("[ExportController] Nœud beam_plane_strain introuvable") + print("[DisplacementExporter] Nœud beam_plane_strain introuvable") return dofs = strain.getObject("dofs") if dofs is None: - print("[ExportController] MechanicalObject 'dofs' introuvable") + print("[DisplacementExporter] MechanicalObject 'dofs' introuvable") return # positions déformées @@ -88,12 +77,12 @@ def onAnimateEndEvent(self, event): uy = disp[i, 1] f.write(f"{x:.6f} {y:.6f} {ux:.8f} {uy:.8f}\n") - print(f"[ExportController] ✓ Exporté {len(pos_def)} nœuds → {OUTPUT_FILE}") + print(f"[DisplacementExporter] ✓ Exporté {len(pos_def)} nœuds → {OUTPUT_FILE}") self.exported = True root.animate.value = False except Exception as e: - print(f"[ExportController] Erreur d'écriture: {e}") + print(f"[DisplacementExporter] Erreur d'écriture: {e}") def _build_initial_grid(nx, ny, y_offset=0.0): @@ -132,7 +121,7 @@ def createScene(rootNode): displayFlags="showBehaviorModels showForceFields showWireframe") # ── Contrôleur d'export - rootNode.addObject(ExportController(name="exporter", node=rootNode)) + rootNode.addObject(DisplacementExporter(name="exporter", node=rootNode)) # ── Nœud beam_plane_strain (Vec3d) strain = rootNode.addChild("beam_plane_strain") @@ -210,4 +199,4 @@ def createScene(rootNode): if __name__ == "__main__": # Pour exécution directe (non-SOFA) print(f"Paramètres: NX={NX}, NY={NY}, output={OUTPUT_FILE}") - print(f"Total nœuds: {NX*NY}") \ No newline at end of file + print(f"Total nœuds: {NX*NY}")