Skip to content

Commit 532ab27

Browse files
Fimacheth-skam
Fimache
authored and
th-skam
committed
Add FreeFem++ elasticity simulation files
1 parent 0d99ead commit 532ab27

File tree

5 files changed

+373
-0
lines changed

5 files changed

+373
-0
lines changed

‎examples/Freefem/3d-cube.edp‎

Lines changed: 78 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,78 @@
1+
load "msh3"
2+
load "medit"
3+
load "iovtk"
4+
5+
6+
// Pb elasticite dans un cube
7+
mesh3 Th = cube(10, 10, 10);
8+
9+
// Vérifier les labels du cube
10+
cout << "Labels bords : " << endl;
11+
int[int] labs = labels(Th);
12+
for(int i = 0; i < labs.n; i++)
13+
cout << " label " << labs[i] << endl;
14+
15+
real E = 210e9;
16+
real sigma = 0.3;
17+
real rho = 7800.;
18+
real gravity = -9.81 * rho;
19+
real mu = E / (2.*(1.+sigma));
20+
real lambda = E*sigma / ((1.+sigma)*(1.-2.*sigma));
21+
22+
// Zmin/zmax
23+
real zmin = 1e100, zmax = -1e100;
24+
for (int i = 0; i < Th.nv; i++) {
25+
zmin = min(zmin, Th(i).z);
26+
zmax = max(zmax, Th(i).z);
27+
}
28+
cout << "zmin=" << zmin << " zmax=" << zmax << endl;
29+
30+
fespace Vh(Th, [P1,P1,P1]);
31+
Vh [ux,uy,uz],[vx,vy,vz];
32+
33+
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
34+
macro Div(u1,u2,u3) (dx(u1)+dy(u2)+dz(u3)) // EOM
35+
36+
37+
problem Elasticite3D([ux,uy,uz],[vx,vy,vz], solver=sparsesolver) =
38+
int3d(Th)(
39+
lambda * Div(vx,vy,vz) * Div(ux,uy,uz)
40+
+ 2.*mu * (Epsilon(vx,vy,vz)' * Epsilon(ux,uy,uz))
41+
)
42+
- int3d(Th)(gravity * vz)
43+
+ on(1, ux=0, uy=0, uz=0); // encastrement
44+
45+
Elasticite3D;
46+
47+
cout << "Deplacement max uz = " << uz[].linfty << " m" << endl;
48+
cout << "Theorique uz_max = " << rho*9.81*(zmax-zmin)^2/(2.*E) << " m" << endl;
49+
50+
fespace Wh(Th, P1);
51+
Wh normu = sqrt(ux^2 + uy^2 + uz^2);
52+
medit("Norme", Th, normu, wait=1);
53+
medit("uz", Th, uz, wait=1);
54+
55+
real umax = max(ux[].linfty, max(uy[].linfty, uz[].linfty));
56+
real coef = 0.1*(zmax-zmin)/umax;
57+
cout << "Coefficient amplification = " << coef << endl;
58+
mesh3 Thd = movemesh3(Th, transfo=[x+coef*ux, y+coef*uy, z+coef*uz]);
59+
//savemesh(Thd, "cube_deforme.mesh");
60+
medit("Cube deforme", Thd, wait=1);
61+
62+
63+
64+
65+
66+
67+
68+
// Exporter la solution en format VTK
69+
int[int] order = [1, 1, 1]; // Ordre P1
70+
savevtk("cube_deplacement.vtu", Th, [ux, uy, uz], dataname="Deplacement", order=order);
71+
savevtk("cube_norme.vtu", Th, normu, dataname="Norme", order=order);
72+
savevtk("cube_uz.vtu", Th, uz, dataname="Uz", order=order);
73+
74+
// Exporter le maillage déformé
75+
mesh3 Thd = movemesh3(Th, transfo=[x + 1e6*ux, y + 1e6*uy, z + 1e6*uz]);
76+
savevtk("cube_deforme.vtu", Thd, [ux, uy, uz], dataname="Deplacement", order=order);
77+
78+
cout << "Fichiers VTK exportés. Ouvrez-les avec ParaView." << endl;

‎examples/Freefem/barre-1d-solution-man.edp‎

Lines changed: 44 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,44 @@
1+
// Solution manufacturée normalisée
2+
// On résout : -u'' = pi^2*sin(pi*x) EA simplifie =1
3+
real L = 1.0;
4+
real E = 1.0;
5+
real A = 1.0;
6+
real pi = 3.14159265358979;
7+
8+
int[int] nvals = [2, 4, 8, 16, 32, 64];
9+
10+
ofstream file("convergence_grossier.txt");
11+
file << "# n\th\terrL2\terrH1" << endl;
12+
13+
for (int idx = 0; idx < nvals.n; idx++) {
14+
int n = nvals[idx];
15+
real h = L / n;
16+
17+
meshL Th = segment(n, [x * L]);
18+
fespace Vh(Th, P1);
19+
Vh u, v;
20+
21+
problem Traction(u, v) =
22+
int1d(Th)( dx(u) * dx(v) )
23+
- int1d(Th)( pi * pi * sin(pi * x) * v )
24+
+ on(1, u = 0)
25+
+ on(2, u = 0);
26+
27+
Traction;
28+
29+
Vh uexact = sin(pi * x);
30+
31+
// Erreur sur les noeuds uniquement
32+
real errL2 = 0, errH1 = 0;
33+
34+
errL2 = sqrt(int1d(Th)( (u - sin(pi*x))^2 ));
35+
errH1 = sqrt(int1d(Th)( (dx(u) - pi*cos(pi*x))^2 ));
36+
37+
38+
cout << "n = " << n
39+
<< ", h = " << h
40+
<< ", errL2 = " << errL2
41+
<< ", errH1 = " << errH1 << endl;
42+
43+
file << n << "\t" << h << "\t" << errL2 << "\t" << errH1 << endl;
44+
}

‎examples/Freefem/barre_exacte.edp‎

Lines changed: 42 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,42 @@
1+
// === Élasticité linéaire - BARRE en traction ===
2+
// Version pour étude de convergence
3+
4+
real L = 1.0;
5+
real E = 1e6;
6+
real A = 0.001;
7+
real F = 1000;
8+
9+
// Différentes valeurs de n pour l'étude de convergence
10+
11+
int[int] nvals = [2, 4, 8, 16, 32, 64];
12+
ofstream file("convergence_data.txt");
13+
file << "# n\th\terrL2\terrH1" << endl;
14+
15+
for (int idx = 0; idx < nvals.n; idx++) {
16+
int n = nvals[idx];
17+
real h = L/n;
18+
19+
meshL Th = segment(n, [x*L]);
20+
fespace Vh(Th, P1);
21+
Vh u, v;
22+
23+
problem Traction(u, v) =
24+
int1d(Th)( E * A * dx(u) * dx(v) )
25+
- int0d(Th, 2)( F * v )
26+
+ on(1, u=0);
27+
28+
Traction;
29+
30+
Vh uexact = (F/(E*A)) * x;
31+
32+
// Erreur L2
33+
real errL2 = sqrt(int1d(Th)((u - uexact)^2));
34+
35+
// Erreur H1 (semi-norme)
36+
real errH1 = sqrt(int1d(Th)((dx(u) - dx(uexact))^2));
37+
38+
cout << "n = " << n << ", h = " << h << ", errL2 = " << errL2 << ", errH1 = " << errH1 << endl;
39+
40+
// Sauvegarde
41+
file << n << "\t " << h << "\t " << errL2 << "\t " << errH1 << endl;
42+
}

‎examples/Freefem/beam-2d.edp‎

Lines changed: 104 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,104 @@
1+
// Beam under its own weight
2+
3+
// Parameters
4+
real E = 21.5;
5+
real sigma = 0.29;
6+
real gravity = -0.05;
7+
int n=140;
8+
int m=35;
9+
10+
11+
// Maillage non structure " avec courbes parametriques "
12+
border a(t=2, 0){x=0; y=t ;label=1;} // gauche cote de la poutre à encastre
13+
border b(t=0, 10){x=t; y=0 ;label=2;} // bas (libre)
14+
border c(t=0, 2){x=10; y=t ;label=3;} // droite (libre)
15+
border d(t=0, 10){x=10-t; y=2; label=4;} // haut (libre)
16+
mesh th = buildmesh(b(n) + c(m) + d(n) + a(m)); // un maillage equilibre previlige celui là
17+
// n=100 unite de long pour les cotes haut et bas, et n=25 unite de haut pour les cote gauche et droite
18+
//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
19+
// Fespace element finis espace P1 ;
20+
plot (th,cmm="maillage-poutre", wait=1);
21+
fespace Vh(th, [P1, P1]);
22+
Vh [uu, vv], [w, s];
23+
24+
// Macro
25+
real sqrt2 = sqrt(2.);
26+
macro epsilon(u1, u2) [dx(u1), dy(u2), (dy(u1) + dx(u2))/sqrt2] // EOM
27+
macro div(u, v) (dx(u) + dy(v)) // EOM
28+
29+
// Coefficients de Lamé
30+
real mu = E/(2*(1 + sigma));
31+
real lambda = E*sigma/((1 + sigma)*(1 - 2*sigma));
32+
33+
// Problème
34+
solve Elasticity ([uu, vv], [w, s], solver=sparsesolver)
35+
= int2d(th)(
36+
lambda*div(w, s)*div(uu, vv)
37+
+ 2.*mu*(epsilon(w, s)' * epsilon(uu, vv))
38+
)
39+
- int2d(th)(gravity*s)
40+
+ on(1, uu=0, vv=0);
41+
42+
// Résultats
43+
cout << "Max displacement = " << uu[].linfty << ", " << vv[].linfty << endl;
44+
45+
// Visualisation
46+
plot([uu, vv], wait=1, cmm="Deplacement");
47+
48+
49+
mesh th1 = movemesh(th, [x + uu, y + vv]);
50+
plot( th1, wait=1, cmm="Maillage initial et deformation");
51+
52+
// Espace pour les contraintes
53+
fespace Wh(th, P1);
54+
Wh sigmaxx, sigmayy, sigmaxy, sigmavm;
55+
56+
// Calcul des contraintes
57+
sigmaxx = lambda * (dx(uu) + dy(vv)) + 2*mu*dx(uu);
58+
sigmayy = lambda * (dx(uu) + dy(vv)) + 2*mu*dy(vv);
59+
sigmaxy = mu*(dy(uu) + dx(vv));
60+
61+
// Contrainte de von Mises pour plasticité
62+
sigmavm = sqrt(sigmaxx^2 + sigmayy^2 - sigmaxx*sigmayy + 3*sigmaxy^2);
63+
64+
// Visualisation
65+
plot(sigmaxx, wait=1, fill=1, value=1, cmm="Contrainte sigma_xx (MPa)");
66+
plot(sigmavm, wait=1, fill=1, value=1, cmm="Contrainte von Mises (MPa)");
67+
68+
// Valeurs extrêmes
69+
cout << "sigma_xx min = " << sigmaxx[].min << ", max = " << sigmaxx[].max << endl;
70+
cout << "sigma_vm max = " << sigmavm[].max << endl;
71+
72+
73+
74+
// Calcul des réactions par intégration sur le bord encastré (label=1)
75+
real reactionx = int1d(th, 1)( (sigmaxx * N.x + sigmaxy * N.y) );
76+
real reactiony = int1d(th, 1)( (sigmaxy * N.x + sigmayy * N.y) );
77+
78+
// Poids total
79+
real poidstotal = abs(gravity) * int2d(th)(1);
80+
81+
// Affichage des résultats
82+
83+
cout << "Reaction horizontale a l'encastrement : " << reactionx << endl;
84+
cout << "Reaction verticale a l'encastrement : " << reactiony << endl;
85+
cout << "Poids total de la poutre : " << poidstotal << endl;
86+
cout << "Rapport reaction verticale / poids : " << reactiony / poidstotal << endl;
87+
88+
// Test d'équilibre
89+
real erreurequilibre = abs(reactiony/poidstotal - 1) * 100;
90+
cout << "Erreur d'equilibre : " << erreurequilibre << " %" << endl;
91+
92+
if (erreurequilibre < 5)
93+
cout << "equilibre verifie " << endl;
94+
else
95+
cout << " Attention : desequilibre detecte !" << endl;
96+
97+
// Vérification du moment à l'encastrement
98+
real moment = int1d(th, 1)( (y - 1) * (sigmaxy * N.x + sigmayy * N.y) - (x) * (sigmaxx * N.x + sigmaxy * N.y) );
99+
cout << "Moment a l'encastrement : " << moment << endl;
100+
101+
102+
103+
104+

‎examples/Freefem/elasticity_3D.edp‎

Lines changed: 105 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,105 @@
1+
load "msh3"
2+
load "medit"
3+
load "iovtk"
4+
// remarque : la comparaison theorique qu'on a obtenue n'est pas vraiment valide
5+
// car on compare une geometrie tres difficle à une poutre, mais on reste tjrs dans le cadre du respect de la loi d la physique
6+
mesh3 Th("C:/Program Files (x86)/FreeFem++/examples/3d/dodecaedre01.mesh");
7+
8+
// Calculer zmin
9+
real zmin = 1e100, zmax = -1e100;
10+
for (int i = 0; i < Th.nv; i++) {
11+
zmin = min(zmin, Th(i).z);
12+
zmax = max(zmax, Th(i).z);
13+
}
14+
cout << "zmin=" << zmin << " zmax=" << zmax << endl;
15+
16+
// Changer le label des faces du bas (proches de zmin)
17+
Th = change(Th, flabel=(abs(z - zmin) < 0.1) ? 1 : 0);
18+
19+
// Vérifier les labels modifiés
20+
cout << "Nouveaux labels :" << endl;
21+
for (int i = 0; i < Th.nbe; i++) {
22+
if (Th.be(i).label == 1)
23+
cout << "Face " << i << " label = " << Th.be(i).label << endl;
24+
}
25+
26+
// Vérifier tous les labels présents
27+
int[int] labs = labels(Th);
28+
cout << "Labels presents dans le maillage :" << endl;
29+
for(int i = 0; i < labs.n; i++)
30+
cout << " label " << labs[i] << endl;
31+
32+
// Paramètres physiques
33+
real E = 210e9;
34+
real sigma = 0.3;
35+
real rho = 7800.;
36+
real gravity = -9.81 * rho;
37+
real mu = E / (2.*(1.+sigma));
38+
real lambda = E*sigma / ((1.+sigma)*(1.-2.*sigma));
39+
40+
// Espace éléments finis
41+
fespace Vh(Th, [P1,P1,P1]);
42+
Vh [ux,uy,uz],[vx,vy,vz];
43+
44+
// Macros (attention aux parenthèses superflues)
45+
macro Epsilon(u1,u2,u3) [dx(u1), dy(u2), dz(u3),
46+
(dy(u3)+dz(u2))/sqrt(2.),
47+
(dz(u1)+dx(u3))/sqrt(2.),
48+
(dx(u2)+dy(u1))/sqrt(2.)] // EOM
49+
macro Div(u1,u2,u3) (dx(u1) + dy(u2) + dz(u3)) // EOM
50+
51+
// Problème d'élasticité on fixe les faces label=1
52+
problem Elasticite3D([ux,uy,uz],[vx,vy,vz], solver=sparsesolver) =
53+
int3d(Th)(
54+
lambda * Div(vx,vy,vz) * Div(ux,uy,uz)
55+
+ 2.*mu * (Epsilon(vx,vy,vz)' * Epsilon(ux,uy,uz))
56+
)
57+
- int3d(Th)(gravity * vz)
58+
+ on(1, ux=0, uy=0, uz=0);
59+
60+
// Résolution
61+
Elasticite3D;
62+
63+
// Résultats
64+
cout << "Deplacement max |ux| = " << ux[].linfty << " m" << endl;
65+
cout << "Deplacement max |uy| = " << uy[].linfty << " m" << endl;
66+
cout << "Deplacement max |uz| = " << uz[].linfty << " m" << endl;
67+
68+
// Création d'un espace scalaire P1
69+
fespace Wh(Th, P1);
70+
Wh normu;
71+
72+
73+
normu = sqrt(ux*ux + uy*uy + uz*uz);
74+
75+
cout << "Deplacement max total = " << normu[].linfty << " m" << endl;
76+
77+
78+
79+
cout << "Theorique (poutre) = " << rho*9.81*(zmax-zmin)^2/(2.*E) << " m" << endl;
80+
81+
// Visualisation
82+
medit("Norme du deplacement", Th, normu, wait=1);
83+
medit("Deplacement vertical uz", Th, uz, wait=1);
84+
85+
// Maillage déformé
86+
real umax = max(ux[].linfty, max(uy[].linfty, uz[].linfty));
87+
real coef = 0.1*(zmax-zmin)/umax;
88+
mesh3 Thd = movemesh3(Th, transfo=[x+coef*ux, y+coef*uy, z+coef*uz]);
89+
medit("Dodecaedre deforme", Thd, wait=1);
90+
91+
92+
93+
94+
95+
96+
// Exporter la solution en format VTK
97+
//int[int] order = [1, 1, 1];
98+
//savevtk("cube_deplacement.vtu", Th, [ux, uy, uz], dataname="Deplacement", order=order);
99+
//savevtk("cube_norme.vtu", Th, normu, dataname="Norme", order=order);
100+
//savevtk("cube_uz.vtu", Th, uz, dataname="Uz", order=order);
101+
102+
// Exporter le maillage déformé
103+
//mesh3 Thd = movemesh3(Th, transfo=[x + 1e6*ux, y + 1e6*uy, z + 1e6*uz]);
104+
//savevtk("cube_deforme.vtu", Thd, [ux, uy, uz], dataname="Deplacement", order=order);
105+

0 commit comments

Comments
 (0)