See

discussions, stats, and author profiles for this publication at: https://www.researchgate.net/publication/260002309

The Finite Element Method with FreeFem++ for

beginners

Article in The Electronic Journal of Mathematics and Technology · June 2013

CITATION READS

1 4,008

2 authors:

Roberto Font Francisco Periago

Biometric Vox Universidad Politécnica de Cartagena
9 PUBLICATIONS 8 CITATIONS 32 PUBLICATIONS 237 CITATIONS

SEE PROFILE SEE PROFILE

Some of the authors of this publication are also working on these related projects:

Control of Quantum Wave Equations View project

All content following this page was uploaded by Roberto Font on 03 February 2014.

The user has requested enhancement of the downloaded file. All in-text references underlined in blue are added to the original document
and are linked to publications on ResearchGate, letting you access and read them immediately.
 The Electronic Journal of Mathematics and Technology, Volume 7, Number 4, ISSN 1933-2823

The Finite Element Method with FreeFem++ for

beginners

Roberto Font and Francisco Periago

e-mail: [email protected] - [email protected]
Department of Applied Mathematics and Statistics
Universidad Politécnica de Cartagena
30202 Cartagena
Spain

Abstract
These notes are concerned with the numerical resolution of Partial Differential Equations
(PDE) by the Finite Element Method (FEM). Emphasis is placed in the practical numerical res-
olution of this type of problems by using the free software FreeFem++. Theoretical background
is briefly reviewed (without entering in technical details) and a number of examples are given,
along with the code. The goal is to introduce in a simple way this difficult (but very interesting
in real-world applications) topic to undergraduate students who do not have a solid background
both in variational methods and in numerical analysis for PDE.

1 Introduction
Mathematical modeling is the art of representing a physical phenomenon in mathematical terms. Most
of these models (equations of linear elasticity, Navier-Stokes equations of fluid mechanics, Maxwell
equations of electromagnetism, etc...) are written as a partial differential equation (PDE) or a system
of PDE defined on a suitable domain of Rn , n = 1, 2, 3. Boundary (and possibly initial) conditions
complete the model. The interest in Industry of solving numerically these problems is nowadays
unquestionable. Therefore, these aspects of Mathematics should be considered as an essential part of
education in schools of mathematics, physics and engineering.
In these notes, we focus on the numerical resolution of PDE problems by the Finite Element
Method (FEM). Our goal is twofold: (i) first, we aim at introducing FEM to undergraduate students
who have received some very preliminary courses on linear algebra and calculus, and (ii) we shall
illustrate the power of this method by solving a number of specific examples with the free software
FreeFem++. Of course, there is a number of wonderful books (e.g., [6]) dealing with the numerical
resolution of PDE and there is also a tutorial for FreeFem++ [3]. These notes do not aim at replacing
or completing those references. We just will try to collect the main aspects of this topic in a simple
way in order to make it accessible to students without a solid background on the theoretical ingredients
underlying to FEM.
 The Electronic Journal of Mathematics and Technology, Volume 7, Number 4, ISSN 1933-2823

The rest of the paper is organized as follows. In Section 2 we briefly review on the FEM. Then, in
Section 3 we follow the Learning by Examples method to illustrate how to use FreeFem++ to solve a
number of interesting PDE problems by the FEM.

2 A review on the Finite Element Method

Next, we shall describe the FEM through a very simple one-dimensional model. Part of the material
of this section has been adapted from [4, 5].

2.1 A very simple one-dimensional model in linear elasticity

Consider a perfectly elastic and flexible string stretched along the segment [0, L] . Assume that on the
string acts a vertical force f = f (x) and that the string is fixed at both endpoints.
Hooke’s law on linear elasticity and the momentum conservation law lead to the equation − (κu0 )0 =
f in (0, L) , where κ = κ (x) ≥ κ0 > 0 depends on the physical properties of the string. Thus, the
mathematical model for this problem is

− (κu0 )0 = f

in (0, L)
(SP)
u (0) = u (L) = 0.

Assume now that the external force is located at x = L/2. It is clear that physically we would
observe something like in Figure 1.

x L/2
0 L

u(x)

Figure 1: String fixed at its ends carrying a load at x = L/2.

But a function as given in the preceding picture is not differentiable at x = L/2. Hence, such a
function cannot be a solution of a second-order differential equation − (κu0 )0 = f in (0, L) , at least
in the classical sense, that is, in the sense of being u twice differentiable at each point of the open
interval (0, L) and satisfying the differential equation point-wise in that interval. This simple example
shows that it is necessary to look for a new concept of a solution of a differential equation.

290
 The Electronic Journal of Mathematics and Technology, Volume 7, Number 4, ISSN 1933-2823

2.2 Variational formulation of a differential equation

Instead of looking at the conservation of momentum principle, we may start from the virtual work
principle. This way, if the force f produces a virtual displacement v (x) at a point x, then the work
produced by f (x) is f (x) v (x) so that the work produced by f along the whole string is given by
Z L
f (x) v (x) dx.
0

Proceeding in the same way in the left-hand side of the equation in problem (SP) and integrating by
parts, Z L Z L
0 0
− (κ (x) u (x)) v (x) dx = κ (x) u0 (x) v 0 (x) dx,
0 0

since v (0) = v (L) = 0. This process transforms the differential equation − (ku0 )0 = f into the
integral equation Z L Z L
0 0
κ (x) u (x) v (x) dx = f (x) v (x) dx, (1)
0 0
which as indicated above reflects virtual work principle in Physics. Hence, we may think in a solution
of problem (SP) as a function u : [0, L] → R, with u (0) = u (L) = 0, which satisfies (1) for all
possible displacements v such that v (0) = v (L) = 0. One of the main advantages of this approach
is that we have reduced the requirements that the function u ought to satisfy. Nevertheless, there are
still two main points which are not clear for the present moment:

(1) u should be at least once differentiable, but this is not the case for a charge located at a point.

(2) if the charge is located at x = L/2 (that is, f (x) = 0 for all 0 ≤ x ≤ L, x 6= L/2), then
RL
0
f (x) v (x) dx = 0 so that u = 0 satisfies (1), in contradiction with physical experience (see
Figure 1).

The above shows that it is necessary to understand better the way of dealing mathematically with
the concept of a charge located at a point. As a first approach, let us consider that the charge is
distributed in a small portion around x = L/2, that is,
L
− ε ≤ x ≤ L2 + ε

1/ (2ε) , 2
fε (x) =
0, otherwise

The work produced by this charge with a displacement v is now given by

L
Z L Z +ε
2 v (x)
fε (x) v (x) dx = dx = v (ξε ) ,
0 L
2
−ε 2ε

where L2 − ε ≤ ξε ≤ L2 + ε, and the equality being a consequence of the mean value theorem for
integrals. Now letting ε → 0 and assuming v is continuous it is concluded that, in this sense, the
work produced by a single charge located at L/2 and producing a displacement v = v (x) is v (L/2) .
This reasoning shows that we can deal mathematically with a single charge located at a point x0 as

291
 The Electronic Journal of Mathematics and Technology, Volume 7, Number 4, ISSN 1933-2823

a mapping, say δx0 , which acts on a certain class of functions and produce numbers by following the
rule
δx0 : v 7→< δx0 , v >=def v (x0 ) .
All of this makes a rigorous sense in the framework of the theory of distributions developed by the
French mathematician L. Schwartz in 1946, and this is basically the mathematical background we
need to formulate precisely our boundary value problem. Next we briefly review some of the main
aspects of this theory.
Given an open interval I ⊆ R, D (I) stands for the test functions space composed of all functions
v : I → R which are indefinitely differentiable and have compact support1 in I. Now given a sequence
{vn }∞
n=1 ⊆ D (I) it is said that this sequence converges to v ∈ D (I) if (i) there exists a compact set
dm
K ⊆ I such that supp vn ⊆ K for all n ∈ N, and (ii) for each m ∈ N, dxm vn converges uniformly to
dm
dxm
v.
A distribution u is defined as a linear and continuous (with respect to the above notion of conver-
gence) mapping u : D (I) → R. Let us now show the two typical examples of distributions.

• If f : I → R is locally integrable, then a distribution uf may be associated with f by means of

the mapping Z
def
uf : D (I) → R, v 7→< uf , v >= f (x) v (x) dx.
I
It is very common to denote by f the distribution uf . In what follows we use this notation.

• Given x0 ∈ I, it is not hard to show that the mapping

δx0 : D (I) → R, v 7→< δx0 , v >=def v (x0 )

is a distribution. This distribution is called the Dirac delta (or Dirac mass).

Next, we wish to do calculus with distributions. In particular, we wonder if it is possible to

define the derivative of a distribution. Assume first that f : I → R is a differentiableR function
and take v ∈ D (I) . From the distributions point of view, the derivative of f is given by f 0 v, but
integrating
R 0 Rby parts and taking into account that v vanishes at the extremes of I it is concluded that
0
f v = − f v , that is, in the language of distributions

< f 0 , v >= − < f, v 0 > .

In general, given a distribution u, its distributional derivative is defined as the distribution u0 given by

< u0 , v >= − < u, v 0 > for all v ∈ D (I) . (2)

Notice that, contrary to what happens with functions, a distribution may be derived infinitely times.
Once we have at one’s disposal distributions, let us go back to (SP) or in particular to the integral
equation (1). Let L2 (0, L) denote the Hilbert space of functions f : (0, L) → R of square integrable
R 1/2
L 2
endowed with the norm kf k2 = 0 f (x) dx . In order to the right-hand side of (1) to make

1
the support of a continuous function v is defined as supp v = {x ∈ I : v (x) 6= 0}.

292
 The Electronic Journal of Mathematics and Technology, Volume 7, Number 4, ISSN 1933-2823

sense, by the Cauchy-Schwarz inequality, all we need is to have f, v ∈ L2 . As for the left-hand side
of (1), if we assume κ to be bounded, then we just need u0 , v 0 ∈ L2 . As usual, we denote by

H01 (0, L) = u ∈ L2 (0, L) : u0 ∈ L2 (0, L) , u (0) = u (L) = 0 ,



1/2
which endowed with the norm kukH 1 = kuk22 + ku0 k22 is a Hilbert space. Of course, the deriva-
0
tives in this space are understood in the sense of distributions. It can be proved that the test functions
space D (0, L) is k·kH 1 −dense in H01 (0, L) . From this it follows that the Dirac delta can be extended
0
to a linear and continuous form δx0 : H01 (0, L) → R, acting as < δx0 , u >= u (x0 ) , with u ∈ H01 .
Functions f of L2 may be also considered as linear and continuous forms on H01 just by defining it as
Z
1 def
f : H0 → R, v 7→< f, v >= f v.

These new ideas lead to a new formulation of (SP). Precisely, given a linear and continuous form
f : H01 → R and the bilinear form
Z L
1 1
a : H0 × H0 → R, (u, v) 7→ a (u, v) = κ (x) u0 (x) v 0 (x) dx
0

a function u ∈ H01 is said to be a weak solution of (SP) if the identity

a (u, v) =< f, v >

holds for all v ∈ H01 . The existence and uniqueness of a weak solution is usually obtained through a
Lax-Milgram type theorem [2, p. 297].

2.3 The Finite Element Method

We are now in a position to describe the FEM. We already know that the solution of our problem
lives in a Hilbert space like H01 . The main difficulty is that this space is too big. Precisely, it is
infinite-dimensional. The basic idea of the FEM consists of approximating the big space H01 by some
appropriate finite-dimensional spaces Hh , where h is a positive parameter, satisfying the following
conditions:

(a) Hh ⊆ H01 ,
(b) in Hh we can easily solve the variational problem and hence obtain a solution uh , and
(c) when h & 0, Hh % H01 .

This procedure leads to a numerical approximation uh of the weak solution u of our original
problem. Let us show how this method works for the case of the (SP) problem.
Construction of Hh
Let n ∈ N and h = 1/ (n + 1) . The interval [0, L] can be decomposed as
n
[
[0, L] = [ci, ci+1 ] , ci = ih, 0 ≤ i ≤ n.
i=0

293
 The Electronic Journal of Mathematics and Technology, Volume 7, Number 4, ISSN 1933-2823

Now set
n o
Hh = v : [0, L] → R continuous, v (0) = v (L) = 0 and v|[ci ,ci+1 ] ∈ P1

where P1 stands for the space of polynomials of degree less than or equal to 1. It is not hard to show
that Hh ⊆ H01 . As for the dimension of Hh , the family of functions
1 − |x−c i|

, ci−1 ≤ x ≤ ci+1
φi (x) = h
0, otherwise
for 1 ≤ i ≤ n, is a basis of Hh . Therefore, dim (Hh ) = n.
The construction of the space Hh is not unique. Other options are possible as well. Anyway, we
are considering here the simplest case.
The variational problem in Hh  
Since Hh ⊆ H01 is finite-dimensional, Hh , k·kH 1 is a Hilbert space. Thus, the variational
0

problem (SP) in Hh is formulated as follows: find uh = ni=1 uih φi ∈ Hh such that the identity
P
Z L
a (uh , vh ) =def
κu0h vh0 dx = < f, vh > (3)
0
RL
holds for all vh ∈ Hh . Remember that if f ∈ L2 , then < f, vh > = 0 f vh ; and if f is a Dirac delta
at x0 , then < f, vh > = vh (x0 ) .
It is evident that if (3) holds for all vh ∈ Hh , then it is held for all φi , 1 ≤ i ≤ n; and conversely,
since Hh is a vector space, if (3) holds for all φi , 1 ≤ i ≤ n, then it is held for all vh ∈ Hh . Hence,
the integral equation (3) transforms into the linear system
n
X
uih a (φi , φj ) =< f, φj >, 1 ≤ j ≤ n,
i=1

and taking into account the properties of the functions φi , this system is written as


 a (φ1 , φ1 ) u1h + a (φ2 , φ1 ) u2h + 0 = < f, φ1 >
 ··············· ··· ···


0 + a (φi−1 , φi ) ui−1
h + a (φ i , φ i ) u i
h + a (φi+1 , φi ) ui+1
h + 0 = < f, φi >
··············· ··· ···




0 + a (φn−1 , φn ) un−1 + a (φ , φ ) u n
= < f, φn > .

h n n h

Note that since the matrix of this system is three-diagonal, a numerical method for linear systems such
as the Cholesky method is suitable to solve it. As for the computation of the integrals in a (φi , φj ) and
in < f, φi >, this may be done by the mid-point rule. In higher dimensions a number of well-known
methods may be used for the numerical computation of these integrals (see for instance [8]).

2.4 Extension to higher dimensions

In order to illustrate the FEM in higher dimensions we consider the following Poisson equation as the
test model. 
−∇ · (c∇u) + au = f in Ω
∂u (4)
c ∂n + qu = g on ∂Ω.

294
 The Electronic Journal of Mathematics and Technology, Volume 7, Number 4, ISSN 1933-2823

Here, Ω ⊂ RN is a bounded domain with smooth boundary ∂Ω  (N = 2, 3 in real-world

 applications),
∂u ∂u
x = (x1 , · · · , xN ) ∈ Ω, u = u (x) is the unknown , ∇u = ∂x1 , · · · , ∂xN is the gradient of u,
c = c (x), a = a (x) and f = f (x) are given functions defined in Ω, q = q (x) and g = g (x) are
∂u
given functions defined on ∂Ω, and finally ∂n = ∇u · n is the directional derivative of u with respect
to the outward unit normal vector n to ∂Ω. Also we recall that ∇· (also denoted by div) denotes the
divergence operator: given a vector field F = (F1 , · · · , FN ) ,
∂F1 ∂FN
div (F ) ≡ ∇ · F = + ··· + .
∂x1 ∂xN
Finally, we also recall the definition of the Laplacian operator: given an scalar function u = u (x) ,
the Laplacian of u, denoted by ∆u (also by ∇2 u) is given by
∂ 2u ∂ 2u
∇2 u ≡ ∆u = div(∇u) = + · · · + .
∂x21 ∂x2N
To obtain a variational formulation of (4) we proceed as in the one-dimensional case. First we
multiply the PDE by some function v and then integrate in Ω. We obtain
Z Z
[−∇ · (c∇u) + au] v dx = f v dx.
Ω Ω

Next, we apply Green theorem (integration by parts formula in higher dimensions) to get
Z Z Z
(c∇u · ∇v + auv) dx − (n · c∇u) v dΓ = f v dx,
Ω ∂Ω Ω

where the second integral in this expression is a line integral (N = 2) or a surface integral (N = 3).
Finally, we impose the boundary condition. Thus,
Z Z Z
(c∇u · ∇v + auv) dx − (−qu + g) v dΓ − f v dx = 0 ∀v. (5)
Ω ∂Ω Ω

As in the 1D case, both u and v belong to a suitable function space, but it is not necessary to enter in
technical details in this first contact with FEM.
Also, a discretization procedure as in the 1D case applies. The generation of a good mesh and the
definition of the corresponding finite element space may be, in practise, a very difficult task. For more
details on this passage see [3, 5, 6]. Finally, we notice that we do not consider in this introductory
paper the analysis of the order of convergence in the FEM. We refer the reader to [1], for instance.

3 Learning FreeFem++ by examples

FreeFem++ is a free software to solve PDE using the FEM in dimensions 2 and 3 (also 1D problems
may be solved from 2D). It has been (and it is being) developed by the Laboratoire Jacques-Louis
Lions (Université Pierre et Marie Curie, France). FreeFem++ runs on Windows, Linux and Mac
OS. The download page is http : //www.freefem.org. A detailed description of the possibilities of
FreeFem++ as well as a tutorial and a number of slides of courses may be found in this web site.

295
 The Electronic Journal of Mathematics and Technology, Volume 7, Number 4, ISSN 1933-2823

In order to make accessible this wonderful software to students (mainly of a master level, but
also to undergraduate students in Physics, Mathematics or Engineering) who have not had a previous
experience with the numerical resolution of PDE, next we show a number of illustrative examples.
Some of the examples that follow have been taken (with a few small changes) from FreeFem++
tutorial. They run under FreeFem++-cs 12.4 32 version.

3.1 Poisson’s equation on the unit disk

As a first example, let us consider the Poisson equation on the unit disk with zero Dirichlet boundary
conditions. In its classical form, the problem reads as

−∆u = f in Ω = {(x, y) ∈ R2 : x2 + y 2 < 1}
(6)
u=0 on ∂Ω,

where u = u (x, y) is the unknown and f is the source term.

This problem is an acceptable mathematical model for a number of physical problems. For in-
stance, (6) models the shape of a soap film which is glued to a ring in a plane. The shape is described
as the graph (x, y, u(x, y)) of the vertical displacement u(x, y) under a vertical pressure p in terms of
force per unit area and an initial tension µ in terms of force per unit length. This way, f = p/µ in
(6). The boundary condition u = 0 on ∂Ω in (6) indicates that the soap film is glued to the ring. An-
other possible application of (6) is in electrostatics where u is an electrostatic potential and f = ρ/,
with ρ the space charge density and  the dielectric coefficient. Similar applications appear in fluid
mechanics or stationary heat conduction.
The variational formulation of (6) reduces to find u ∈ H01 (Ω) such that
Z Z
∇u · ∇v dxdy = f v dxdy ∀v ∈ H01 (Ω) .
Ω Ω

Next we include the code in FreeFem++ to solve this problem and comment on it. For comparison
1−(x2 +y 2 )
reasons we take f = 1. Hence, (6) has an exact solution which is given by u (x, y) = 4
. We
2
may use this expression to analyze the evolution of the error (in the L −norm) with respect to the
type of finite elements and/or the number of nodes of the mesh. All scripts files in FreeFem++ must
be saved with extension .edp.
Figure 2 shows the picture of the solution obtained with the code below. As for the error, we
obtain ku − uh kL2 (Ω) = 5.19123 × 10−6 .

// Poisson’s equation on the 2D unit disk

border C(t=0,2*pi){x=cos(t); y=sin(t);} //boundary of the domain
mesh Th=buildmesh(C(1000)); // mesh with 100 points on the boundary
savemesh(Th,’’mesh_malla1.msh’’); // to save the mesh data
plot(Th,ps=’’malla.eps’’); // to plot and save the mesh
fespace Vh(Th,P1); // P1 Lagrange finite elements
Vh uh,vh; // uh,vh belong to Vh
func f=1;//source term
func u=0.25*(1-xˆ2-yˆ2); // exact solution

296
 The Electronic Journal of Mathematics and Technology, Volume 7, Number 4, ISSN 1933-2823

real L2error;
// solving the variational formulation of the problem
solve Poisson(uh,vh)= int2d(Th)(dx(uh)*dx(vh)+dy(uh)*dy(vh))
-int2d(Th)(f*vh)
+on(C,uh=0);
real [int] viso=[0.,0.05,0.1,0.15,0.2,0.25];
plot(uh,value=true,fill=true,ShowAxes=0,ColorScheme =1,
viso=viso,ps="fig_poisson2d.eps"); //to plot and save the solution
L2error=sqrt(int2d(Th)((u-uh)ˆ2)); // error in Lˆ2-norm
cout<<’’L2error = ’’<<L2error<<endl; // to show the error

Figure 2: Picture of the solution of problem (6) obtained with the above code.

3.2 Poisson’s equation on the 3D unit cube

In this example, we also consider the Poisson equation but in 3D and with a Robin-type boundary
condition. There are different ways of building a 3D mesh in FreeFem++. Here we consider the
function buildlayers in which a 3D domain is generated from a 2D one by making a sequence
of horizontal layers. In our case, we take the 3D unit cube [0, 1]3 which is generated from the 2D unit
cube [0, 1]2 . The problem is

in Ω3D = [0, 1]2 × [0, 1]


 −∆u = f
∂u
= α(ue − u) on face number 2 (7)
 ∂n∂u
∂n
=0 on the rest of the boundary,
where now u = u (x, y, z) and face number 2 is the face {x = 1}. In its weak or variational form, the
problem reads as
Z Z Z Z
∇u · ∇v dxdydx + αuv dΓ − αue v dΓ − f v dxdydz = 0 ∀v.
Ω face no 2 face no 2 Ω

297
 The Electronic Journal of Mathematics and Technology, Volume 7, Number 4, ISSN 1933-2823

Note that now v = v (x, y, z) does not vanish on the boundary.

// Poisson’s equation on the 3D unit cube

load ’’msh3’’ // we need msh3 in 3D
mesh Th2=square(10,10);//2D mesh
int ncapas=10; // number of horizontal layers on the z-axis
real zmin=0, zmax=1;
mesh3 Th=buildlayers(Th2,ncapas,zbound=[zmin,zmax]); //3D mesh
real alpha=0.25, ue=25;
func 1000.*(x>0.4)*(x<0.6)*(y>0.4)*(y<0.6)*(z>0.4)*(z<0.6);
fespace Vh(Th,P13d); //P1 Lagrange finite elements in 3D
Vh u,v; // u,v belong to Vh
// and now we use macros to simplify expressions
macro grad(u) [dx(u),dy(u),dz(u)] //EOM

solve lapace3d(u,v)=int3d(Th)(grad(u)’*grad(v))
+int2d(Th,2)(alpha*u*v)
-int3d(Th)(f*v)
-int2d(Th,2)(alpha*ue*v);
plot(u);

A possible physical interpretation of this problem is the following: A body that occupies the region
Ω has an internal heat source f . All the faces of the cube Ω are isolated except the face x = 1 where
∂u
there is a convective heat flux modelled by Newton’s law of cooling: heat transfer (κ ∂n ) takes place
at a rate proportional (−α) to the difference of temperature between the solid (u) and the surrounding
(ue ). The thermal conductivity of the material, κ, has been normalized to 1. u(x, y, z) represents the
temperature at the point (x, y, z) at the stationary state. Figure 3 displays a cut plane picture of the
solution. It is observed that temperature is higher on the internal region where heat is generated an
decreases as we approach the face x = 1 where there exists a heat transfer to the surrounding.

3.3 Transient heat equation

FreeFem++ has two building blocks: the elliptic PDE system solver problem solve, and the convec-
tion module convect. However, there is also the possibility of implementing our own algorithms
like the following example shows. Precisely, we shall solve the following evolution problem for the
heat equation:  ∂u
 ∂t − ∇ · (κ∇u) = 0 in Ω × (0, T )
u (x, y, 0) = u0 (x, y) in Ω (8)
 ∂u
κ ∂n + α (u − ue ) = 0 on ∂Ω × (0, T ) .
Here u = u (x, y, t) stands for the temperature at position (x, y) and time t, and u0 (x, y) is the initial
temperature. For simplicity, we take Ω = (0, 6) × (0, 1) . The thermal conductivity is given by

1 if y < 0.5
κ (x, y) =
0.2 else.

298
 The Electronic Journal of Mathematics and Technology, Volume 7, Number 4, ISSN 1933-2823

Figure 3: Picture of the solution of problem (7) obtained with the above code.

As in the preceding example, the boundary condition reflects Newton’s law of cooling. To obtain the
variational formulation of (8) we proceed as in the stationary case. To begin with, we multiply the
PDE by some test function v = v (x, y) . Then, we integrate over Ω and apply integration by parts
formula. We get Z Z Z
∂u ∂u
v dxdy − κ v dΓ + κ∇u · ∇v dxdy = 0.
Ω ∂t ∂Ω ∂n Ω
Finally, we impose the boundary condition to obtain
Z Z Z
∂u
v dxdy + α (u − ue ) v dΓ + κ∇u · ∇v dxdy = 0 ∀v. (9)
Ω ∂t ∂Ω Ω

Next, we discretize the time variable by using an implicit Euler scheme. For the convergence proper-
ties of this scheme we refer the reader to [6]. Let us denote by dt = T /N, with N ∈ N fixed , and by
un = u (tn ) , with tn = ndt, 0 ≤ n ≤ N. Then, we have
∂u un − un−1
≈ .
∂t dt
Substituting this expression into (9),
un − un−1
Z Z Z
v dxdy + α (un − ue ) v dΓ + κ∇un · ∇v dxdy = 0 ∀v,
Ω dt ∂Ω Ω

which is the discrete weak formulation that we implement in FreeFem++ next.

299
 The Electronic Journal of Mathematics and Technology, Volume 7, Number 4, ISSN 1933-2823

// transient heat equation

func u0=10+90*x/6; //initial temperature
func k=1.8*(y<0.5)+0.2; // thermal conductivity
real ue=25; //outside temperature
real alpha=0.25;
real T=0.7; //final time
real dt=0.1; //time step
mesh Th=square(30,5,[6*x,y]); // domain
fespace Vh(Th,P1); // Lagrange P1 finite elements
Vh u=u0,v,uold;
Vh[int] sol((T/dt)+1);
// variational formulation. Implicit Euler scheme
problem calor(u,v)=int2d(Th)(u*v/dt+k*(dx(u)*dx(v)+dy(u)*dy(v)))
+int1d(Th,1,2,3,4)(alpha*u*v)
-int1d(Th,1,2,3,4)(alpha*ue*v)
-int2d(Th)(uold*v/dt);
ofstream ff(’’calor.dat’’); // to generate a data file
for(int j=0;j<T/(dt);j=j+1)
{
uold=u;
calor;
sol[j+1]=u; // to store the solution at different times
ff<<u(3,0.5)<<endl;// calor.dat saves temperature at the point (3,0.5)
}
plot(sol[0],value=true,ShowAxes=0,ColorScheme = 1,ps="sol_t0.eps");
plot(sol[5],value=true,ShowAxes=0,ColorScheme = 1,ps="sol_t05.eps");
plot(sol[10],value=true,ShowAxes=0,ColorScheme = 1,ps="sol_t1.eps");
plot(sol[20],value=true,ShowAxes=0,ColorScheme = 1,ps="sol_t2.eps");

Figure 4 shows the iso-thermal lines of the solution of problem (8) at some discrete times. These
pictures simulate how air cools a rectangular plate which has an initial temperature u0 (x, y) and a
discontinuous thermal conductivity coefficient. Notice that at time t = 0, the temperature at the right-
hand edge is 100 and the temperature at left-hand edge is 10. Then, as time evolves there is a heat
flux from the heat region to the cold one.

3.4 Lamé system in 2D

In this example we move on a system of PDE, precisely, the Lamé system of linear elasticity in the
plane. Consider a plate as shown in Figure 5 (top). Assume that the left vertical side of the plate is
fixed and that there is no boundary force on the lower, upper, right sides and on the boundary of the
interior circle. The body force is given by f = (0, −1) . When the displacement u = (u1 , u2 ) is small

300
 The Electronic Journal of Mathematics and Technology, Volume 7, Number 4, ISSN 1933-2823

Figure 4: From left to right and from top to bottom, pictures of the FEM solution of problem (8) at
times t = 0 (top left), t = 0.5, t = 1 and t = 2 (bottom right).

and the material is (linear) elastic, in the stationary case the following conditions are satisfied

 −∇ · σ = f in Ω
u=0 on left side (10)
σ·n=0 on the rest of the boundary,


with
σij (u) = λδij trace (ε) + 2µεi,j , 1 ≤ i, j ≤ 2,
the stress tensor, where δij = 1 if i = j, 0 otherwise, λ, µ are the Lamé coefficients which describe
the mechanical properties of the material, and
 
1 ∂ui ∂uj
εij (u) = + , 1 ≤ i, j ≤ 2
2 ∂xj ∂xi

is the linearized strain tensor. As usual n is the outward unit normal vector. The variational form of
(10) is Z Z
σ (u) : ε (v) dxdy − f · v dxdy = 0 ∀v such that v = 0 on left side.
Ω Ω
P
The notation σ (u) : ε (v) = i,j σij εij stands for the contraction of σ and ε.

301
 The Electronic Journal of Mathematics and Technology, Volume 7, Number 4, ISSN 1933-2823

// Linear elasticity in 2D

//Domain and mesh

border a1(t=0,20){x=t;y=-1;label=1;}; //Lower side
border a2(t=-1, 1){x=20;y=t;label=2;}; //Right side
border a3(t=0,20){x=20-t;y=1;label=3;}; //Upper side
border a4(t=-1,1){x=0;y=-t;label=4;}; //Left side
border C(t=0, 2*pi){x=10+0.5*cos(t); y=0.5*sin(t);}; //Interior circle
//We define the meh from a1,a2,a3,a4 and C
mesh Th=buildmesh(a1(100)+a2(10)+a3(100)+a4(10)+C(-10));
plot(Th,wait=true); //Plots the mesh

fespace Vh(Th,P2); //Finite elements (P2) space definition

Vh u,v,uu,vv; //Variables difinition

//Some constants and physical parameters

real sqrt2=sqrt(2.);
real E = 21e5, nu = 0.28; //E: Young’s modulus, nu: Poisson’s ratio
real lambda = E*nu/((1+nu)*(1-2*nu)), mu= E/(2*(1+nu)); //Lame coefficients
real f = -1; //Body force

//We can define macros to make programming easier

macro epsilon(u,v) [dx(u),dy(v),(dy(u)+dx(v))/sqrt2] // EOM
macro div(u,v) (dx(u)+dy(v)) // EOM

//Solving the problem in variational form

solve lame([u,v],[uu,vv])= int2d(Th)(
lambda*div(u,v)*div(uu,vv)
+2.*mu*( epsilon(u,v)’*epsilon(uu,vv) ) )
- int2d(Th)(f*vv)
+ on(4,u=0,v=0);

//Ploting the result

real coef=100; //A coefficient of amplification
mesh Thd = movemesh(Th, [x+u*coef, y+v*coef]);
plot(Thd,wait=1); //Result as a deformed mesh
real dxmin = u[].min; //Minimum of displacement in x direction
real dymin = v[].min; //Minimum of displacement in y direction

Sometimes it is interesting to obtain magnitudes that are not the solution of our PDE but can be
evaluated from this solution. In the case of linear elasticity, it might be interesting to obtain, for
example, the value of the von Mises stress. The von Mises stress is frequently used in engineering as
a scalar representation of a 3–dimensional load state. This way, the yielding under any condition can
be predicted from the results of simple uniaxial tensile tests.

302
 The Electronic Journal of Mathematics and Technology, Volume 7, Number 4, ISSN 1933-2823

Figure 5: Initial mesh (top) and solution as a deformed mesh (bottom).

For a 2D problem the von Mises stress is

q
2 2 2
σvm = σ11 − σ11 σ22 + σ22 + 3σ12 . (11)

The next piece of code shows how to obtain and plot this value. The result is plotted in Figure 6.

//Obtaining the von Mises stress

fespace Wh(Th,P1); //We can define a new finite elements space
Wh Sigmavm;

//Stress tensor (since it is symmetric, it is enough with 3 elements)

macro Sigma(u,v) [2*mu*dx(u)+lambda*(dx(u)+dy(v)),
2*mu*dy(v)+lambda*(dx(u)+dy(v)),mu*(dy(u)+dx(v))]//EOM

//Von Mises stress

Sigmavm = sqrt(Sigma(u,v)[0]*Sigma(u,v)[0]-Sigma(u,v)[0]*Sigma(u,v)[1]
+Sigma(u,v)[1]*Sigma(u,v)[1]+3*Sigma(u,v)[2]*Sigma(u,v)[2]);
plot(Sigmavm); //Ploting the result

real Sigmavmmax = Sigmavm[].max; //Max. von Mises stress

cout << " Max von Mises Stress = " << Sigmavmmax << endl;

3.5 Skin effect in AC power electromagnetics

Next example shows the skin effect when AC current is carried by a wire of copper with circular cross
section. The propagation of plane electromagnetics waves may be modeled by a complex Helmholtz’s

303
 The Electronic Journal of Mathematics and Technology, Volume 7, Number 4, ISSN 1933-2823

Figure 6: Von Mises stress. We can see that the maximum stress concentrates on the lower left corner.

equation for the electric field, namely,

(  
1
−∇ · µ ∇Ec + (iωσ − ω 2 ε)Ec = 0 in Ω = {(x, y) ∈ R2 : x2 + y 2 < 0.1}
(12)
Ec = 1/σ on ∂Ω,
where Ec = Ec (x, y) is the complex electric field, σ is the conductivity, µ is magnetic permeability,
ε is the coefficient of dielectricity, ω is angular frequency, and i stands for the pure imaginary number
(also denoted by j in Electricity). The boundary condition indicates that, due to induction, the current
density in the interior of the conductor is smaller than at the outer surface where it is set to 1.
A weak solution of (12) is a complex function Ec = u + 1/σ, with u ∈ H01 (Ω; C) , such that u
satisfies
Z Z Z
1 2
∇u · ∇v dxdy + (iωσ − ω ε)u · v dxdy = −(iωσ − ω 2 ε)1/σv dxdy ∀v ∈ H01 (Ω; C) .
Ω µ Ω Ω

The FreeFem++ code for this problem is the following.

// physical parameters
real sigma=57*1e+6;
real mu=4*pi*1e-7;
real omega=2*pi*50;
//real omega=2*pi*1; a lower value for angular frequency
real epsilon=8.8*1e-12;
// domain and mesh
border C(t=0,2*pi){x=0.1*cos(t);y=0.1*sin(t);};
mesh Th=buildmesh(C(50));
// space of complex finite elements P1
fespace Vh(Th,P1);
Vh<complex> u,v;
// solving the variational formulation
solve electro(u,v)=int2d(Th)((1/mu)*(dx(u)*dx(v)+dy(u)*dy(v))
+(omega*sigma*1i-omega*omega*epsilon)*u*v)
+on(C,u=1/sigma);
plot(u,value=true,fill=true,ShowAxes=0,ColorScheme = 1,ps="fig_electro.eps");

304
 The Electronic Journal of Mathematics and Technology, Volume 7, Number 4, ISSN 1933-2823

Figure 7 displays the solution obtained with this code. It is observed that the largest values of the
electric field are concentrated on the boundary. This is the so-called skin effect. We suggest the reader
to play with different values of ω in the code to show how decreasing the frequency of the alternating
current results in a decreasing of the skin effect.

Figure 7: Picture of the modulus of Ec , solution of problem (12). It is observed that the largest values
of this solution are concentrated on the boundary. This is the so-called skin effect.

3.6 Reynolds’ equation in hydrodynamic lubrication

Consider the problem of lubricating the friction between a fixed rigid body presenting some abrupt
edges, and a regular surface in movement by using an incompressible fluid in the separating region.
This type of problems appears in real applications such as in feed-box or shaft-bearing systems. This
problem may be modeled by the following Reynolds equation in lubrication theory (see [7])
(    
U0 h h3 V0 h h3
2
− 12µ px + 2 − 12µ py = 0 in Ω
x y (13)
p = p0 on ∂Ω,

where Ω = [0, 1]2 , p = p (x, y) is the fluid pressure, µ is the fluid viscosity, h = h (x, y) represents
the distance between the surface and the solid, and (U0 , V0 ) is the velocity of the surface.
A weak solution of (13) is a function p = u + p0 , with u ∈ H01 (Ω) such that
h3
Z Z
h
∇u · ∇v dxdy = (U0 , V0 ) · ∇v dxdy ∀v ∈ H01 (Ω) .
Ω 12µ Ω 2
We take a discontinuous h to model an abrupt edge on the body, e.g.,

0.2, 0 ≤ x ≤ 0.5
h (x, y) =
0.8, 0.5 < x ≤ 1.
Since in the left-side term of the PDE we compute the partial derivative of h with respect to x,
we get a Dirac mass type right-hand side term for the equation. This implies that the solution is not

305
 The Electronic Journal of Mathematics and Technology, Volume 7, Number 4, ISSN 1933-2823

differentiable at x = 0.5. Consequently, in order to properly capture this singularity is very convenient
to refine the mesh in this region. We use the FreeFem++ function adaptmesh for this purpose.

real mu=0.391; // viscosity

mesh Th=square(100,100); // surface
plot(Th);
func h=0.2+0.6*(x>0.5); // separation between surface and solid
real U0=10, V0=5; // surface velocity
fespace Vh(Th,P1); // P1 elements
Vh p,v; // p=pressure
// variational formulation
problem lubrica(p,v)=int2d(Th)(((hˆ3)/(12*mu))*(dx(p)*dx(v)+dy(p)*dy(v)))
-int2d(Th)((h/2)*(U0*dx(v)+V0*dy(v)))
+on(1,2,3,4,p=1);
lubrica;
plot(p);
Th=adaptmesh(Th,p); // refine mesh according to singularities of p
lubrica;
plot(Th,p,value=true,fill=true,ColorScheme = 1,ps="fig_lubrica.eps");

Figure 8 shows the pictures of both the refined mesh and the solution of problem (13). It is
observed that the largest values of p are located at the same place where there is the discontinuity of
the function h, i.e. x = 0.5, which models the separation between the surface and the solid. Also,
Figure 8 seems to indicated that p is not differentiable at x = 0.5 in agreement with the theoretical
results in [7].

Figure 8: Picture of pressure p(x, y), solution of (13).

Acknowledgements. The authors express their gratitude to the FreeFem++ team in Paris IV for
providing us with this wonderful free software. Also, the second author is very grateful to Professor

306
 The Electronic Journal of Mathematics and Technology, Volume 7, Number 4, ISSN 1933-2823

Miguel Sanz Alix (University of Valencia) who introduced him variational methods for PDE, and
to José Antonio Vallejo for his kind invitation to the Fourth Summer School UVEG-UASLP held at
Universidad Autónoma of San Luı́s Potosı́ (Mexico) where a part of this paper was written. This
work was partially supported by project 08720/PI/08 from Fundación Séneca (Agencia de Ciencia y
Tecnologı́a de la Región de Murcia, Spain, II PCTRM 2007-10).

References
[1] Babuška, I., Error-bounds for finite element method. Numer. Math. (1970/1971) 16, 322-333.

[2] Evans, L. C., Partial Differential Equations, 1998 (Graduate Studies in Mathematics 19 American
Mathematical Society).

[3] F. Hecht, A. Le Hyaric, O. Pironneau, K. Ohtsuka, Tutorial FreeFem++ . Available at

http://www.freefem.org/

[4] Periago, F., A first step towards variational methods in engineering, Int. J. Math. Educ. Sci. Tech-
nol. (2003) 34, no. 4, 549-559

[5] Periago, F., Modelización matemática y simulación numérica en Ingenierı́a, Notes of a Master
Course (in Spanish). Available at the web site http://filemon.upct.es/˜fperiago/

[6] Raviart, P. A., Thomas, J. M., Introduction à l’Analyse Numérique des Équations aux Dérivées
Partielles, 1988 (Masson).

[7] Tello, J. I., Regularity of solutions to a lubrication problem with discontinuous separation data.
Nonlinear Anal. 53 (2003), no. 7-8, 1167-1177.

[8] Xia, X., Computing multiple integrals by MatLab, The Electronic Journal of Mathematics and
Technology, (2012) 6, n 2, 159-174.

307

View publication stats