FreeFEM Documentation on GitHub

stars - forks

Mesh Generation

Square mesh

 1 mesh Th0 = square(10 ,10);
 2 
 3 mesh Th1 = square(4, 5);
 4 
 5 real x0 = 1.2;
 6 real x1 = 1.8;
 7 real y0 = 0;
 8 real y1 = 1;
 9 int n = 5;
10 real m = 20;
11 mesh Th2 = square(n, m, [x0+(x1-x0)*x, y0+(y1-y0)*y]);
12 
13 for (int i = 0; i < 5; ++i){
14    int[int] labs = [11, 12, 13, 14];
15    mesh Th3 = square(3, 3, flags=i, label=labs, region=10);
16    plot(Th3, wait=1, cmm="square flags = "+i );
17 }

Mesh adaptation

 1 // Parameters
 2 real eps = 0.0001;
 3 real h = 1;
 4 real hmin = 0.05;
 5 func f = 10.0*x^3 + y^3 + h*atan2(eps, sin(5.0*y)-2.0*x);
 6 
 7 // Mesh
 8 mesh Th = square(5, 5, [-1+2*x, -1+2*y]);
 9 
10 // Fespace
11 fespace Vh(Th,P1);
12 Vh fh = f;
13 plot(fh);
14 
15 // Adaptmesh
16 for (int i = 0; i < 2; i++){
17    Th = adaptmesh(Th, fh);
18    fh = f; //old mesh is deleted
19    plot(Th, fh, wait=true);
20 }
MeshAdaptation1

Fig. 186 Initial mesh

MeshAdaptation2

Fig. 187 Adapted mesh

Mesh adaptation

Mesh adaptation for the Poisson’s problem

 1 // Parameters
 2 real error = 0.1;
 3 
 4 // Mesh
 5 border ba(t=0, 1){x=t; y=0; label=1;}
 6 border bb(t=0, 0.5){x=1; y=t; label=1;}
 7 border bc(t=0, 0.5){x=1-t; y=0.5; label=1;}
 8 border bd(t=0.5, 1){x=0.5; y=t; label=1;}
 9 border be(t=0.5, 1){x=1-t; y=1; label=1;}
10 border bf(t=0, 1){x=0; y=1-t; label=1;}
11 mesh Th = buildmesh(ba(6) + bb(4) + bc(4) + bd(4) + be(4) + bf(6));
12 
13 // Fespace
14 fespace Vh(Th, P1);
15 Vh u, v;
16 
17 // Function
18 func f = 1;
19 
20 // Problem
21 problem Poisson(u, v, solver=CG, eps=1.e-6)
22    = int2d(Th)(
23         dx(u)*dx(v)
24       + dy(u)*dy(v)
25    )
26    - int2d(Th)(
27         f*v
28    )
29    + on(1, u=0);
30 
31 // Adaptmesh loop
32 for (int i = 0; i < 4; i++){
33    Poisson;
34    Th = adaptmesh(Th, u, err=error);
35    error = error/2;
36 }
37 
38 // Plot
39 plot(u);
MeshAdaptationPoisson1

Fig. 188 Initial mesh

MeshAdaptationPoisson2

Fig. 189 Adapted mesh

MeshAdaptationPoissonU

Fig. 190 Solution on adapted mesh

Mesh adaptation (Poisson)

Uniform mesh adaptation

1 mesh Th = square(2, 2); // The initial mesh
2 plot(Th, wait=true);
3 
4 Th = adaptmesh(Th, 1./30., IsMetric=1, nbvx=10000);
5 plot(Th, wait=true);
6 
7 Th = adaptmesh(Th, 1./30., IsMetric=1, nbvx=10000); // More than one time due to the
8 Th = adaptmesh(Th, 1./30., IsMetric=1, nbvx=10000); // adaptation bound `maxsubdiv=`
9 plot(Th, wait=true);
UniformMeshAdaptation1

Fig. 191 Initial mesh

UniformMeshAdaptation2

Fig. 192 Adapted mesh

Uniform mesh adaptation

Borders

 1 {
 2    int upper = 1;
 3    int others = 2;
 4    int inner = 3;
 5 
 6    border C01(t=0, 1){x=0; y=-1+t; label=upper;}
 7    border C02(t=0, 1){x=1.5-1.5*t; y=-1; label=upper;}
 8    border C03(t=0, 1){x=1.5; y=-t; label=upper;}
 9    border C04(t=0, 1){x=1+0.5*t; y=0; label=others;}
10    border C05(t=0, 1){x=0.5+0.5*t; y=0; label=others;}
11    border C06(t=0, 1){x=0.5*t; y=0; label=others;}
12    border C11(t=0, 1){x=0.5; y=-0.5*t; label=inner;}
13    border C12(t=0, 1){x=0.5+0.5*t; y=-0.5; label=inner;}
14    border C13(t=0, 1){x=1; y=-0.5+0.5*t; label=inner;}
15 
16    int n = 10;
17    plot(C01(-n) + C02(-n) + C03(-n) + C04(-n) + C05(-n)
18       + C06(-n) + C11(n) + C12(n) + C13(n), wait=true);
19 
20    mesh Th = buildmesh(C01(-n) + C02(-n) + C03(-n) + C04(-n) + C05(-n)
21       + C06(-n) + C11(n) + C12(n) + C13(n));
22 
23    plot(Th, wait=true);
24 
25    cout << "Part 1 has region number " << Th(0.75, -0.25).region << endl;
26    cout << "Part 2 has redion number " << Th(0.25, -0.25).region << endl;
27 }
28 
29 {
30    border a(t=0, 2*pi){x=cos(t); y=sin(t); label=1;}
31    border b(t=0, 2*pi){x=0.3+0.3*cos(t); y=0.3*sin(t); label=2;}
32    plot(a(50) + b(30)); //to see a plot of the border mesh
33    mesh Thwithouthole = buildmesh(a(50) + b(30));
34    mesh Thwithhole = buildmesh(a(50) + b(-30));
35    plot(Thwithouthole);
36    plot(Thwithhole);
37 }
38 
39 {
40    real r=1;
41    border a(t=0, 2*pi){x=r*cos(t); y=r*sin(t); label=1;}
42    r=0.3;
43    border b(t=0, 2*pi){x=r*cos(t); y=r*sin(t); label=1;}
44    //  mesh Thwithhole = buildmesh(a(50) + b(-30)); // do not do this because the two
45    // circles have the same radius = $0.3$
46 }
Borders1

Fig. 193 Mesh with two regions

Borders2

Fig. 194 Mesh without a hole

Borders3

Fig. 195 Mesh with a hole

Borders

Change

 1 verbosity=3;
 2 
 3 // Mesh
 4 mesh Th1 = square(10, 10);
 5 mesh Th2 = square(20, 10, [x+1, y]);
 6 
 7 int[int] r1=[2, 0];
 8 plot(Th1, wait=true);
 9 
10 Th1 = change(Th1, label=r1); // Change edges' label from 2 to 0
11 plot(Th1, wait=true);
12 
13 int[int] r2=[4, 0];
14 Th2 = change(Th2, label=r2); // Change edges' label from 4 to 0
15 plot(Th2, wait=true);
16 
17 mesh Th = Th1 + Th2; // 'gluing together' Th1 and Th2 meshes
18 cout << "nb lab = " << int1d(Th1,1,3,4)(1./lenEdge)+int1d(Th2,1,2,3)(1./lenEdge)
19    << " == " << int1d(Th,1,2,3,4)(1./lenEdge) << " == " << ((10+20)+10)*2 << endl;
20 plot(Th, wait=true);
21 
22 fespace Vh(Th, P1);
23 Vh u, v;
24 
25 macro Grad(u) [dx(u),dy(u)] // Definition of a macro
26 
27 solve P(u, v)
28    = int2d(Th)(
29         Grad(u)'*Grad(v)
30    )
31    -int2d(Th)(
32         v
33    )
34    + on(1, 3, u=0)
35    ;
36 
37 plot(u, wait=1);
../_images/Change.jpg

Fig. 196 Result

Cube

 1 load "msh3"
 2 
 3 int[int] l6 = [37, 42, 45, 40, 25, 57];
 4 int r11 = 11;
 5 mesh3 Th = cube(4, 5, 6, [x*2-1, y*2-1, z*2-1], label=l6, flags =3, region=r11);
 6 
 7 cout << "Volume = " << Th.measure << ", border area = " << Th.bordermeasure << endl;
 8 
 9 int err = 0;
10 for(int i = 0; i < 100; ++i){
11    real s = int2d(Th,i)(1.);
12    real sx = int2d(Th,i)(x);
13    real sy = int2d(Th,i)(y);
14    real sz = int2d(Th,i)(z);
15 
16    if(s){
17       int ix = (sx/s+1.5);
18       int iy = (sy/s+1.5);
19       int iz = (sz/s+1.5);
20       int ii = (ix + 4*(iy+1) + 16*(iz+1) );
21       //value of ix,iy,iz => face min 0, face max 2, no face 1
22       cout << "Label = " << i << ", s = " << s << " " << ix << iy << iz << " : " << ii << endl;
23       if( i != ii ) err++;
24    }
25 }
26 real volr11 = int3d(Th,r11)(1.);
27 cout << "Volume region = " << 11 << ": " << volr11 << endl;
28 if((volr11 - Th.measure )>1e-8) err++;
29 plot(Th, fill=false);
30 cout << "Nb err = " << err << endl;
31 assert(err==0);
../_images/Cube.jpg

Fig. 197 Cube

Empty mesh

 1 {
 2    border a(t=0, 2*pi){x=cos(t); y=sin(t); label=1;}
 3    mesh Th = buildmesh(a(20));
 4    Th = emptymesh(Th);
 5    plot(Th);
 6 }
 7 {
 8    mesh Th = square(10, 10);
 9    int[int] ssd(Th.nt);
10    // Builds the pseudo region numbering
11    for(int i = 0; i < ssd.n; i++){
12       int iq = i/2; // Because we have 2 triangles per quad
13       int ix = iq%10;
14       int iy = iq/10;
15       ssd[i] = 1 + (ix>=5) + (iy>=5)*2;
16    }
17    // Builds an emtpy mesh with all edges that satisfy e=T1 cap T2 and ssd[T1] != ssd[T2]
18    Th = emptymesh(Th, ssd);
19    // Plot
20    plot(Th);
21 }
EmptyMesh1

Fig. 198 Empty square

EmptyMesh2

Fig. 199 Empty diamond

Empty mesh

3 points

 1 // Square for Three-Point Bend Specimens fixed on Fix1, Fix2
 2 // It will be loaded on Load
 3 real a = 1, b = 5, c = 0.1;
 4 int n = 5, m = b*n;
 5 border Left(t=0, 2*a){x=-b; y=a-t;}
 6 border Bot1(t=0, b/2-c){x=-b+t; y=-a;}
 7 border Fix1(t=0, 2*c){x=-b/2-c+t; y=-a;}
 8 border Bot2(t=0, b-2*c){x=-b/2+c+t; y=-a;}
 9 border Fix2(t=0, 2*c){x=b/2-c+t; y=-a;}
10 border Bot3(t=0, b/2-c){x=b/2+c+t; y=-a;}
11 border Right(t=0, 2*a){x=b; y=-a+t;}
12 border Top1(t=0, b-c){x=b-t; y=a;}
13 border Load(t=0, 2*c){x=c-t; y=a;}
14 border Top2(t=0, b-c){x=-c-t; y=a;}
15 
16 mesh Th = buildmesh(Left(n) + Bot1(m/4) + Fix1(5) + Bot2(m/2)
17    + Fix2(5) + Bot3(m/4) + Right(n) + Top1(m/2) + Load(10) + Top2(m/2));
18 plot(Th, bw=true);
../_images/3Points.jpg

Fig. 200 3 Points

Bezier

 1 // A cubic Bezier curve connecting two points with two control points
 2 func real bzi(real p0, real p1, real q1, real q2, real t){
 3    return p0*(1-t)^3 + q1*3*(1-t)^2*t + q2*3*(1-t)*t^2 + p1*t^3;
 4 }
 5 
 6 real[int] p00 = [0, 1], p01 = [0, -1], q00 = [-2, 0.1], q01 = [-2, -0.5];
 7 real[int] p11 = [1,-0.9], q10 = [0.1, -0.95], q11=[0.5, -1];
 8 real[int] p21 = [2, 0.7], q20 = [3, -0.4], q21 = [4, 0.5];
 9 real[int] q30 = [0.5, 1.1], q31 = [1.5, 1.2];
10 border G1(t=0, 1){
11    x=bzi(p00[0], p01[0], q00[0], q01[0], t);
12    y=bzi(p00[1], p01[1], q00[1], q01[1], t);
13 }
14 border G2(t=0, 1){
15    x=bzi(p01[0], p11[0], q10[0], q11[0], t);
16    y=bzi(p01[1], p11[1], q10[1], q11[1], t);
17 }
18 border G3(t=0, 1){
19    x=bzi(p11[0], p21[0], q20[0], q21[0], t);
20    y=bzi(p11[1], p21[1], q20[1], q21[1], t);
21 }
22 border G4(t=0, 1){
23    x=bzi(p21[0], p00[0], q30[0], q31[0], t);
24    y=bzi(p21[1], p00[1], q30[1], q31[1], t);
25 }
26 int m = 5;
27 mesh Th = buildmesh(G1(2*m) + G2(m) + G3(3*m) + G4(m));
28 plot(Th, bw=true);
../_images/Bezier.jpg

Fig. 201 Bezier

Build layer mesh

 1 load "msh3"
 2 load "tetgen"
 3 load "medit"
 4 
 5 // Parameters
 6 int C1 = 99;
 7 int C2 = 98;
 8 
 9 // 2D mesh
10 border C01(t=0, pi){x=t; y=0; label=1;}
11 border C02(t=0, 2*pi){ x=pi; y=t; label=1;}
12 border C03(t=0, pi){ x=pi-t; y=2*pi; label=1;}
13 border C04(t=0, 2*pi){ x=0; y=2*pi-t; label=1;}
14 
15 border C11(t=0, 0.7){x=0.5+t; y=2.5; label=C1;}
16 border C12(t=0, 2){x=1.2; y=2.5+t; label=C1;}
17 border C13(t=0, 0.7){x=1.2-t; y=4.5; label=C1;}
18 border C14(t=0, 2){x=0.5; y=4.5-t; label=C1;}
19 
20 border C21(t=0, 0.7){x=2.3+t; y=2.5; label=C2;}
21 border C22(t=0, 2){x=3; y=2.5+t; label=C2;}
22 border C23(t=0, 0.7){x=3-t; y=4.5; label=C2;}
23 border C24(t=0, 2){x=2.3; y=4.5-t; label=C2;}
24 
25 mesh Th = buildmesh(C01(10) + C02(10) + C03(10) + C04(10)
26    + C11(5) + C12(5) + C13(5) + C14(5)
27    + C21(-5) + C22(-5) + C23(-5) + C24(-5));
28 
29 mesh Ths = buildmesh(C01(10) + C02(10) + C03(10) + C04(10)
30    + C11(5) + C12(5) + C13(5) + C14(5));
31 
32 // Construction of a box with one hole and two regions
33 func zmin = 0.;
34 func zmax = 1.;
35 int MaxLayer = 10;
36 
37 func XX = x*cos(y);
38 func YY = x*sin(y);
39 func ZZ = z;
40 
41 int[int] r1 = [0, 41], r2 = [98, 98, 99, 99, 1, 56];
42 int[int] r3 = [4, 12]; // Change upper surface mesh's triangles labels
43 // generated by the 2D mesh's triangles Th
44 // from label 4 to label 12
45 int[int] r4 = [4, 45]; // Change lower surface mesh's triangles labels
46 // generated by the 2D mesh's triangles Th
47 // from label 4 to label 45
48 
49 mesh3 Th3 = buildlayers(Th, MaxLayer, zbound=[zmin, zmax], region=r1,
50    labelmid=r2, labelup=r3, labeldown=r4);
51 medit("box 2 regions 1 hole", Th3);
52 
53 // Construction of a sphere with TetGen
54 func XX1 = cos(y)*sin(x);
55 func YY1 = sin(y)*sin(x);
56 func ZZ1 = cos(x);
57 
58 real[int] domain = [0., 0., 0., 0, 0.001];
59 string test = "paACQ";
60 cout << "test = " << test << endl;
61 mesh3 Th3sph = tetgtransfo(Ths, transfo=[XX1, YY1, ZZ1],
62    switch=test, nbofregions=1, regionlist=domain);
63 medit("sphere 2 regions", Th3sph);
BuildLayerMesh1

Fig. 202 Box with a hole

BuildLayerMesh2

Fig. 203 Sphere

Build layer mesh

Sphere

 1 // Parameter
 2 real hh = 0.1;
 3 
 4 // Mesh 2D
 5 mesh Th = square(10, 20, [x*pi-pi/2, 2*y*pi]); // ]-pi/2, pi/2[X]0, 2pi[
 6 // A parametrization of a sphere
 7 func f1 = cos(x)*cos(y);
 8 func f2 = cos(x)*sin(y);
 9 func f3 = sin(x);
10 // Partial derivative of the parametrization DF
11 func f1x = sin(x)*cos(y);
12 func f1y = -cos(x)*sin(y);
13 func f2x = -sin(x)*sin(y);
14 func f2y = cos(x)*cos(y);
15 func f3x = cos(x);
16 func f3y = 0;
17 //M = DF^t DF
18 func m11 = f1x^2 + f2x^2 + f3x^2;
19 func m21 = f1x*f1y + f2x*f2y + f3x*f3y;
20 func m22 = f1y^2 + f2y^2 + f3y^2;
21 
22 // Periodic condition
23 func perio = [[4, y], [2, y], [1, x], [3, x]];
24 
25 // Mesh adaptation
26 real vv = 1/square(hh);
27 Th = adaptmesh(Th, m11*vv, m21*vv, m22*vv, IsMetric=1, inquire=1, periodic=perio);
28 Th = adaptmesh(Th, m11*vv, m21*vv, m22*vv, IsMetric=1, periodic=perio);
29 Th = adaptmesh(Th, m11*vv, m21*vv, m22*vv, IsMetric=1, periodic=perio);
30 Th = adaptmesh(Th, m11*vv, m21*vv, m22*vv, IsMetric=1, periodic=perio);
31 
32 // Sphere
33 mesh3 Th3 = movemesh23(Th, transfo=[f1, f2, f3]);
34 plot(Th3);
Sphere1

Fig. 204 Initial mesh

Sphere2

Fig. 205 Sphere

Sphere

Table of content