freefem-code-linear elasticity

examples/Freefem/edp/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;

examples/Freefem/edp/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;
+}

examples/Freefem/edp/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;
+}

examples/Freefem/edp/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;
+
+
+
+
+

examples/Freefem/edp/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);
+