FreeFEM Documentation on GitHub

stars - forks

# Developers

## FFT

 1load "dfft"
2
3// Parameters
4int nx = 32;
5real ny = 16;
6real N = nx*ny;
7func f1 = cos(2*x*2*pi)*cos(3*y*2*pi);
8
9// Mesh
10//warning: the fourier space is not exactly the unit square due to periodic condition
11mesh Th = square(nx-1, ny-1, [(nx-1)*x/nx, (ny-1)*y/ny]);
12//warning: the numbering of the vertices (x,y) is
13//given by i = x/nx + nx*y/ny
14
15// Fespace
16fespace Vh(Th,P1);
17Vh<complex> u = f1, v;
18Vh w = f1;
19Vh ur, ui;
20
21// FFT
22//in dfft the matrix n, m is in row-major order and array n, m is
23//store j + m*i (the transpose of the square numbering)
24v[] = dfft(u[], ny, -1);
25u[] = dfft(v[], ny, +1);
26cout << "||u||_\infty " << u[].linfty << endl;
27
28u[] *= 1./N;
29cout << "||u||_\infty " << u[].linfty << endl;
30
31ur = real(u);
32
33// Plot
34plot(w, wait=1, value=1, cmm="w");
35plot(ur, wait=1, value=1, cmm="u");
36v = w - u;
37cout << "diff = " << v[].max << " " << v[].min << endl;
38assert( norm(v[].max) < 1e-10 && norm(v[].min) < 1e-10);
39
40// Other example
41//FFT Lapacian
42//-\Delta u = f with biperiodic condition
43func f = cos(3*2*pi*x)*cos(2*2*pi*y);
44func ue = (1./(square(2*pi)*13.))*cos(3*2*pi*x)*cos(2*2*pi*y); //the exact solution
45Vh<complex> ff = f;
46Vh<complex> fhat;
47Vh<complex> wij;
48
49// FFT
50fhat[] = dfft(ff[],ny,-1);
51
52//warning in fact we take mode between -nx/2, nx/2 and -ny/2, ny/2
53//thanks to the operator ?:
54wij = square(2.*pi)*(square(( x<0.5?x*nx:(x-1)*nx)) + square((y<0.5?y*ny:(y-1)*ny)));
55wij[][0] = 1e-5; //to remove div / 0
56fhat[] = fhat[] ./ wij[];
57u[] = dfft(fhat[], ny, 1);
58u[] /= complex(N);
59ur = real(u); //the solution
60w = real(ue); //the exact solution
61
62// Plot
63plot(w, ur, value=1, cmm="ue", wait=1);
64
65// Error
66w[] -= ur[];
67real err = abs(w[].max) + abs(w[].min);
68cout << "err = " << err << endl;
69assert(err < 1e-6);
70
71fftwplan p1 = plandfft(u[], v[], ny, -1);
72fftwplan p2 = plandfft(u[], v[], ny, 1);
73real ccc = square(2.*pi);
74cout << "ny = " << ny << endl;
75map(wij[], ny, ccc*(x*x+y*y));
76wij[][0] = 1e-5;
77plot(wij, cmm="wij");


## Complex

 1real a = 2.45, b = 5.33;
2complex z1 = a + b*1i, z2 = a + sqrt(2.)*1i;
3
4func string pc(complex z){
5    string r = "(" + real(z);
6    if (imag(z) >= 0) r = r + "+";
7    return r + imag(z) + "i)";
8}
9
10func string toPolar(complex z){
11    return "";//abs(z) + "*(cos(" + arg(z) + ")+i*sin(" + arg(z) + "))";
12}
13
14cout << "Standard output of the complex " << pc(z1) << " is the pair: " << z1 << endl;
15cout << pc(z1) << " + " << pc(z2) << " = " << pc(z1+z2) << endl;
16cout << pc(z1) << " - " << pc(z2) << " = " << pc(z1-z2) << endl;
17cout << pc(z1) << " * " << pc(z2) << " = " << pc(z1*z2) << endl;
18cout << pc(z1) << " / " << pc(z2) << " = " << pc(z1/z2) << endl;
19cout << "Real part of " << pc(z1) << " = " << real(z1) << endl;
20cout << "Imaginary part of " << pc(z1) << " = " << imag(z1) << endl;
21cout << "abs(" << pc(z1) << ") = " << abs(z1) << endl;
22cout << "Polar coordinates of " << pc(z2) << " = " << toPolar(z2) << endl;
23cout << "de Moivre formula: " << pc(z2) << "^3 = " << toPolar(z2^3) << endl;
24cout << " and polar(" << abs(z2) << ", " << arg(z2) << ") = " << pc(polar(abs(z2), arg(z2))) << endl;
25cout << "Conjugate of " <<pc(z2) << " = " << pc(conj(z2)) <<endl;
26cout << pc(z1) << " ^ " << pc(z2) << " = " << pc(z1^z2) << endl;


Output of this script is:

 1Standard output of the complex (2.45+5.33i) is the pair: (2.45,5.33)
2(2.45+5.33i) + (2.45+1.41421i) = (4.9+6.74421i)
3(2.45+5.33i) - (2.45+1.41421i) = (0+3.91579i)
4(2.45+5.33i) * (2.45+1.41421i) = (-1.53526+16.5233i)
5(2.45+5.33i) / (2.45+1.41421i) = (1.692+1.19883i)
6Real part of (2.45+5.33i) = 2.45
7Imaginary part of (2.45+5.33i) = 5.33
8abs((2.45+5.33i)) = 5.86612
9Polar coordinates of (2.45+1.41421i) =
10de Moivre formula: (2.45+1.41421i)^3 =
11 and polar(2.82887, 0.523509) = (2.45+1.41421i)
12Conjugate of (2.45+1.41421i) = (2.45-1.41421i)
13(2.45+5.33i) ^ (2.45+1.41421i) = (8.37072-12.7078i)


## String

 1// Concatenation
2string tt = "toto1" + 1 + " -- 77";
3
4// Append
5string t1 = "0123456789";
6t1(4:3) = "abcdefghijk-";
7
8// Sub string
9string t55 = t1(4:14);
10
11cout << "tt = " << tt << endl;
12
13cout << "t1 = " << t1 << endl;
14cout << "t1.find(abc) = " << t1.find("abc") << endl;
15cout << "t1.rfind(abc) = " << t1.rfind("abc") << endl;
16cout << "t1.find(abc, 10) = " << t1.find("abc",10) << endl;
17cout << "t1.ffind(abc, 10) = " << t1.rfind("abc",10) << endl;
18cout << "t1.length = " << t1.length << endl;
19
20cout << "t55 = " << t55 << endl;


The output of this script is:

1tt = toto11 -- 77
2t1 = 0123abcdefghijk-456789
3t1.find(abc) = 4
4t1.rfind(abc) = 4
5t1.find(abc, 10) = -1
6t1.ffind(abc, 10) = 4
7t1.length = 22
8t55 = abcdefghijk


## Elementary function

 1real b = 1.;
2real a = b;
3func real phix(real t){
4    return (a+b)*cos(t) - b*cos(t*(a+b)/b);
5}
6func real phiy(real t){
7    return (a+b)*sin(t) - b*sin(t*(a+b)/b);
8}
9
10border C(t=0, 2*pi){x=phix(t); y=phiy(t);}
11mesh Th = buildmesh(C(50));
12plot(Th);


## Array

  1real[int] tab(10), tab1(10); //2 array of 10 real
2//real[int] tab2; //bug: array with no size
3
4tab = 1.03; //set all the array to 1.03
5tab[1] = 2.15;
6
7cout << "tab: " << tab << endl;
8cout << "min: " << tab.min << endl;
9cout << "max: " << tab.max << endl;
10cout << "sum: " << tab.sum << endl;
11
12tab.resize(12); //change the size of array tab to 12 with preserving first value
13tab(10:11) = 3.14; //set values 10 & 11
14cout << "resized tab: " << tab << endl;
15
16tab.sort ; //sort the array tab
17cout << "sorted tab:" << tab << endl;
18
19real[string] tt; //array with string index
20tt["+"] = 1.5;
21cout << "tt[\"a\"] = " << tt["a"] << endl;
22cout << "tt[\"+\"] = " << tt["+"] << endl;
23
24real[int] a(5), b(5), c(5), d(5);
25a = 1;
26b = 2;
27c = 3;
28a[2] = 0;
29d = ( a ? b : c ); //for i = 0, n-1 : d[i] = a[i] ? b[i] : c[i]
30cout << " d = ( a ? b : c ) is " << d << endl;
31d = ( a ? 1 : c ); //for i = 0, n-1: d[i] = a[i] ? 1 : c[i]
32d = ( a ? b : 0 ); //for i = 0, n-1: d[i] = a[i] ? b[i] : 0
33d = ( a ? 1 : 0 ); //for i = 0, n-1: d[i] = a[i] ? 0 : 1
34
35int[int] ii(0:d.n-1); //set array ii to 0, 1, ..., d.n-1
36d = -1:-5; //set d to -1, -2, ..., -5
37
38sort(d, ii); //sort array d and ii in parallel
39cout << "d: " << d << endl;
40cout << "ii: " << ii << endl;
41
42
43{
44    int[int] A1(2:10); //2, 3, 4, 5, 6, 7, 8, 9, 10
45    int[int] A2(2:3:10); //2, 5, 8
46    cout << "A1(2:10): " << A1 << endl;
47    cout << "A2(2:3:10): " << A1 << endl;
48    A1 = 1:2:5;
49    cout << "1:2:5 => " << A1 << endl;
50}
51{
52    real[int] A1(2:10); //2, 3, 4, 5, 6, 7, 8, 9, 10
53    real[int] A2(2:3:10); //2, 5, 8
54    cout << "A1(2:10): " << A1 << endl;
55    cout << "A2(2:3:10): " << A1 << endl;
56    A1 = 1.:0.5:3.999;
57    cout << "1.:0.5:3.999 => " << A1 << endl;
58}
59{
60    complex[int] A1(2.+0i:10.+0i); //2, 3, 4, 5, 6, 7, 8, 9, 10
61    complex[int] A2(2.:3.:10.); //2, 5, 8
62    cout << " A1(2.+0i:10.+0i): " << A1 << endl;
63    cout << " A2(2.:3.:10.)= " << A2 << endl;
64    cout << " A1.re real part array: " << A1.re << endl ;
65    // he real part array of the complex array
66    cout << " A1.im imag part array: " << A1.im << endl ;
67    //the imaginary part array of the complex array
68}
69
70// Integer array operators
71{
72    int N = 5;
73    real[int] a(N), b(N), c(N);
74    a = 1;
75    a(0:4:2) = 2;
76    a(3:4) = 4;
77    cout << "a: " << a << endl;
78    b = a + a;
79    cout <<"b = a + a: " << b << endl;
80    b += a;
81    cout <<"b += a: " << b << endl;
82    b += 2*a;
83    cout <<"b += 2*a: " << b << endl;
84    b /= 2;
85    cout <<" b /= 2: " << b << endl;
86    b .*= a; // same as b = b .* a
87    cout << "b .*= a: " << b << endl;
88    b ./= a; //same as b = b ./ a
89    cout << "b ./= a: " << b << endl;
90    c = a + b;
91    cout << "c = a + b: " << c << endl;
92    c = 2*a + 4*b;
93    cout << "c = 2*a + 4b: " << c << endl;
94    c = a + 4*b;
95    cout << "c = a + 4b: " << c << endl;
96    c = -a + 4*b;
97    cout << "c = -a + 4b: " << c << endl;
98    c = -a - 4*b;
99    cout << "c = -a - 4b: " << c << endl;
100    c = -a - b;
101    cout << "c = -a -b: " << c << endl;
102
103    c = a .* b;
104    cout << "c = a .* b: " << c << endl;
105    c = a ./ b;
106    cout << "c = a ./ b: " << c << endl;
107    c = 2 * b;
108    cout << "c = 2 * b: " << c << endl;
109    c = b * 2;
110    cout << "c = b * 2: " << c << endl;
111
112    //this operator do not exist
113    //c = b/2;
114    //cout << "c = b / 2: " << c << endl;
115
116    //Array methods
117    cout << "||a||_1 = " << a.l1 << endl;
118    cout << "||a||_2 = " << a.l2 << endl;
119    cout << "||a||_infty = " << a.linfty << endl;
120    cout << "sum a_i = " << a.sum << endl;
121    cout << "max a_i = " << a.max << " a[ " << a.imax << " ] = " << a[a.imax] << endl;
122    cout << "min a_i = " << a.min << " a[ " << a.imin << " ] = " << a[a.imin] << endl;
123
124    cout << "a' * a = " << (a'*a) << endl;
125    cout << "a quantile 0.2 = " << a.quantile(0.2) << endl;
126
127    //Array mapping
128    int[int] I = [2, 3, 4, -1, 3];
129    b = c = -3;
130    b = a(I); //for (i = 0; i < b.n; i++) if (I[i] >= 0) b[i] = a[I[i]];
131    c(I) = a; //for (i = 0; i < I.n; i++) if (I[i] >= 0) C(I[i]) = a[i];
132    cout << "b = a(I) : " << b << endl;
133    cout << "c(I) = a " << c << endl;
134    c(I) += a; //for (i = 0; i < I.n; i++) if (I[i] >= 0) C(I[i]) += a[i];
135    cout << "b = a(I) : " << b << endl;
136    cout << "c(I) = a " << c << endl;
137
138}
139
140{
141    // Array versus matrix
142    int N = 3, M = 4;
143
144    real[int, int] A(N, M);
145    real[int] b(N), c(M);
146    b = [1, 2, 3];
147    c = [4, 5, 6, 7];
148
149    complex[int, int] C(N, M);
150    complex[int] cb = [1, 2, 3], cc = [10i, 20i, 30i, 40i];
151
152    b = [1, 2, 3];
153
154    int [int] I = [2, 0, 1];
155    int [int] J = [2, 0, 1, 3];
156
157    A = 1; //set all the matrix
158    A(2, :) = 4; //the full line 2
159    A(:, 1) = 5; //the full column 1
160    A(0:N-1, 2) = 2; //set the column 2
161    A(1, 0:2) = 3; //set the line 1 from 0 to 2
162
163    cout << "A = " << A << endl;
164
165    //outer product
166    C = cb * cc';
167    C += 3 * cb * cc';
168    C -= 5i * cb * cc';
169    cout << "C = " << C << endl;
170
171    //this transforms an array into a sparse matrix
172    matrix B;
173    B = A;
174    B = A(I, J); //B(i, j) = A(I(i), J(j))
175    B = A(I^-1, J^-1); //B(I(i), J(j)) = A(i,j)
176
177    //outer product
178    A = 2. * b * c';
179    cout << "A = " << A << endl;
180    B = b*c'; //outer product B(i, j) = b(i)*c(j)
181    B = b*c'; //outer product B(i, j) = b(i)*c(j)
182    B = (2*b*c')(I, J); //outer product B(i, j) = b(I(i))*c(J(j))
183    B = (3.*b*c')(I^-1,J^-1); //outer product B(I(i), J(j)) = b(i)*c(j)
184    cout << "B = (3.*b*c')(I^-1,J^-1) = " << B << endl;
185
186    //row and column of the maximal coefficient of A
187    int i, j, ii, jj;
188    ijmax(A, ii, jj);
189
190    i = A.imax;
191    j = A.jmax;
192
193    cout << "Max " << i << " " << j << ", = " << A.max << endl;
194
195    //row and column of the minimal coefficient of A
196    ijmin(A, i, j);
197
198    ii = A.imin;
199    jj = A.jmin;
200
201    cout << "Min " << ii << " " << jj << ", = " << A.min << endl;
202}


The output os this script is:

  1tab: 10
2    1.03    2.15    1.03    1.03    1.03
3    1.03    1.03    1.03    1.03    1.03
4
5min: 1.03
6max: 2.15
7sum: 11.42
8resized tab: 12
9    1.03    2.15    1.03    1.03    1.03
10    1.03    1.03    1.03    1.03    1.03
11    3.14    3.14
12sorted tab:12
13    1.03    1.03    1.03    1.03    1.03
14    1.03    1.03    1.03    1.03    2.15
15    3.14    3.14
16tt["a"] = 0
17tt["+"] = 1.5
18 d = ( a ? b : c ) is 5
19      2   2   3   2   2
20
21d: 5
22     -5  -4  -3  -2  -1
23
24ii: 5
25      4   3   2   1   0
26
27A1(2:10): 9
28      2   3   4   5   6
29      7   8   9  10
30A2(2:3:10): 9
31      2   3   4   5   6
32      7   8   9  10
331:2:5 => 3
34      1   3   5
35A1(2:10): 9
36      2   3   4   5   6
37      7   8   9  10
38A2(2:3:10): 9
39      2   3   4   5   6
40      7   8   9  10
411.:0.5:3.999 => 6
42      1 1.5   2 2.5   3
43    3.5
44 A1(2.+0i:10.+0i): 9
45    (2,0)   (3,0)   (4,0)   (5,0)   (6,0)
46    (7,0)   (8,0)   (9,0)   (10,0)
47 A2(2.:3.:10.)= 3
48    (2,0)   (5,0)   (8,0)
49 A1.re real part array: 9
50      2   3   4   5   6
51      7   8   9  10
52 A1.im imag part array: 9
53      0   0   0   0   0
54      0   0   0   0
55a: 5
56      2   1   2   4   4
57
58b = a + a: 5
59      4   2   4   8   8
60
61b += a: 5
62      6   3   6  12  12
63
64b += 2*a: 5
65     10   5  10  20  20
66
67 b /= 2: 5
68      5 2.5   5  10  10
69
70b .*= a: 5
71     10 2.5  10  40  40
72
73b ./= a: 5
74      5 2.5   5  10  10
75
76c = a + b: 5
77      7 3.5   7  14  14
78
79c = 2*a + 4b: 5
80     24  12  24  48  48
81
82c = a + 4b: 5
83     22  11  22  44  44
84
85c = -a + 4b: 5
86     18   9  18  36  36
87
88c = -a - 4b: 5
89    -22 -11 -22 -44 -44
90
91c = -a -b: 5
92     -7 -3.5     -7 -14 -14
93
94c = a .* b: 5
95     10 2.5  10  40  40
96
97c = a ./ b: 5
98    0.4 0.4 0.4 0.4 0.4
99
100c = 2 * b: 5
101     10   5  10  20  20
102
103c = b * 2: 5
104     10   5  10  20  20
105
106||a||_1 = 13
107||a||_2 = 6.40312
108||a||_infty = 4
109sum a_i = 13
110max a_i = 4 a[ 3 ] = 4
111min a_i = 1 a[ 1 ] = 1
112a' * a = 41
113a quantile 0.2 = 2
114b = a(I) : 5
115      2   4   4  -3   4
116
117c(I) = a 5
118     -3  -3   2   4   2
119
120b = a(I) : 5
121      2   4   4  -3   4
122
123c(I) = a 5
124     -3  -3   4   9   4
125
126A = 3 4
127       1   5   2   1
128       3   3   3   1
129       4   5   2   4
130
131C = 3 4
132     (-50,-40) (-100,-80) (-150,-120) (-200,-160)
133     (-100,-80) (-200,-160) (-300,-240) (-400,-320)
134     (-150,-120) (-300,-240) (-450,-360) (-600,-480)
135
136A = 3 4
137       8  10  12  14
138      16  20  24  28
139      24  30  36  42
140
141B = (3.*b*c')(I^-1,J^-1) = # Sparse Matrix (Morse)
142# first line: n m (is symmetic) nbcoef
143# after for each nonzero coefficient:   i j a_ij where (i,j) \in  {1,...,n}x{1,...,m}
1443 4 0  12
145        1         1 10
146        1         2 12
147        1         3 8
148        1         4 14
149        2         1 15
150        2         2 18
151        2         3 12
152        2         4 21
153        3         1 5
154        3         2 6
155        3         3 4
156        3         4 7


## Block matrix

 1// Parameters
2real f1 = 1.;
3real f2 = 1.5;
4
5// Mesh
6mesh Th1 = square(10, 10);
7mesh Th2 = square(10, 10, [1+x, -1+y]);
8plot(Th1, Th2);
9
10// Fespace
11fespace Uh1(Th1, P1);
12Uh1 u1;
13
14fespace Uh2(Th2, P2);
15Uh2 u2;
16
17// Macro
18macro grad(u) [dx(u), dy(u)] //
19
20// Problem
21varf vPoisson1 (u, v)
22    = int2d(Th1)(
23          grad(u)' * grad(v)
24    )
25    - int2d(Th1)(
26          f1 * v
27    )
28    + on(1, 2, 3, 4, u=0)
29    ;
30
31varf vPoisson2 (u, v)
32    = int2d(Th2)(
33          grad(u)' * grad(v)
34    )
35    - int2d(Th2)(
36          f1 * v
37    )
38    + on(1, 2, 3, 4, u=0)
39    ;
40matrix<real> Poisson1 = vPoisson1(Uh1, Uh1);
41real[int] Poisson1b = vPoisson1(0, Uh1);
42
43matrix<real> Poisson2 = vPoisson2(Uh2, Uh2);
44real[int] Poisson2b = vPoisson2(0, Uh2);
45
46//block matrix
47matrix<real> G = [[Poisson1, 0], [0, Poisson2]];
48set(G, solver=sparsesolver);
49
50//block right hand side
51real[int] Gb = [Poisson1b, Poisson2b];
52
53// Solve
54real[int] sol = G^-1 * Gb;
55
56// Dispatch
57[u1[], u2[]] = sol;
58
59// Plot
60plot(u1, u2);


## Matrix operations

 1// Mesh
2mesh Th = square(2, 1);
3
4// Fespace
5fespace Vh(Th, P1);
6Vh f, g;
7f = x*y;
8g = sin(pi*x);
9
10Vh<complex> ff, gg; //a complex valued finite element function
11ff= x*(y+1i);
12gg = exp(pi*x*1i);
13
14// Problem
15varf mat (u, v)
16    = int2d(Th)(
17          1*dx(u)*dx(v)
18        + 2*dx(u)*dy(v)
19        + 3*dy(u)*dx(v)
20        + 4*dy(u)*dy(v)
21    )
22    + on(1, 2, 3, 4, u=1)
23    ;
24
25varf mati (u, v)
26    = int2d(Th)(
27         1*dx(u)*dx(v)
28        + 2i*dx(u)*dy(v)
29        + 3*dy(u)*dx(v)
30        + 4*dy(u)*dy(v)
31    )
32    + on(1, 2, 3, 4, u=1)
33    ;
34
35matrix A = mat(Vh, Vh);
36matrix<complex> AA = mati(Vh, Vh); //a complex sparse matrix
37
38// Operations
39Vh m0; m0[] = A*f[];
40Vh m01; m01[] = A'*f[];
41Vh m1; m1[] = f[].*g[];
42Vh m2; m2[] = f[]./g[];
43
44// Display
45cout << "f = " << f[] << endl;
46cout << "g = " << g[] << endl;
47cout << "A = " << A << endl;
48cout << "m0 = " << m0[] << endl;
49cout << "m01 = " << m01[] << endl;
50cout << "m1 = "<< m1[] << endl;
51cout << "m2 = "<< m2[] << endl;
52cout << "dot Product = "<< f[]'*g[] << endl;
53cout << "hermitien Product = "<< ff[]'*gg[] << endl;
54cout << "outer Product = "<< (A=f[]*g[]') << endl;
55cout << "hermitien outer Product = "<< (AA=ff[]*gg[]') << endl;
56
57// Diagonal
58real[int] diagofA(A.n);
59diagofA = A.diag; //get the diagonal of the matrix
60A.diag = diagofA ; //set the diagonal of the matrix
61
62// Sparse matrix set
63int[int] I(1), J(1);
64real[int] C(1);
65
66[I, J, C] = A; //get the sparse term of the matrix A (the array are resized)
67cout << "I = " << I << endl;
68cout << "J = " << J << endl;
69cout << "C = " << C << endl;
70
71A = [I, J, C]; //set a new matrix
72matrix D = [diagofA]; //set a diagonal matrix D from the array diagofA
73cout << "D = " << D << endl;


The output of this script is:

  1f = 6
2      0   0   0   0 0.5
3      1
4g = 6
5      0   1 1.224646799e-16   0   1
6    1.224646799e-16
7A = # Sparse Matrix (Morse)
8# first line: n m (is symmetic) nbcoef
9# after for each nonzero coefficient:   i j a_ij where (i,j) \in  {1,...,n}x{1,...,m}
106 6 0  24
11        1         1 1.0000000000000000199e+30
12        1         2 0.49999999999999994449
13        1         4 0
14        1         5 -2.5
15        2         1 0
16        2         2 1.0000000000000000199e+30
17        2         3 0.49999999999999994449
18        2         5 0.49999999999999977796
19        2         6 -2.5
20        3         2 0
21        3         3 1.0000000000000000199e+30
22        3         6 0.49999999999999977796
23        4         1 0.49999999999999977796
24        4         4 1.0000000000000000199e+30
25        4         5 0
26        5         1 -2.5
27        5         2 0.49999999999999977796
28        5         4 0.49999999999999994449
29        5         5 1.0000000000000000199e+30
30        5         6 0
31        6         2 -2.5
32        6         3 0
33        6         5 0.49999999999999994449
34        6         6 1.0000000000000000199e+30
35
36m0 = 6
37    -1.25   -2.25   0.5   0 5e+29
38    1e+30
39m01 = 6
40    -1.25   -2.25     0 0.25    5e+29
41    1e+30
42m1 = 6
43      0   0   0   0 0.5
44    1.224646799e-16
45m2 = 6
46    -nan      0   0 -nan    0.5
47    8.165619677e+15
48dot Product = 0.5
49hermitien Product = (1.11022e-16,2.5)
50outer Product = # Sparse Matrix (Morse)
51# first line: n m (is symmetic) nbcoef
52# after for each nonzero coefficient:   i j a_ij where (i,j) \in  {1,...,n}x{1,...,m}
536 6 0  8
54        5         2 0.5
55        5         3 6.1232339957367660359e-17
56        5         5 0.5
57        5         6 6.1232339957367660359e-17
58        6         2 1
59        6         3 1.2246467991473532072e-16
60        6         5 1
61        6         6 1.2246467991473532072e-16
62
63hermitien outer Product = # Sparse Matrix (Morse)
64# first line: n m (is symmetic) nbcoef
65# after for each nonzero coefficient:   i j a_ij where (i,j) \in  {1,...,n}x{1,...,m}
666 6 0  24
67        2         1 (0,0.5)
68        2         2 (0.5,3.0616169978683830179e-17)
69        2         3 (6.1232339957367660359e-17,-0.5)
70        2         4 (0,0.5)
71        2         5 (0.5,3.0616169978683830179e-17)
72        2         6 (6.1232339957367660359e-17,-0.5)
73        3         1 (0,1)
74        3         2 (1,6.1232339957367660359e-17)
75        3         3 (1.2246467991473532072e-16,-1)
76        3         4 (0,1)
77        3         5 (1,6.1232339957367660359e-17)
78        3         6 (1.2246467991473532072e-16,-1)
79        5         1 (0.5,0.5)
80        5         2 (0.5,-0.49999999999999994449)
81        5         3 (-0.49999999999999994449,-0.50000000000000011102)
82        5         4 (0.5,0.5)
83        5         5 (0.5,-0.49999999999999994449)
84        5         6 (-0.49999999999999994449,-0.50000000000000011102)
85        6         1 (1,1)
86        6         2 (1,-0.99999999999999988898)
87        6         3 (-0.99999999999999988898,-1.000000000000000222)
88        6         4 (1,1)
89        6         5 (1,-0.99999999999999988898)
90        6         6 (-0.99999999999999988898,-1.000000000000000222)
91
92I = 8
93      4   4   4   4   5
94      5   5   5
95J = 8
96      1   2   4   5   1
97      2   4   5
98C = 8
99    0.5 6.123233996e-17 0.5 6.123233996e-17   1
100    1.224646799e-16   1 1.224646799e-16
101  -- Raw Matrix    nxm  =6x6 nb  none zero coef. 8
102  -- Raw Matrix    nxm  =6x6 nb  none zero coef. 6
103D = # Sparse Matrix (Morse)
104# first line: n m (is symmetic) nbcoef
105# after for each nonzero coefficient:   i j a_ij where (i,j) \in  {1,...,n}x{1,...,m}
1066 6 1  6
107        1         1 0
108        2         2 0
109        3         3 0
110        4         4 0
111        5         5 0.5
112        6         6 1.2246467991473532072e-16


Warning

Due to Fortran indices starting at one, the output of a diagonal matrix D is indexed from 1. but in FreeFEM, the indices start from 0.

## Matrix inversion

 1load "lapack"
2load "fflapack"
3
4// Matrix
5int n = 5;
6real[int, int] A(n, n), A1(n, n), B(n,n);
7for (int i = 0; i < n; ++i)
8    for (int j = 0; j < n; ++j)
9        A(i, j) = (i == j) ? n+1 : 1;
10cout << A << endl;
11
12// Inversion (lapack)
13A1 = A^-1; //def in "lapack"
14cout << A1 << endl;
15
16B = 0;
17for (int i = 0; i < n; ++i)
18    for (int j = 0; j < n; ++j)
19        for (int k = 0; k < n; ++k)
20            B(i, j) += A(i,k)*A1(k,j);
21cout << B << endl;
22
23// Inversion (fflapack)
24inv(A1); //def in "fflapack"
25cout << A1 << endl;


The output of this script is:

 15 5
2       6   1   1   1   1
3       1   6   1   1   1
4       1   1   6   1   1
5       1   1   1   6   1
6       1   1   1   1   6
7
85 5
9     0.18 -0.02 -0.02 -0.02 -0.02
10     -0.02 0.18 -0.02 -0.02 -0.02
11     -0.02 -0.02 0.18 -0.02 -0.02
12     -0.02 -0.02 -0.02 0.18 -0.02
13     -0.02 -0.02 -0.02 -0.02 0.18
14
155 5
16       1 1.040834086e-17 1.040834086e-17 1.734723476e-17 2.775557562e-17
17     3.469446952e-18   1 -1.734723476e-17 1.734723476e-17 2.775557562e-17
18     2.428612866e-17 -3.122502257e-17   1 1.734723476e-17 2.775557562e-17
19     2.081668171e-17 -6.938893904e-17 -3.469446952e-17   1   0
20     2.775557562e-17 -4.163336342e-17 -2.775557562e-17   0   1
21
225 5
23       6   1   1   1   1
24       1   6   1   1   1
25       1   1   6   1   1
26       1   1   1   6   1
27       1   1   1   1   6


Tip

To compile lapack.cpp and fflapack.cpp, you must have the lapack library on your system and compile the plugin with the command:

1ff-c++ lapack.cpp -llapack     ff-c++ fflapack.cpp -llapack


## FE array

 1// Mesh
2mesh Th = square(20, 20, [2*x, 2*y]);
3
4// Fespace
5fespace Vh(Th, P1);
6Vh u, v, f;
7
8// Problem
9problem Poisson (u, v)
10    = int2d(Th)(
11          dx(u)*dx(v)
12        + dy(u)*dy(v)
13    )
14    + int2d(Th)(
15        - f*v
16    )
17    + on(1, 2, 3, 4, u=0)
18    ;
19
20Vh[int] uu(3); //an array of FE function
21// Solve problem 1
22f = 1;
23Poisson;
24uu[0] = u;
25// Solve problem 2
26f = sin(pi*x)*cos(pi*y);
27Poisson;
28uu[1] = u;
29// Solve problem 3
30f = abs(x-1)*abs(y-1);
31Poisson;
32uu[2] = u;
33
34// Plot
35for (int i = 0; i < 3; i++)
36    plot(uu[i], wait=true);


Finite element array

## Loop

 1for (int i = 0; i < 10; i=i+1)
2    cout << i << endl;
3
4real eps = 1.;
5while (eps > 1e-5){
6    eps = eps/2;
7    if (i++ < 100)
8        break;
9    cout << eps << endl;
10}
11
12for (int j = 0; j < 20; j++){
13    if (j < 10) continue;
14    cout << "j = " << j << endl;
15}


## Implicit loop

 1real [int, int] a(10, 10);
2real [int] b(10);
3
4for [i, bi : b]{
5    bi = i+1;
6    cout << i << " " << bi << endl;
7}
8cout << "b = " << b << endl;
9
10for [i, j, aij : a]{
11    aij = 1./(2+i+j);
12    if (abs(aij) < 0.2) aij = 0;
13}
14cout << "a = " << a << endl;
15
16matrix A = a;
17string[string] ss; //a map
18ss["1"] = 1;
19ss["2"] = 2;
20ss["3"] = 5;
21for [i, bi : ss]
22    bi = i + 6 + "-dddd";
23cout << "ss = " << ss << endl;
24
25int[string] si;
26si[1] = 2;
27si[50] = 1;
28for [i, vi : si]{
29    cout << " i " << setw(3) << i << " " << setw(10) << vi << endl;
30    vi = atoi(i)*2;
31}
32cout << "si = " << si << endl;
33
34for [i, j, aij : A]{
35    cout << i << " " << j << " " << aij << endl;
36    aij = -aij;
37}
38cout << A << endl;


The output of this script is:

 10 1
21 2
32 3
43 4
54 5
65 6
76 7
87 8
98 9
109 10
11b = 10
12      1   2   3   4   5
13      6   7   8   9  10
14
15a = 10 10
16     0.5 0.3333333333 0.25 0.2   0   0   0   0   0   0
17     0.3333333333 0.25 0.2   0   0   0   0   0   0   0
18     0.25 0.2   0   0   0   0   0   0   0   0
19     0.2   0   0   0   0   0   0   0   0   0
20       0   0   0   0   0   0   0   0   0   0
21       0   0   0   0   0   0   0   0   0   0
22       0   0   0   0   0   0   0   0   0   0
23       0   0   0   0   0   0   0   0   0   0
24       0   0   0   0   0   0   0   0   0   0
25       0   0   0   0   0   0   0   0   0   0
26
27ss = 1 1
282 2
293 5
30
31 i   1          2
32 i  50          1
33si = 1 2
3450 100
35
360 0 0.5
370 1 0.333333
380 2 0.25
390 3 0.2
401 0 0.333333
411 1 0.25
421 2 0.2
432 0 0.25
442 1 0.2
453 0 0.2
46# Sparse Matrix (Morse)
47# first line: n m (is symmetic) nbcoef
48# after for each nonzero coefficient:   i j a_ij where (i,j) \in  {1,...,n}x{1,...,m}
4910 10 0  10
50        1         1 -0.5
51        1         2 -0.33333333333333331483
52        1         3 -0.25
53        1         4 -0.2000000000000000111
54        2         1 -0.33333333333333331483
55        2         2 -0.25
56        2         3 -0.2000000000000000111
57        3         1 -0.25
58        3         2 -0.2000000000000000111
59        4         1 -0.2000000000000000111


## I/O

 1int i;
2cout << "std-out" << endl;
3cout << " enter i = ?";
4cin >> i;
5
6{
7    ofstream f("toto.txt");
8    f << i << "hello world'\n";
9} //close the file f because the variable f is delete
10
11{
12    ifstream f("toto.txt");
13    f >> i;
14}
15
16{
17    ofstream f("toto.txt", append);
18    //to append to the existing file "toto.txt"
19    f << i << "hello world'\n";
20} //close the file f because the variable f is delete
21
22cout << i << endl;


## File stream

 1int where;
2real[int] f = [0, 1, 2, 3, 4, 5];
3real[int] g(6);
4
5{
6    ofstream file("f.txt", binary);
7    file.precision(16);
8    file << f << endl;
9    where = file.tellp();
10    file << 0.1 ;
11
12    cout << "Where in file " << where << endl;
13    file << " # comment bla bla ... 0.3 \n";
14    file << 0.2 << endl;
15    file.flush; //to flush the buffer of file
16}
17
18//Function to skip comment starting with # in a file
19func ifstream skipcomment(ifstream &ff){
20    while(1){
21        int where = ff.tellg(); //store file position
22        string comment;
23        ff >> comment;
24        if (!ff.good()) break;
25        if (comment(0:0)=="#"){
26            getline(ff, comment);
27            cout << " -- #" << comment << endl;
28        }
29        else{
30            ff.seekg(where); //restore file position
31            break;
32        }
33    }
34    return ff;
35}
36
37{
38    real xx;
39    ifstream file("f.txt", binary);
40    cout << "Where " << file.seekg << endl;
41    file.seekg(where);
42    file >> xx;
43    cout << " xx = " << xx << " good ? " << file.good() << endl;
44    assert(xx == 0.1);
45    skipcomment(file) >> xx;
46    assert(xx == 0.2);
47    file.seekg(0); //rewind
48    cout << "Where " << file.tellg() << " " << file.good() << endl;
49    file >> g;
50}


## Command line arguments

When using the command:

1FreeFem++ script.edp arg1 arg2


The arguments can be used in the script with:

1for (int i = 0; i < ARGV.n; i++)
2    cout << ARGV[i] << endl;


When using the command:

1FreeFem++ script.edp -n 10 -a 1. -d 42.


The arguments can be used in the script with:

1include "getARGV.idp"
2
3int n = getARGV("-n", 1);
4real a = getARGV("-a", 1.);
5real d = getARGV("-d", 1.);


## Macro

 1// Macro without parameters
2macro xxx() {
3    real i = 0;
4    int j = 0;
5    cout << i << " " << j << endl;
6}//
7
8xxx
9
10// Macro with parameters
11macro toto(i) i //
12
13toto({real i = 0; int j = 0; cout << i << " " << j << endl;})
14
15// Macro as parameter of a macro
16real[int,int] CC(7, 7), EE(6, 3), EEps(4, 4);
17
18macro VIL6(v, i) [v(1,i), v(2,i), v(4,i), v(5,i), v(6,i)] //
19macro VIL3(v, i) [v(1,i), v(2,i)] //
20macro VV6(v, vv) [
21    v(vv,1), v(vv,2),
22    v(vv,4), v(vv,5),
23    v(vv,6)] //
24macro VV3(v, vv) [v(vv,1), v(vv,2)] //
25
26func C5x5 = VV6(VIL6, CC);
27func E5x2 = VV6(VIL3, EE);
28func Eps = VV3(VIL3, EEps);
29
30// Macro concatenation
31mesh Th = square(2, 2);
32fespace Vh(Th, P1);
33Vh Ux=x, Uy=y;
34
35macro div(V) (dx(V#x) + dy(V#y)) //
36
37cout << int2d(Th)(div(U)) << endl;
38
39// Verify the quoting
40macro foo(i, j, k) i j k //
41foo(, , )
42foo({int[}, {int] a(10}, {);})
43
44//NewMacro - EndMacro
45NewMacro grad(u) [dx(u), dy(u)] EndMacro
46cout << int2d(Th)(grad(Ux)' * grad(Uy)) << endl;
47
48// IFMACRO - ENDIFMACRO
49macro AA CAS1 //
50
51IFMACRO(AA,CAS1 )
52cout << "AA = " << Stringification(AA) << endl;
53macro CASE file1.edp//
54ENDIFMACRO
55IFMACRO(AA, CAS2)
56macro CASE file2.edp//
57ENDIFMACRO
58
59cout << "CASE = " << Stringification(CASE) << endl;
60
61IFMACRO(CASE)
62include Stringification(CASE)
63ENDIFMACRO
64
65// FILE - LINE
66cout << "In " << FILE << ", line " << LINE << endl;


The output script generated with macros is:

 11 : // Macro without parameters
22 :  macro xxx {
33 :     real i = 0;
44 :     int j = 0;
55 :     cout << i << " " << j << endl;
66 : }//
77 :
88 :
91 :
102 :
113 :
124 :  {
131 :     real i = 0;
142 :     int j = 0;
153 :     cout << i << " " << j << endl;
164 : }
179 :
1810 : // Macro with parameters
1911 :  macro toto(i )   i //
2012 :
2113 :                    real i = 0; int j = 0; cout << i << " " << j << endl;
2214 :
2315 : // Macro as parameter of a macro
2416 : real[int,int] CC(7, 7), EE(6, 3), EEps(4, 4);
2517 :
2618 :   macro VIL6(v,i )   [v(1,i), v(2,i), v(4,i), v(5,i), v(6,i)] //
2719 :   macro VIL3(v,i )   [v(1,i), v(2,i)] //
2820 :   macro VV6(v,vv )   [
2921 :    v(vv,1), v(vv,2),
3022 :    v(vv,4), v(vv,5),
3123 :    v(vv,6)] //
3224 :   macro VV3(v,vv )   [v(vv,1), v(vv,2)] //
3325 :
3426 : func C5x5 =
351 :
362 :
373 :       [
381 :             [ CC(1,1),  CC(2,1),  CC(4,1),  CC(5,1),  CC(6,1)] ,         [ CC(1,2),  CC(2,2),  CC(4,2),  CC(5,2),  CC(6,2)] ,
392 :             [ CC(1,4),  CC(2,4),  CC(4,4),  CC(5,4),  CC(6,4)] ,         [ CC(1,5),  CC(2,5),  CC(4,5),  CC(5,5),  CC(6,5)] ,
403 :             [ CC(1,6),  CC(2,6),  CC(4,6),  CC(5,6),  CC(6,6)] ] ;
4127 : func E5x2 =
421 :
432 :
443 :       [
451 :          [ EE(1,1),  EE(2,1)] ,      [ EE(1,2),  EE(2,2)] ,
462 :          [ EE(1,4),  EE(2,4)] ,      [ EE(1,5),  EE(2,5)] ,
473 :          [ EE(1,6),  EE(2,6)] ] ;
4828 : func Eps =      [     [ EEps(1,1),  EEps(2,1)] ,      [ EEps(1,2),  EEps(2,2)] ] ;
4929 :
5030 : // Macro concatenation
5131 : mesh Th = square(2, 2);
5232 : fespace Vh(Th, P1);
5333 : Vh Ux=x, Uy=y;
5434 :
5535 :  macro div(V )   (dx(V#x) + dy(V#y)) //
5636 :
5737 : cout << int2d(Th)(     (dx(Ux) + dy(Uy)) ) << endl;
5838 :
5939 : // Verify the quoting
6040 :    macro foo(i,j,k )   i j k //
6141 :
6242 :         int[ int] a(10 );
6343 :
6444 : //NewMacro - EndMacro
6545 :  macro grad(u )   [dx(u), dy(u)]
6646 : cout << int2d(Th)(    [dx(Ux), dy(Ux)] ' *     [dx(Uy), dy(Uy)] ) << endl;
6747 :
6848 : // IFMACRO - ENDIFMACRO
6949 :   macro AACAS1 //
7050 :
7151 :
721 : cout << "AA = " << Stringification( CAS1 ) << endl;
732 :   macro CASEfile1.edp//
743 :
7552 :
7653 :
7754 : cout << "CASE = " << Stringification(file1.edp) << endl;
7855 :
7956 :
801 : include Stringification(file1.edp)cout << "This is the file 1" << endl;
812 :
822 :
8357 :
8458 : // FILE - LINE
8559 : cout << "In " << FILE << ", line " << LINE << endl;


The output os this script is:

1AA = CAS1
2CASE = file1.edp
3This is the file 1
4In Macro.edp, line 59


## Basic error handling

1real a;
2try{
3    a = 1./0.;
4}
5catch (...) //all exceptions can be caught
6{
7    cout << "Catch an ExecError" << endl;
8    a = 0.;
9}


The output of this script is:

11/0 : d d d
2  current line = 3
3Exec error :  Div by 0
4   -- number :1
5Catch an ExecError


## Error handling

 1// Parameters
2int nn = 5;
3func f = 1; //right hand side function
4func g = 0; //boundary condition function
5
6// Mesh
7mesh Th = square(nn, nn);
8
9// Fespace
10fespace Vh(Th, P1);
11Vh uh, vh;
12
13// Problem
14real cpu = clock();
15problem laplace (uh, vh, solver=Cholesky, tolpivot=1e-6)
16    = int2d(Th)(
17          dx(uh)*dx(vh)
18        + dy(uh)*dy(vh)
19    )
20    + int2d(Th)(
21        - f*vh
22    )
23    ;
24
25try{
26    cout << "Try Cholesky" << endl;
27
28    // Solve
29    laplace;
30
31    // Plot
32    plot(uh);
33
34    // Display
35    cout << "laplacian Cholesky " << nn << ", x_" << nn << " : " << -cpu+clock() << " s, max = " << uh[].max << endl;
36}
37catch(...) { //catch all error
38    cout << " Catch cholesky PB " << endl;
39}


The output of this script is:

1Try Cholesky
2ERREUR choleskypivot (35)= -6.43929e-15 < 1e-06
3  current line = 29
4Exec error : FATAL ERREUR dans ./../femlib/MatriceCreuse_tpl.hpp
5cholesky line:
6   -- number :688
7 catch an erreur in  solve  =>  set  sol = 0 !!!!!!!
8 Catch cholesky PB