Solving PDEs numerically: the finite element method (2)

Today, let’s continue with the finite element method and focus on how to implement it on a computer. How do we discretize our domain? How do we actually compute the stuff that’s in our variational formulation – efficiently, that is? How do we store and visualize our solution?

Last time, my article was quite theoretical and ‘mathematicalish’. Even if I tried not to dive too much into the details, I wanted to give some context and to present the different objects we are working with when we use the finite element method.

However, I only talked about mathematical stuff, so far. In particular, I mentioned space discretization, weak formulations and linear system solving. All of these things are very complex and interesting concepts that would require their own series of article. Here, let’s just analyze some of it and assume people have developed efficient ways to deal with the rest.

Cutting down space!

As we saw before, a basic idea of the finite element method is to consider a discrete equivalent of our continuous world that only contains a finite set of points, linked together by squares or triangles.

So, for example, instead of having a infinitely precise 1×1 square, we can have a much simpler geometry like this one:

A very simple discretization of the 1×1 square: we take 5 points and 4 triangles

You see that these 5 points and the 4 matching triangles sort of reconstruct our square (very crudely, I admit, but still!). Of course, the more points you use to discretize your domain, the closer you can be of the real thing, but the more memory you use up on your computer and the longer the computations will be.

There are usually two ways of discretizing a domain: either with a custom triangulation or with a regular grid. Depending on the type of geometry you are working with, each discretization type has its advantages and drawbacks:

  • regular grids are great for domains with a simple ‘homogeneous’ geometry, meaning a part of space that can be roughly approximated the same way everywhere; they are quite light to store in your memory since you can reconstruct the position of each points just by knowing its index and the discretization step
  • triangulations are better for ‘funky’ shapes: if you have a more complex geometry and/or if you need to add more precision to your approximation in a particular corner, for example; one of the most common triangulations is the Delaunay triangulation because it maps most geometries well and, in particular, it avoids flat triangles as much as possible

Like we saw last time, in 2D, regular grids are tilings made of squares and triangulations are tilings made of triangles. In 3D, they are made of cubes and tetrahedrons.

If you want to get a feel of how a square or a cube can be cut down into small pieces, there is a little program embedded in the last article to play with. Check it out!

Turning our equations into matrices

Why matrices?

In the last article, I mentioned weak formulations as a way of ‘averaging’ our problem and solving a simpler one that approximates our real equations well enough. From a mathematical point of view, and very grossly speaking, this means replacing our variables and their derivatives with integrals.

If you happen to be familiar with the notion of distribution, variational formulations actually rely on multiply our equation by a test function, then integrating the resulting expression over our domain. Thus, for example, the following equation:

Δu + k2u = f

could be transformed into something like:

-∫Ω ∇u . ∇v + k2Ω u.v = ∫Ω f.v

where v is the aforementioned test function. We usually denote a(.,.) the 2-variables function on the left side of the equation that depends on the solution u and v, and l(.) the 1-variable function on the right side that only depends on the problem’s parameters and v.

But as we’ve seen extensively, the main problem when solving physical problems on a computer is representing your data, because you try to work with infinite things: both continuous space and integrals mean an infinite number of points! So, now that we have taken care of space, we see that we also need to change these integrals into something a computer understands… such as a matrix.

Basically, just like we represent our continuous space by a finite set of points, we only consider our integrals on these points and therefore get tables of values, a.k.a. matrices or vectors, for each of them. We thus transform our mathematical expression filled with integrals into a linear system, easier to solve on a computer.

What matrices?

The resulting integrals and therefore matrices depend both on the initial physical equation and on the various boundary conditions you apply on your domain.

In the physical equations we focus on, there can be 3 types of terms:

  • things that involve the values of our solution
  • things that involve the values of the derivative of our solution
  • things that do not involve our solution but rather the rest of the functions and parameters specific to our physical situation

Each of these types is associated with a matrix or a vector, respectively: the mass matrix, the stiffness matrix and the right-hand side vector.

In the first post, I briefly talked about boundary conditions; the main idea is to “force” your solution to have certain properties at certain points of your domain to insure it is unique. There are 3 big types of boundary conditions:

  • Dirichlet conditions: you set the value of your solution
  • Neumann conditions: you set the value of a flux, so you impose the value of the gap between the value of your solution inside of your domain and outside of your domain at this point
  • Robin-Fourier conditions: they mix the two above and determine a flux that depends both on the value of the solution itself and the value of its derivative on the border of the domain

Testing it for yourself

This little program show you how these conditions influence the final linear system, i.e. what matrices we get in the end from our initial equation (here, we consider an inhomogeneous Helmholtz equation):

[code_embed src=”” height=700]

Computing the matrices

We have identified the matrices we need to compute to translate our initial equation into stuff that is understandable by a computer. Two questions remain, however:

  • how do you compute all these matrices efficiently?
  • how do you solve the final linear system quickly to obtain your solution?

Remember we said that we have one value for each matrix in our system, for each point in our discretized domain. If you decide to cut down space in a very thin way, to have a good approximation of the continuous domain, then you may be faced with a huge number of values. This means both creating your matrices and solving the linear system can be very long, if you do it with a naïve method.

The first idea that comes to mind is to simply go over each of the points in our domain and compute each coefficient of our matrices. However, this is a very inefficient method. Instead, the finite element method uses matrix assembly to make this computation quicker. The goal is to go over the triangles in our domain rather than the points and to compute their individual contributions, then sum them to get the global matrices.

To do so, you need to be able to match local and global indices for your points, because as you iterate over your triangles, you only get access to their three local vertices… but when you compute the coefficient, you want to put it at the right place in the global matrix, which is given by the global index of the vertex. This is why we often define a Loc2Glob function that holds this index mapping.

Take a look at this script to see how local and global indices relate to each other:

[code_embed src=””]

You can click on a triangle to zoom in and see the mapping between the global and local indices of its 3 vertices; click in a corner of the frame to zoom out.

The pseudo-code algorithm for computing the matrices is therefore:

A = 0; // initial null matrix
B = 0; // initial null vector
// iterate over the triangles
For p = 1:nb_triangles
   // iterate over the 3 vertices of the triangle
   For i = 1:3
      // global index of the i-th vertex
      I = Loc2Glob(p, i)
      // re-iterate over the 3 vertices
      For j = 1:3
         // global index of the j-th vertex
         J = Loc2Glob(p, j)
         // value of the left-hand side of the
         // equation at this point
         A[I,J] += a(.,.)
      // value of the right-hand side of the
      // equation at this point
      B[I] += l(.)

This script shows how the contribution of each triangle in our mesh is ‘assembled’ into our final mass matrix:

[code_embed src=”” height=600]

We do the same for the stiffness matrix and sum it with the mass matrix to have our left-hand side, then we do the same for the right-hand side.

The last step is the application of Dirichlet conditions if there are any in our problem. As we’ve seen last time, contrary to Neumann and Fourier-Robin conditions that are automatically integrated into our variational formulation, Dirichlet conditions must be treated on their own by reassigning the right value to our solution at specific locations on our mesh. To do so, for each discretization point of index i touched by the condition, we must:

  • set the value of the i-th coefficient of the right-hand side to the value of the Dirichlet condition
  • set the i-th row of the left-hand side matrix to zero
  • set the (i, i) coefficient of the matrix to one

After that, we get the linear system we have to solve.

To solve it efficiently, as mentioned in the previous article, you usually use specific libraries that rely on the BLAS routines to perform linear operations very quickly.

The end result?

After solving the linear system, you end up with a vector that holds as many values as there are points in your domain discretization. This the approximated discrete solution of the weakened problem – so an approximation for the approximated problem. Even if it seems like we have simplified a lot of things to suit us, results are actually quite accurate!

Checking how far off they are from the actual solution is the whole goal of finite element analysis which I’ll talk more about in the third and last article of this series.

Visualizing your solution – and visualizing it properly – is important because it helps you compare and interpret the results more easily. The finite element method is often used with 2D or 3D space. In this context, it is sometimes worth wondering how to correctly show the values of your solution: given that a screen is already (only) a 2D surface, how can we display so many pieces of information at once? Although there is no absolute right answer and it depends on the problem you are solving, there are a few common ways of plotting and displaying your solution:

  • colormaps: you display the values of your solution with a set of colors; this allows you to use the x-, y- and z-positions for the space coordinates and the additional channel of information provided by colors to easily show values
  • altitude and warping: if you are working on a 2D domain, sometimes you can use the 3rd coordinate (the z axis) to show the value of your solution at each (x, y) point: this is like an elevation map that takes a 2D set of altitudes and ‘throws’ it into a higher dimension-space

Here are a couple of images that use both these techniques:

Visualization with a colormap: each color corresponds to a value of the solution (see the range on the right)

Visualization with warping (+ colormap): in addition to color mapping, we translate the points along the z-axis according to the value of the solution at this location

These visualizations can be created thanks to Python libraries like matplotlib or software like Paraview.

What next?

In the last article of this series, I’ll present some additional notions and further developments of the finite element method. In particular, I’ll talk about finite element analysis, some methods for solving time-dependent problems and optimization ideas.

    1. My teacher’s website:
    2. Joseph Flaherty’s course on Finite Element Analysis (see ‘Course Notes’):
    3. Paraview’s website:
    4. I. Wikimedia Foundation, “Finite element method” (, December 2018. [Online; last access 09-January-2019].
    5. I. Wikimedia Foundation, “Boundary value problem” (, February 2018. [Online; last access 13-January-2019].
    6. I. Wikimedia Foundation, “Delaunay triangulation” (, January 2019. [Online; last access 10-February-2019].
    7. I. Wikimedia Foundation, “Helmholtz equation” (, February 2019. [Online; last access 10-February-2019].

Leave a Reply

Your email address will not be published.