FreeFem finite element language

Firstly, this is the color scheme of FreeFem++ for Gedit (tested on Ubuntu 12.04 and 14.04) : ffpp.lang This file should be placed in /usr/share/gtksourceview-3.0/language-specs/ (sudo).

This language is based on C++. It adds many features like:

  • Finite element spaces definition.

  • Integrals problems definition.

  • Data plot (curves, mesh, vectors, …​).

  • And many more…​

A description of this language is available on the FreeFem web site.

The presentation of this language will be directed toward the definition of a Stokes or Navier-Stokes problem.

Data types

  • int : integers

  • real : floats

  • complex : complex numbers

  • bool : booleans

  • string : characters string

  • matrix : matrix

  • border : mesh borders

  • mesh et mesh3 : 2D and 3D meshes

  • fespace : finite element spaces

  • func : functions

  • varf : variational forms

  • problem et solve : integrals problems

Vectors are equally defined with int[int], real[int], int[string] or real[string].

Meshes

The first step of the finite element method consists to create a consistent mesh in 2D or 3D. We will also start with this step. We will create a mesh of a pipe in 2D, then in 3D.

We start with the 2D mesh. The design principle is to define borders analitycally, to guide them and to create a mesh inside.

border gamma1(l=0, L){x=l; y=0.0; label = walls;};
border gamma2(l=0, D){x=L; y=l; label = out;};
border gamma3(l=0, L){x=L-l; y=D; label = walls;};
border gamma4(l=0, D){x=0; y=D-l; label = in;};
mesh Th=buildmesh(gamma1(L*nn) + gamma2(D*nn) + gamma3(L*nn) + gamma4(D*nn));

We define four border with the keyword border. They must be no hole at the edges. We declare our mesh with the keyword mesh. We use the buildmesh function to create the mesh inside borders, we define the number of nodes of each border (nn). Borders must be oriented in the trigonometric way, else the mesh is constructed outside the borders! The result is shown here. Lets move on to the 3D mesh. We will create here a base circle that we will extrude to create a cylinder.

border C(t=0,2*pi) { x = (D/2.0)*cos(t); y=(D/2.0)*sin(t); label=1;}
mesh Baseh = buildmesh(C(pi*D*nn));
int[int] rup=[0,out],  rdown=[0,in], rmid=[1,walls];
func zmin= 0;
func zmax= L;
mesh3 Th=buildlayers(
    Baseh,
    L*nn,
    zbound=[zmin,zmax],
    labelmid=rmid,
    reffaceup = rup,
    reffacelow = rdown
    );

We define the 2D mesh of a circle. After, we use the buildlayers function to create a 3D mesh, declared with the keyword mesh3. This function takes parameters:

  • Base mesh

  • Number of node on the extruded part

  • Bounds of extrusion

  • The reference of the extruded part(old reference, new reference)

  • The reference of the upper part

  • The reference of the lower part

Result is visible here

Problem definition

We will define the stationary non-compressible Stokes problem:

\[\begin{eqnarray} -\mu\Delta\mathbf{u} + \nabla p &= \textbf{f} &\mbox{ in }\Omega\label{eqn:stokes1}\\ \mathrm{div}(\mathbf{u}) &= 0 &\mbox{ in }\Omega\label{eqn:stokes2}\\ \mathbf{u} &= \textbf{g}_1 &\mbox{ on }\Gamma_D\label{eqn:stokes3}\\ \mu\frac{\partial\mathbf{u}}{\partial\mathbf{n}} &= \textbf{g}_2 &\mbox{ on }\Gamma_N\label{eqn:stokes4} \end{eqnarray}\]

All parameters and variables are defined in Stokes steady case section.

The (mixed) variational form of this problem is done by:

\[\begin{eqnarray} \mu\int_{\Omega}{\nabla\mathbf{u}\cdot\nabla\mathbf{v}} - \int_{\Omega}{p\mathrm{div}(\mathbf{v})} - \int_{\Gamma_N}{\textbf{g}_2\cdot\mathbf{v}} &=& \int_{\Omega}{\textbf{f}\cdot\mathbf{v}}\\ \int_{\Omega}{\mathrm{div}(\mathbf{u})q} &=& 0 \end{eqnarray}\]

That gives in FreeFem++:

mesh Th = readmesh(...);
fespace Xh(Th, P2);
fespace Vh(Th, P1);
Xh u1 = 0, v1;
Xh u2 = 0, v2;
Xh up1, up2;
Vh p = 0, q;
problem S ([u1, u2, p],[v1, v2, q])
    = int2d(Th)(
          mu*(dx(u1)*dx(v1) + dy(u1)*dy(v1) + dx(u2)*dx(v2) + dy(u2)*dy(v2))
        - p*q*1.e-9
        - p*dx(v1) - p*dy(v2)
        - dx(u1)*q - dy(u2)*q
    )
    - int2d(Th)(f1*v1 + f2*v2)
    + on (3, u1=0, u2=0)
    + on (1, u1=u1E, u2=u2E)
    ;

We start by import the mesh with the keyword readmesh. We next define the element finite spaces \(\mathbb{P}^1\) and \(\mathbb{P}^2\) (inf-sup condition).

We next define the variational form by using int2d function to define integrals and on to impose Dirichlet condition. To solve the problem, we write:

S;

Simply! And the Stokes problem is solved!

Plot and Save results

If the ratings of the previous algorithm is preserved, we display the solution with the command line:

plot([u1,u2]);
plot(p);

Plot function parameters :

  • value=bool (true ou false), to display or not the scale

  • ps=string, to save the picture in .ps format

  • coef=real, to define the length of vectors

  • nbiso=int, to define the number of iso-lines

  • And others present in the documentation…​

To save results, there are many solutions. We can save them in a file using a C++ method (write in a file). Or it exists functions already coded by the FreeFem contributors like iovtk that save solution in VTk format (viewable with Paraview). To use this format, you must first load the module in the FreeFem++ script with the command load "iovtk", then you can use:

int[int] fforder = [1, 1]
savevtk("NavierStokes.vtk",Th,[u1,u2,0],p,dataname="Velocity Pressure", order=fforder);

Parameters in this example are:

  • The name of the file to save

  • The mesh

  • The solution (in 3D so for an example in 2D the last component must be 0).

  • The data name (optional)

  • The order of the solution (optional)

Others functionality are explained in the Stokes steady case part.