FreeFem++ note

Posted on 17/09/2018, in Infomation Technology, PhD.

I use this note for the FreeFem++ programming language. I have used it since the first year of my PhD thesis and I droped it 2 year laters because of some unexpected troubles. I then decide to use matlab instead. However, I keep using FreeFem++ in some simple cases just for comparing the computation with the codes done with matlab.

Work with PhD

• This is very interesting: numerical modeling transport problem with freefem - Florian De Vuyst.pdf

Miscellaneous

• verbosity = 0; at the beginning of file: only display the most useful
• Name of a file .edp must not contain a white space.
• Line break: just type enter as usual.
• Normal derivative

macro dn(u) (N.x*dx(u)+N.y*dy(u) ) //  def the normal derivative

• Find the area of each triangle K

mesh Th=square(10,10);
fespace Wh(Th,P0);
Wh etaK;
varf varea(unused,chiK) = int2d(Th,levelset=phi)(chiK);
etaK[] = varea(0,Wh);
cout << "etaK[]= " << etaK[] << endl;
//etaK[i] = area the ith triangle

real b(Vh.ndof);
b = a(0,Vh);

• Include cpp file inside a edp

• Normal vector’s direction

int1d(Th, levelset=phi)(dn(u)) // n=(1,0); //for Omega1 (on the left)
int1d(Th, levelset=-phi)(dn(u)) // n=(-1,0); //for Omega2 (on the right)

• Include another edp file into a edp file

include "<name>.edp"; // location is not important

• ff++ never uses label on the vertices

FreeFem++ with Visual Studio Code

• Install FreeFem++ as usual
• Terminal > Configure Default Build Task (choose the folder stick to this task)
• A task.json file opened and fill it with

{
// for the documentation about the tasks.json format
"version": "2.0.0",
{
"label": "Run FreeFem++",
"type": "shell",
"command": "FreeFem++ \${file}",
"group": {
"kind": "build",
"isDefault": true
},
}
]
}

• Press Ctrl + Shift + B.
• Assign keyboard shortcut to this task: Files > Preferences > Keyboard Shortcuts > Click on link “keybindings.json” and add following codes (change “ctrl+r” with the ones you want)

{
"key": "ctrl+r",
"args": "Run FreeFem++"
}

The order working with FreeFem++

1. Define a boundary

mesh Th= square(); // [0,1] x [0,1]

or an arbitary square [x0, x1] x [y0, y1]

real x0=0, x1=2;
real y0=0, y1=1;
int n=40, m=20; // number of fragments
mesh Th=square( n, m, [x0 + (x1-x0)*x, y0 + (y1-y0)*y] );

or

border <name>(t=a,b){x=f(t); y=g(t);};
mesh Th = buildmesh( <name>(nf) + <name>(nf) );

2. FE space

fespace Vh(Th,P1); // P0
Vh u = h(x,y);

3. Starting to work!

Plot

plot(Th, wait=1,fill=1,ps="SmileyFace.eps",bw=1,value=1, cmm="");
• wait: wait for click
• fill: fill figure
• ps: save plot to ps file
• bw: black and white
• value: display isoline
• cmm: insert comment for figure

Modify iso color range

real[int] viso(31);
for (int i=0; i<viso.n; i++){
viso[i] = i*0.3;
}

Mesh

Mesh connectivity

• index starts from 0.
• Vh.nt : number of elements/triangles
• Vh.ndof : number of dof (degree of freedom) = dim of space Vh = u.n
• Vh.ndofK : number of dof on 1 element.
• u.n : number of element of u
• Wh(k,i) : gives the number of ith degrees of freedom of element k.
• Th.nt : number of triangles
• Th.nv : number of vertices
• Th(i) : return the vextex i of Th
• Th(i).x: vertex i’s x-coor
• Th(i).y: vertex i’s y-coor
• Th[k] : return the triangle k of Th
• Th[i][j] : ith triangle’s jth vertex
• Th[i][j].x :
• Th[i][j].y :
• Th[i][j].label : node’s label
• Write

for (int i=0;i<nbtriangles;i++){ //all triangles
for (int j=0; j <3; j++){ //all vertices of the triangle i
cout <<i<<" "<<j<<" Th[i][j] = "<<Th[i][j]<<
" x = "<<Th[i][j].x<< " , y= "<<Th[i][j].y <<
", label=" << Th[i][j].label << endl;
}
}

• Th(0.55,0.6).nuTriangle : triangle’s number of some triangle.
• Th[it00].area : triangle’s area
• Th[it00].region: triangle’s region
• Th[it00].label: triangle’s label
• Diameter of a triangle

Ph h = hTriangle
h[].max

Structure of .msh file

• Save mesh
savemesh(Th, '<name>.msh');

• See the description of the structure of .msh at section 5.1.4 Data Structures and Read/Write Statements for a Mesh.
• They also used (like in matlab): the number of nodes, triangles and edges on boundary .
• The vertices are numbered in a counterclockwise orientation like in malab.
• Using getMeshFromFF.m to read the file .msh from ff++ (downloaded from here).

Region

• Insert into mesh

mesh Th2 = square(nn,nn,[2*x,3*y],region=2);        // RHS
mesh Th1 = square(nn,nn,[-x,3*y],region=1); // LHS
mesh Th = Th1 + Th2;
plot(Th,cmm="Th",wait=1);

• How to use? Th(0,0).region

Labels

for (int i=0;i<5;++i)
{
int[int] labs=[11,12,13,14];
mesh Th=square(3,3,flags=i,label=labs,region=10);
plot(Th,wait=1,cmm=" square flags = "+i );
}

P0, P1

• P0: Vh.nt = u.n = Vh.ndof
• P1: u.n = Vh.dof = dim(V), Vh.ndofK = 3

Connect 2 meshes

mesh Th2 = square(nn,nn,[2*x,3*y],region=2);        // RHS
mesh Th1 = square(nn,nn,[-x,3*y],region=1); // LHS
mesh Th = Th1 + Th2;
plot(Th,cmm="Th",wait=1);

intalledges, jump, mean

varf Ans(u,v)=
int2d(Th)(dx(u)*dx(v)+dy(u)*dy(v)  )
+ intalledges(Th)(//  loop on all  edge of all triangle
// the edge are see nTonEdge times so we / nTonEdge
// remark: nTonEdge =1 on border edge and =2 on internal
// we are in a triange th normal is the exterior normal
// def: jump = external - internal value; on border exter value =0
//      mean = (external + internal value)/2, on border just internal value
( jump(v)*mean(dn(u)) -  jump(u)*mean(dn(v))
+ pena*jump(u)*jump(v) ) / nTonEdge
)

Matrix

• In ff++ code, matrix starts from A(0,0) but in the input, it says (1,1)
• Matrix is sparse in ff++
• A.n dimension of A
• Dimension of matrix obtained from int2d(Th, levelset=phi) has the same dimension without levelset=phi (normal matrix obtained from int2d)
• cout << A gives

4 4 0 5 // row colum is-symmetric value
1 2 1 //value at row 1, column 2 is 1

• A matrix created by A = [[ 0, 1, 0, 10], [0, 0, 2,0], [0,0,0, 3], [ 4,0 , 0, 0]]; is also a sparse matrix whose A(0,0) doesn’t exist.
• [I,J,C] = A collects I indexes of rows, j indexes of columns, C values

int[int] I(1),J(1);
real[int] C(1);
[I,J,C] = A;
B = A(I,J); // B(i,j)=A(I(i),J(j)) //dim(B)=dim(I)*dim(J)
C = A(I^-1,J^-1); //B(I(i),J(j))=A(i,j)

• A.resize(10,20): resize matrix’s dimension
• Matrix’s diagonal

diagofA = A.diag; // get the diagonal of the matrix
A.diag = diagofA ; //set the diagonal of the matrix

Stiffness matrix

varf a(u,v) = int2d(Th)(dx(u)*dx(v) + dy(u)*dy(v)) + on(C,u=0);
varf b([v],[f]) = int2d(Th)(v*f);
matrix A = a(Vh,Vh); //stiffness matrix
matrix B = b(Vh,Vh);
F[] = B*f[]; //load vector
u[] = A^-1*F[];

Loops & condition

for(int i=1;i<n;i++){
commands;
}

while (int i<=n){
commands;
i=i+1;
}

if(condition){
statements;
}

Define a function

Function of a number

func real gamma(real x, real y){
if(x<x1/4 || x>3*x1/4 || y>y1/2)
return 1;
else
if((x>x1/4) && (x<3*x1/4) && (y<y1/2))
return -1;
else
return 0;
}

Function of an array

func real[int] name(real[int] &x){ //cai nay dung con tro thi phai
real[int] u;
....
return u; //hoac cung co the return x va khong can dinh nghia truoc u
}

Function in FE

func f = x^2 + y^2;
Vh fh = f; //fh is projection of f to Vh
plot(f); //error <-- can't plot
plot(fh); //plot is olny used under read/complex FE-function.

If a function returns real but takes values from x,y then it returns in Vh with diff values w.r.t values of x and y (cf. /nxfem/test/levelset/test_vortex_ffpp.edp)

func real velo(real x, real y, int t){
real val;
// lines of codes
return val;
}
Vh u;
real t = 1;
u = velo(x, y, t);

Array in FreeFem++

• index starts from 0
• tab.min, tab.max
• Change size: tab.resize(12)
• tab.sum
• tab.sort
• Others

real[string] text; //use string as an index
text["+"] = 1.5;

int[int] tt(2:10); // 2 to 10
int[int] t1(2:3:10); // 2 5 8

complex[int] tt(2.+0i:10.+0i); // 2,3,4,5,6,7,8,9,10
complex[int] t1(2.:3.:10.); // 2,5,8,
tt.re << endl ; // the real part array
tt.im << endl ; // the imag part array

Operators

• Logic comparison (used inside if)

&& = and
|| =  or
<,>,<=,>=,!=,==

• Logic inline operators

u = 5-(y>=0.3)*10; //if y>=3 is true then (y>=3)=1 else =0
u = (x>5)?1:2; //if x>5 is true, then u=1 else u=2

y=2+(x=5); // x=5 then y = 2+x = 7

x = 5+(3==2); // x = 5+0=5

Errors

• Comparison

real y = 1;
real y0 = 1/2;
real y1 = 1./2;
if (y0 == y/2) cout << "1"; else cout << "0"; //result 0
if (y0 == 1/2) cout << "1"; else cout << "0"; //result 1
if (y0 == 0.5) cout << "1"; else cout << "0"; //result 0
if (y1 == 0.5) cout << "1"; else cout << "0"; //result 1

• It’s strange

real dt = 0.1;
real t = 0;
real Tmax = 1;
for (int i=1; i< 11; i++){
t += dt;
}
bool ss;
ss = (t == Tmax);
cout << "t= " << t << endl; // 1
cout << "Tmax= " << Tmax << endl; // 1
cout << "ss= " << ss << endl; // False

Level set

int2d(Th,1,levelset=phi)()

The number 1 means region = 1.

mesh Th=square(n,m,[x0+(x1-x0)*x,y0+(y1-y0)*y],region=1);

The order

Vh phi=y-0.3;
int2(Th,levelset=-phi) --> upper (Omega1)
phi(x)>0 --> x in Omega1 (upper);
Top