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);
../_images/ElementaryFunction.png

Fig. 223 Mesh

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);
../_images/BlockMatrix.png

Fig. 224 Result

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);
FEArray1

Fig. 225 First result

FEArray2

Fig. 226 Second result

FEArray3

Fig. 227 Third result

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
Table of content