FreeFEM Documentation on GitHub

stars - forks

Misc

Poisson’s Equation

 1 // Parameters
 2 int nn = 20;
 3 real L = 1.;
 4 real H = 1.;
 5 real l = 0.5;
 6 real h = 0.5;
 7 
 8 func f = 1.;
 9 func g = 0.;
10 
11 int NAdapt = 10;
12 
13 // Mesh
14 border b1(t=0, L){x=t; y=0;};
15 border b2(t=0, h){x=L; y=t;};
16 border b3(t=L, l){x=t; y=h;};
17 border b4(t=h, H){x=l; y=t;};
18 border b5(t=l, 0){x=t; y=H;};
19 border b6(t=H, 0){x=0; y=t;};
20 
21 mesh Th = buildmesh(b1(nn*L) + b2(nn*h) + b3(nn*(L-l)) + b4(nn*(H-h)) + b5(nn*l) + b6(nn*H));
22 
23 // Fespace
24 fespace Vh(Th, P1); // Change P1 to P2 to test P2 finite element
25 Vh u, v;
26 
27 // Macro
28 macro grad(u) [dx(u), dy(u)] //
29 
30 // Problem
31 problem Poisson (u, v, solver=CG, eps=-1.e-6)
32    = int2d(Th)(
33         grad(u)' * grad(v)
34    )
35    + int2d(Th)(
36         f * v
37    )
38    + on(b1, b2, b3, b4, b5, b6, u=g)
39    ;
40 
41 // Mesh adaptation iterations
42 real error = 0.1;
43 real coef = 0.1^(1./5.);
44 for (int i = 0; i < NAdapt; i++){
45    // Solve
46    Poisson;
47 
48    // Plot
49    plot(Th, u);
50 
51    // Adaptmesh
52    Th = adaptmesh(Th, u, inquire=1, err=error);
53    error = error * coef;
54 }
PoissonAssociatedMesh

Fig. 181 Adapted mesh

PoissonAdaptedMesh

Fig. 182 Solution on adapted mesh

Poisson

Poisson’s equation 3D

 1 load "tetgen"
 2 
 3 // Parameters
 4 real hh = 0.1;
 5 func ue = 2.*x*x + 3.*y*y + 4.*z*z + 5.*x*y + 6.*x*z + 1.;
 6 func f= -18.;
 7 
 8 // Mesh
 9 mesh Th = square(10, 20, [x*pi-pi/2, 2*y*pi]); // ]-pi/2, pi/2[X]0,2pi[
10 func f1 = cos(x)*cos(y);
11 func f2 = cos(x)*sin(y);
12 func f3 = sin(x);
13 func f1x = sin(x)*cos(y);
14 func f1y = -cos(x)*sin(y);
15 func f2x = -sin(x)*sin(y);
16 func f2y = cos(x)*cos(y);
17 func f3x = cos(x);
18 func f3y = 0;
19 func m11 = f1x^2 + f2x^2 + f3x^2;
20 func m21 = f1x*f1y + f2x*f2y + f3x*f3y;
21 func m22 = f1y^2 + f2y^2 + f3y^2;
22 func perio = [[4, y], [2, y], [1, x], [3, x]];
23 real vv = 1/square(hh);
24 Th = adaptmesh(Th, m11*vv, m21*vv, m22*vv, IsMetric=1, periodic=perio);
25 Th = adaptmesh(Th, m11*vv, m21*vv, m22*vv, IsMetric=1, periodic=perio);
26 plot(Th);
27 
28 real[int] domain = [0., 0., 0., 1, 0.01];
29 mesh3 Th3 = tetgtransfo(Th, transfo=[f1, f2, f3], nbofregions=1, regionlist=domain);
30 plot(Th3);
31 
32 border cc(t=0, 2*pi){x=cos(t); y=sin(t); label=1;}
33 mesh Th2 = buildmesh(cc(50));
34 
35 // Fespace
36 fespace Vh(Th3, P23d);
37 Vh u, v;
38 Vh uhe = ue;
39 cout << "uhe min: " << uhe[].min << " - max: " << uhe[].max << endl;
40 cout << uhe(0.,0.,0.) << endl;
41 
42 fespace Vh2(Th2, P2);
43 Vh2 u2, u2e;
44 
45 // Macro
46 macro Grad3(u) [dx(u), dy(u), dz(u)] //
47 
48 // Problem
49 problem Lap3d (u, v, solver=CG)
50     = int3d(Th3)(
51           Grad3(v)' * Grad3(u)
52     )
53     - int3d(Th3)(
54           f * v
55     )
56     + on(0, 1, u=ue)
57     ;
58 
59 // Solve
60 Lap3d;
61 cout << "u min: " << u[]. min << " - max: " << u[].max << endl;
62 
63 // Error
64 real err = int3d(Th3)(square(u-ue));
65 cout << int3d(Th3)(1.) << " = " << Th3.measure << endl;
66 Vh d = ue - u;
67 cout << " err = " << err << " - diff l^intfy = " << d[].linfty << endl;
68 
69 // Plot
70 u2 = u;
71 u2e = ue;
72 plot(u2, wait=true);
73 plot(u2, u2e,wait=true);
../_images/poisson_3d.jpg

Fig. 183 Iso-surfaces of the solution

Stokes Equation on a cube

 1 load "msh3"
 2 load "medit" // Dynamically loaded tools for 3D
 3 
 4 // Parameters
 5 int nn = 8;
 6 
 7 // Mesh
 8 mesh Th0 = square(nn, nn);
 9 int[int] rup = [0, 2];
10 int[int] rdown = [0, 1];
11 int[int] rmid = [1, 1, 2, 1, 3, 1, 4, 1];
12 real zmin = 0, zmax = 1;
13 mesh3 Th = buildlayers(Th0, nn, zbound=[zmin, zmax],
14     reffacemid=rmid, reffaceup=rup, reffacelow=rdown);
15 
16 medit("c8x8x8", Th); // 3D mesh visualization with medit
17 
18 // Fespaces
19 fespace Vh2(Th0, P2);
20 Vh2 ux, uz, p2;
21 
22 fespace VVh(Th, [P2, P2, P2, P1]);
23 VVh [u1, u2, u3, p];
24 VVh [v1, v2, v3, q];
25 
26 // Macro
27 macro Grad(u) [dx(u), dy(u), dz(u)] //
28 macro div(u1,u2,u3) (dx(u1) + dy(u2) + dz(u3)) //
29 
30 // Problem (directly solved)
31 solve vStokes ([u1, u2, u3, p], [v1, v2, v3, q])
32     = int3d(Th, qforder=3)(
33           Grad(u1)' * Grad(v1)
34         + Grad(u2)' * Grad(v2)
35         + Grad(u3)' * Grad(v3)
36         - div(u1, u2, u3) * q
37         - div(v1, v2, v3) * p
38         + 1e-10 * q * p
39     )
40     + on(2, u1=1., u2=0, u3=0)
41     + on(1, u1=0, u2=0, u3=0)
42     ;
43 
44 // Plot
45 plot(p, wait=1, nbiso=5); // 3D visualization of pressure isolines
46 
47 // See 10 plan of the velocity in 2D
48 for(int i = 1; i < 10; i++){
49     // Cut plane
50     real yy = i/10.;
51     // 3D to 2D interpolation
52     ux = u1(x,yy,y);
53     uz = u3(x,yy,y);
54     p2 = p(x,yy,y);
55     // Plot
56     plot([ux, uz], p2, cmm="cut y = "+yy, wait= 1);
57 }
Stokes3d

Fig. 184 Solution

Stokes3d-Th

Fig. 185 Associated mesh

Stokes

Cavity

 1 //Parameters
 2 int m = 300;
 3 real L = 1;
 4 real rho = 500.;
 5 real mu = 0.1;
 6 
 7 real uin = 1;
 8 func fx = 0;
 9 func fy = 0;
10 int[int] noslip = [1, 2, 4];
11 int[int] inflow = [3];
12 
13 real dt = 0.1;
14 real T = 50;
15 
16 real eps = 1e-3;
17 
18 //Macros
19 macro div(u) (dx(u#x) + dy(u#y))//
20 macro grad(u) [dx(u), dy(u)]//
21 macro Grad(u) [grad(u#x), grad(u#y)]//
22 
23 //Time
24 real cpu;
25 real tabcpu;
26 
27 //mesh
28 border C1(t = 0, L){ x = t; y = 0; label = 1; }
29 border C2(t = 0, L){ x = L; y = t; label = 2; }
30 border C3(t = 0, L){ x = L-t; y = L; label = 3; }
31 border C4(t = 0, L){ x = 0; y = L-t; label = 4; }
32 mesh th = buildmesh( C1(m) + C2(m) + C3(m) + C4(m) );
33 
34 fespace UPh(th, [P2,P2,P1]);
35 UPh [ux, uy, p];
36 UPh [uhx, uhy, ph];
37 UPh [upx, upy, pp];
38 
39 //Solve
40 varf navierstokes([ux, uy, p], [uhx, uhy, ph])
41   = int2d(th)(
42       rho/dt* [ux, uy]'* [uhx, uhy]
43     + mu* (Grad(u):Grad(uh))
44     - p* div(uh)
45     - ph* div(u)
46     - 1e-10 *p*ph
47     )
48 
49   + int2d(th) (
50       [fx, fy]' * [uhx, uhy]
51     + rho/dt* [convect([upx, upy], -dt, upx), convect([upx, upy], -dt, upy)]'* [uhx, uhy]
52     )
53 
54   + on(noslip, ux=0, uy=0)
55   + on(inflow, ux=uin, uy=0)
56   ;
57 
58 //Initialization
59 [ux, uy, p]=[0, 0, 0];
60 
61 matrix<real> NS = navierstokes(UPh, UPh, solver=sparsesolver);
62 real[int] NSrhs = navierstokes(0, UPh);
63 
64 //Time loop
65 for(int j = 0; j < T/dt; j++){
66   [upx, upy, pp]=[ux, uy, p];
67 
68   NSrhs = navierstokes(0, UPh);
69   ux[] = NS^-1 * NSrhs;
70 
71   plot( [ux,uy], p, wait=0, cmm=j);
72 }
73 
74 //CPU
75 cout << " CPU = " << clock()-cpu << endl ;
76 tabcpu = clock()-cpu;
Table of content