Convert XML scenes to Python scripts and integrate FreeFem simulations with SOFA comparison

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)
+

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)
+
+

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<th.nv; i++){{
+        real xi = th(i).x;
+        real yi = th(i).y;
+        fout << xi << " "
+             << yi << " "
+             << uu(xi, yi) << " "
+             << vv(xi, yi) << " "
+             << sigmaxx(xi, yi) << " "
+             << sigmayy(xi, yi) << " "
+             << sigmaxy(xi, yi) << " "
+             << sigmavm(xi, yi) << endl;
+    }}
+}}
+""".format(E=E, NU=NU, GRAVITY=GRAVITY,
+           N_LONG=N_LONG, N_HAUT=N_HAUT,
+           OUTPUT_FILE=OUTPUT_FILE)
+
+# ========== EXÉCUTION ==========
+print("Lancement FreeFEM...")
+with tempfile.NamedTemporaryFile(mode='w', suffix='.edp', delete=False) 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.stderr:
+        print("Erreurs:")
+        print(result.stderr)
+finally:
+    os.unlink(temp_file)
+
+# ========== LECTURE DES RÉSULTATS ==========
+if os.path.exists(OUTPUT_FILE):
+
+    df = pd.read_csv(OUTPUT_FILE, sep='\s+')
+    print("\n=== RÉSULTATS ===")
+    print(f"Max uy = {df['uy'].max()}")
+    print(f"Max ux = {df['ux'].max()}")
+    print(f"Max sigmaxx = {df['sigmaxx'].max()}")
+    print(f"Max sigmavm = {df['sigmavm'].max()}")
+
+
+else:
+    print(f"Fichier {OUTPUT_FILE} non trouvé")
+
+
+