FreeFEM Documentation on GitHub

stars - forks

Developers

FFT

 1 load "dfft"
 2 
 3 // Parameters
 4 int nx = 32;
 5 real ny = 16;
 6 real N = nx*ny;
 7 func 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
11 mesh 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
16 fespace Vh(Th,P1);
17 Vh<complex> u = f1, v;
18 Vh w = f1;
19 Vh 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)
24 v[] = dfft(u[], ny, -1);
25 u[] = dfft(v[], ny, +1);
26 cout << "||u||_\infty " << u[].linfty << endl;
27 
28 u[] *= 1./N;
29 cout << "||u||_\infty " << u[].linfty << endl;
30 
31 ur = real(u);
32 
33 // Plot
34 plot(w, wait=1, value=1, cmm="w");
35 plot(ur, wait=1, value=1, cmm="u");
36 v = w - u;
37 cout << "diff = " << v[].max << " " << v[].min << endl;
38 assert( norm(v[].max) < 1e-10 && norm(v[].min) < 1e-10);
39 
40 // Other example
41 //FFT Lapacian
42 //-\Delta u = f with biperiodic condition
43 func f = cos(3*2*pi*x)*cos(2*2*pi*y);
44 func ue = (1./(square(2*pi)*13.))*cos(3*2*pi*x)*cos(2*2*pi*y); //the exact solution
45 Vh<complex> ff = f;
46 Vh<complex> fhat;
47 Vh<complex> wij;
48 
49 // FFT
50 fhat[] = 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 ?:
54 wij = square(2.*pi)*(square(( x<0.5?x*nx:(x-1)*nx)) + square((y<0.5?y*ny:(y-1)*ny)));
55 wij[][0] = 1e-5; //to remove div / 0
56 fhat[] = fhat[] ./ wij[];
57 u[] = dfft(fhat[], ny, 1);
58 u[] /= complex(N);
59 ur = real(u); //the solution
60 w = real(ue); //the exact solution
61 
62 // Plot
63 plot(w, ur, value=1, cmm="ue", wait=1);
64 
65 // Error
66 w[] -= ur[];
67 real err = abs(w[].max) + abs(w[].min);
68 cout << "err = " << err << endl;
69 assert(err < 1e-6);
70 
71 fftwplan p1 = plandfft(u[], v[], ny, -1);
72 fftwplan p2 = plandfft(u[], v[], ny, 1);
73 real ccc = square(2.*pi);
74 cout << "ny = " << ny << endl;
75 map(wij[], ny, ccc*(x*x+y*y));
76 wij[][0] = 1e-5;
77 plot(wij, cmm="wij");

Complex

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

Output of this script is:

 1 Standard 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)
 6 Real part of (2.45+5.33i) = 2.45
 7 Imaginary part of (2.45+5.33i) = 5.33
 8 abs((2.45+5.33i)) = 5.86612
 9 Polar coordinates of (2.45+1.41421i) =
10 de Moivre formula: (2.45+1.41421i)^3 =
11  and polar(2.82887, 0.523509) = (2.45+1.41421i)
12 Conjugate 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
 2 string tt = "toto1" + 1 + " -- 77";
 3 
 4 // Append
 5 string t1 = "0123456789";
 6 t1(4:3) = "abcdefghijk-";
 7 
 8 // Sub string
 9 string t55 = t1(4:14);
10 
11 cout << "tt = " << tt << endl;
12 
13 cout << "t1 = " << t1 << endl;
14 cout << "t1.find(abc) = " << t1.find("abc") << endl;
15 cout << "t1.rfind(abc) = " << t1.rfind("abc") << endl;
16 cout << "t1.find(abc, 10) = " << t1.find("abc",10) << endl;
17 cout << "t1.ffind(abc, 10) = " << t1.rfind("abc",10) << endl;
18 cout << "t1.length = " << t1.length << endl;
19 
20 cout << "t55 = " << t55 << endl;

The output of this script is:

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

Elementary function

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

Fig. 221 Mesh

Array

  1 real[int] tab(10), tab1(10); //2 array of 10 real
  2 //real[int] tab2; //bug: array with no size
  3 
  4 tab = 1.03; //set all the array to 1.03
  5 tab[1] = 2.15;
  6 
  7 cout << "tab: " << tab << endl;
  8 cout << "min: " << tab.min << endl;
  9 cout << "max: " << tab.max << endl;
 10 cout << "sum: " << tab.sum << endl;
 11 
 12 tab.resize(12); //change the size of array tab to 12 with preserving first value
 13 tab(10:11) = 3.14; //set values 10 & 11
 14 cout << "resized tab: " << tab << endl;
 15 
 16 tab.sort ; //sort the array tab
 17 cout << "sorted tab:" << tab << endl;
 18 
 19 real[string] tt; //array with string index
 20 tt["+"] = 1.5;
 21 cout << "tt[\"a\"] = " << tt["a"] << endl;
 22 cout << "tt[\"+\"] = " << tt["+"] << endl;
 23 
 24 real[int] a(5), b(5), c(5), d(5);
 25 a = 1;
 26 b = 2;
 27 c = 3;
 28 a[2] = 0;
 29 d = ( a ? b : c ); //for i = 0, n-1 : d[i] = a[i] ? b[i] : c[i]
 30 cout << " d = ( a ? b : c ) is " << d << endl;
 31 d = ( a ? 1 : c ); //for i = 0, n-1: d[i] = a[i] ? 1 : c[i]
 32 d = ( a ? b : 0 ); //for i = 0, n-1: d[i] = a[i] ? b[i] : 0
 33 d = ( a ? 1 : 0 ); //for i = 0, n-1: d[i] = a[i] ? 0 : 1
 34 
 35 int[int] ii(0:d.n-1); //set array ii to 0, 1, ..., d.n-1
 36 d = -1:-5; //set d to -1, -2, ..., -5
 37 
 38 sort(d, ii); //sort array d and ii in parallel
 39 cout << "d: " << d << endl;
 40 cout << "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:

  1 tab: 10
  2     1.03    2.15    1.03    1.03    1.03
  3     1.03    1.03    1.03    1.03    1.03
  4 
  5 min: 1.03
  6 max: 2.15
  7 sum: 11.42
  8 resized 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
 12 sorted 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
 16 tt["a"] = 0
 17 tt["+"] = 1.5
 18  d = ( a ? b : c ) is 5
 19       2   2   3   2   2
 20 
 21 d: 5
 22      -5  -4  -3  -2  -1
 23 
 24 ii: 5
 25       4   3   2   1   0
 26 
 27 A1(2:10): 9
 28       2   3   4   5   6
 29       7   8   9  10
 30 A2(2:3:10): 9
 31       2   3   4   5   6
 32       7   8   9  10
 33 1:2:5 => 3
 34       1   3   5
 35 A1(2:10): 9
 36       2   3   4   5   6
 37       7   8   9  10
 38 A2(2:3:10): 9
 39       2   3   4   5   6
 40       7   8   9  10
 41 1.: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
 55 a: 5
 56       2   1   2   4   4
 57 
 58 b = a + a: 5
 59       4   2   4   8   8
 60 
 61 b += a: 5
 62       6   3   6  12  12
 63 
 64 b += 2*a: 5
 65      10   5  10  20  20
 66 
 67  b /= 2: 5
 68       5 2.5   5  10  10
 69 
 70 b .*= a: 5
 71      10 2.5  10  40  40
 72 
 73 b ./= a: 5
 74       5 2.5   5  10  10
 75 
 76 c = a + b: 5
 77       7 3.5   7  14  14
 78 
 79 c = 2*a + 4b: 5
 80      24  12  24  48  48
 81 
 82 c = a + 4b: 5
 83      22  11  22  44  44
 84 
 85 c = -a + 4b: 5
 86      18   9  18  36  36
 87 
 88 c = -a - 4b: 5
 89     -22 -11 -22 -44 -44
 90 
 91 c = -a -b: 5
 92      -7 -3.5     -7 -14 -14
 93 
 94 c = a .* b: 5
 95      10 2.5  10  40  40
 96 
 97 c = a ./ b: 5
 98     0.4 0.4 0.4 0.4 0.4
 99 
100 c = 2 * b: 5
101      10   5  10  20  20
102 
103 c = b * 2: 5
104      10   5  10  20  20
105 
106 ||a||_1 = 13
107 ||a||_2 = 6.40312
108 ||a||_infty = 4
109 sum a_i = 13
110 max a_i = 4 a[ 3 ] = 4
111 min a_i = 1 a[ 1 ] = 1
112 a' * a = 41
113 a quantile 0.2 = 2
114 b = a(I) : 5
115       2   4   4  -3   4
116 
117 c(I) = a 5
118      -3  -3   2   4   2
119 
120 b = a(I) : 5
121       2   4   4  -3   4
122 
123 c(I) = a 5
124      -3  -3   4   9   4
125 
126 A = 3 4
127        1   5   2   1
128        3   3   3   1
129        4   5   2   4
130 
131 C = 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 
136 A = 3 4
137        8  10  12  14
138       16  20  24  28
139       24  30  36  42
140 
141 B = (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}
144 3 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
 2 real f1 = 1.;
 3 real f2 = 1.5;
 4 
 5 // Mesh
 6 mesh Th1 = square(10, 10);
 7 mesh Th2 = square(10, 10, [1+x, -1+y]);
 8 plot(Th1, Th2);
 9 
10 // Fespace
11 fespace Uh1(Th1, P1);
12 Uh1 u1;
13 
14 fespace Uh2(Th2, P2);
15 Uh2 u2;
16 
17 // Macro
18 macro grad(u) [dx(u), dy(u)] //
19 
20 // Problem
21 varf 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 
31 varf 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     ;
40 matrix<real> Poisson1 = vPoisson1(Uh1, Uh1);
41 real[int] Poisson1b = vPoisson1(0, Uh1);
42 
43 matrix<real> Poisson2 = vPoisson2(Uh2, Uh2);
44 real[int] Poisson2b = vPoisson2(0, Uh2);
45 
46 //block matrix
47 matrix<real> G = [[Poisson1, 0], [0, Poisson2]];
48 set(G, solver=sparsesolver);
49 
50 //block right hand side
51 real[int] Gb = [Poisson1b, Poisson2b];
52 
53 // Solve
54 real[int] sol = G^-1 * Gb;
55 
56 // Dispatch
57 [u1[], u2[]] = sol;
58 
59 // Plot
60 plot(u1, u2);
../_images/BlockMatrix.png

Fig. 222 Result

Matrix operations

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

The output of this script is:

  1 f = 6
  2       0   0   0   0 0.5
  3       1
  4 g = 6
  5       0   1 1.224646799e-16   0   1
  6     1.224646799e-16
  7 A = # 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}
 10 6 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 
 36 m0 = 6
 37     -1.25   -2.25   0.5   0 5e+29
 38     1e+30
 39 m01 = 6
 40     -1.25   -2.25     0 0.25    5e+29
 41     1e+30
 42 m1 = 6
 43       0   0   0   0 0.5
 44     1.224646799e-16
 45 m2 = 6
 46     -nan      0   0 -nan    0.5
 47     8.165619677e+15
 48 dot Product = 0.5
 49 hermitien Product = (1.11022e-16,2.5)
 50 outer 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}
 53 6 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 
 63 hermitien 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}
 66 6 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 
 92 I = 8
 93       4   4   4   4   5
 94       5   5   5
 95 J = 8
 96       1   2   4   5   1
 97       2   4   5
 98 C = 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
103 D = # 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}
106 6 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

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

The output of this script is:

 1 5 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 
 8 5 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 
15 5 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 
22 5 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:

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

FE array

 1 // Mesh
 2 mesh Th = square(20, 20, [2*x, 2*y]);
 3 
 4 // Fespace
 5 fespace Vh(Th, P1);
 6 Vh u, v, f;
 7 
 8 // Problem
 9 problem 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 
20 Vh[int] uu(3); //an array of FE function
21 // Solve problem 1
22 f = 1;
23 Poisson;
24 uu[0] = u;
25 // Solve problem 2
26 f = sin(pi*x)*cos(pi*y);
27 Poisson;
28 uu[1] = u;
29 // Solve problem 3
30 f = abs(x-1)*abs(y-1);
31 Poisson;
32 uu[2] = u;
33 
34 // Plot
35 for (int i = 0; i < 3; i++)
36     plot(uu[i], wait=true);
FEArray1

Fig. 223 First result

FEArray2

Fig. 224 Second result

FEArray3

Fig. 225 Third result

Finite element array

Loop

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

Implicit loop

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

The output of this script is:

 1 0 1
 2 1 2
 3 2 3
 4 3 4
 5 4 5
 6 5 6
 7 6 7
 8 7 8
 9 8 9
10 9 10
11 b = 10
12       1   2   3   4   5
13       6   7   8   9  10
14 
15 a = 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 
27 ss = 1 1
28 2 2
29 3 5
30 
31  i   1          2
32  i  50          1
33 si = 1 2
34 50 100
35 
36 0 0 0.5
37 0 1 0.333333
38 0 2 0.25
39 0 3 0.2
40 1 0 0.333333
41 1 1 0.25
42 1 2 0.2
43 2 0 0.25
44 2 1 0.2
45 3 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}
49 10 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

 1 int i;
 2 cout << "std-out" << endl;
 3 cout << " enter i = ?";
 4 cin >> 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 
22 cout << i << endl;

File stream

 1 int where;
 2 real[int] f = [0, 1, 2, 3, 4, 5];
 3 real[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
19 func 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:

1 FreeFem++ script.edp arg1 arg2

The arguments can be used in the script with:

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

When using the command:

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

The arguments can be used in the script with:

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

Macro

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

The output script generated with macros is:

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

The output os this script is:

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

Basic error handling

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

The output of this script is:

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

Error handling

 1 // Parameters
 2 int nn = 5;
 3 func f = 1; //right hand side function
 4 func g = 0; //boundary condition function
 5 
 6 // Mesh
 7 mesh Th = square(nn, nn);
 8 
 9 // Fespace
10 fespace Vh(Th, P1);
11 Vh uh, vh;
12 
13 // Problem
14 real cpu = clock();
15 problem 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 
25 try{
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 }
37 catch(...) { //catch all error
38     cout << " Catch cholesky PB " << endl;
39 }

The output of this script is:

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