FEM, demathtified


The Finite Element Method (FEM) is a powerful tool for modeling fields, flows, structures, etc. FEM is used in many disciplines, but it is generally associated with engineering and physics, where phenomena are commonly described by sets of differential equations. In our case, we are dealing with oceanography and ecosystem models, but the principles remain largely the same.

There are a number of introductions to FEM, including dozens of books, papers, etc. Unfortunately for my purposes, most of them dive immediately into mathematics, without any motivation or conceptual background. If your math courses (like mine) did not extend past calculus and linear algebra, you may find this approach to be rather mathtifying (:-).

This page was written in an effort to understand the basic ideas underlying FEM, using as little mathematics as possible. Because I don't understand the underlying math and don't care about the implementation details, my explanations are likely to be simplistic or worse.

However, they may still be able to provide a useful conceptual introduction to the things that FEM software can do and the ways it needs to be used. My intent is to provide something like Lindley Darden's "grey box" descriptions of biological mechanisms, coupled with the "useful lies" of Kurt Vonnegut's Bokononism. So, YMMV... That said, let's give it a try...

Note: A gentle introduction to the Finite Element Method, by Francisco-Javier Sayas, contains some of the best introductory material (e.g., helpful text, useful images) I've encountered. If you have a solid background in mathematics and physics, Dr. Sayas's tutorial may be enough to let you get started on FEM. With his kind permission, I am using his example and some of his images below, while skipping over most of the details his notes provide. Any mistakes you find are almost certainly mine; please tell me about them and I'll see what I can do...


Wikipedia's math pages seem to be written by and for mathematicians; this one is a typical mixture of English and jargon:

In mathematics, the finite element method (FEM) is a numerical technique for finding approximate solutions to boundary value problems for differential equations. It uses variational methods (the calculus of variations) to minimize an error function and produce a stable solution. Analogous to the idea that connecting many tiny straight lines can approximate a larger circle, FEM encompasses all the methods for connecting many simple element equations over many small subdomains, named finite elements, to approximate a more complex equation over a larger domain. ...

-- Finite element method (WP)

The last sentence, however, captures the essence of FEM: "connecting many simple element equations over many small subdomains, named finite elements, to approximate a more complex equation over a larger domain". Let's say we have a continuous function f that describes a 3-dimensional surface, e.g.: z = f(x, y).

Because x, y, and z are all real numbers, the number of points on the surface is (uncountably) infinite. So, calculating the values for all of these points is impractical, at least for mortals:

God does not care about our mathematical difficulties. He integrates empirically.

-- Albert Einstein

So, instead of performing an infinite number of calculations, we cheat. Specifically, we create a piecewise approximation of the surface, using small (but still finite) 2-dimensional elements (e.g., rectangles, triangles). To get a good approximation, we may have to use lots (e.g., millions) of elements. However, this sort of calculation is quite feasible, even for personal computers.

Boundary Value Problems

The example we're using (borrowed from Dr. Sayas's tutorial) is a "boundary value problem", where we are trying to calculate the behavior of a differential equation within a region whose boundaries impose certain constraints, i.e.:

... a boundary value problem is a differential equation together with a set of additional restraints, called the boundary conditions. A solution to a boundary value problem is a solution to the differential equation which also satisfies the boundary conditions. ...

-- Boundary value problem (WP)

Although the definition above uses the term "differential equation" multiple times, there is actually no requirement that the equation be differential. Any continuous function (and even some discontinuous functions) can be modeled using FEM.

However, in practice most FEM applications involve differential equations. Aside from the fact that many scientific problems use them, there are significant performance advantages to using FEM with differential equations:

  • The FEM software may be able to traverse the boundary and generate the interior vertices, rather than exploring the entire interior of the domain.

  • Differential equations provide useful information about the function's derivative (i.e., slope) which can help the FEM software in exploring the domain's interior.

Differential Equations

Wikipedia's definition is probably entirely correct, but a bit dry for my taste:

A differential equation is a mathematical equation that relates some function of one or more variables with its derivatives. Differential equations arise whenever a deterministic relation involving some continuously varying quantities (modeled by functions) and their rates of change in space and/or time (expressed as derivatives) is known or postulated. ...

-- Differential equation (WP)

So, let's just say that a differential equation describes how something changes over time, based on its previous value and (uaually) other factors. For example, the speed of a falling raindrop will increase (due to gravity) until it reaches terminal velocity (due to wind resistance).

Boundary Conditions

There are several types of possible boundary conditions, but these two are by far the most frequently used:

A given domain may have any combination of Dirichlet and Neumann boundaries. In this example, we will have one of each.

Approximating the Boundary

The problem domain's actual boundary may be composed of straight sections, curved sections, or both. If there are curved sections, we'll need to approximate them, using sequences of straight sections. In the case of our example domain, this produces a five-sided, irregular polygon:

Sorry about the Greek letters; I'm afraid they come with the territory. The region inside the boundary (i.e., the physical domain) is called Ω (Omega). In our example, the boundary approximation is divided into two parts, called Γ (Gamma) D and N:

  • Gamma D's two sections define a Dirichlet (value-based) boundary condition.

  • Gamma N's three sections define a Neumann (slope-based) boundary condition.

Although our domain has an external, convex boundary, this is not a requirement. Boundaries may be concave, internal to the domain, or any other required shape. The only requirement is that they allow the domain to be approximated by one or more polygons.

Triangulating the Domain

Now that we have a polygonal approximation of the domain, we can triangulate it (i.e., divide it into finite, triangular elements), e.g.:


The image above only shows about two dozen triangles, but this is an extreme simplification: a typical finite element problem could use millions of elements. And, although we're using triangles, this is not a requirement; any polygon can be used as an element, as long as certain rules are followed:

  • Elements may not overlap.

  • Elements may only intersect at a common vertex or along a common full edge.

  • The triangulation must respect the partition of the boundary into Dirichlet and Neumann boundaries.

Calculating the z values

In the image above, we're looking down at the (approximated, triangulated) domain. So, we know nothing about the z value (i.e., height) of the equation at any point. For purposes of simplicity, let's assume that its initial value is always zero.

So, we know the initial x, y, and z values for each vertex. So far, so good. Now, we want to calculate the values for the following time step. To do this, we evaluate an arbitrary function, e.g.: z = f(x, y, ...). The function can use a variety of input values, eg:

  • the current time
  • x and y values (fixed) for the current vertex
  • other values (dynamic) for the current vertex
  • z values for neighboring vertices
  • slopes (e.g., ∂z/∂x) for related elements

The handling of slopes for related elements is a bit arbitrary: do we need to treat the element as an entity (with attributes) or can we simply store its attributes in the related vertices? For example, if we know the x, y, and z values for our own and neighboring vertices, can we simply calculate the relevant slopes? Never mind; implementation detail (:-).

Depending on the function and the values it uses, this approach can model either static or dynamic phenomena. In particular, if the new z value for a vertex calculation is based on the previous z value, we're edging toward differential equations.

Visualizing the Results

If we calculate a non-zero z value for a single vertex, we might see something like this:

One of our vertices has come up from the plane, in order to match the z value of the equation at that point in the domain. By adjusting each vertex, we can approximate any surface which our equations define. And, by repeating the calculations (e.g., for multiple times), we can model dynamic behavior.

This image has only a few dozen elements, so the approximation is very coarse. However, by using small enough elements, we can achieve an arbitrarily close approximation. In some FEM systems, all of the elements are the same size; in others, smaller elements are used in areas where z changes more rapidly. This is simply an optimization, however; the basic approach is unchanged.

Obeying Constraints

Because the elements are polygons, they are 2-D (i.e., flat, planar) objects. So, they have some useful properties:

  • The z value of any point on an edge can be determined from the values at its two neighboring vertices (i.e., endpoints).

  • The slope of the element (e.g., at an edge) can be determined from the values of any three vertices which do not fall on a line.

These properties allow us to ensure that every point along the boundary has the required value (if Dirichlet) or slope (if Neumann), to some specified degree of accuracy.


Inductive Process Modeling (IPM) defines entities (with attributes) and equations (with processes). We'd like to define some number of processes that handle certain situations, then attach those processes to sub-regions of the domain in question.

I am unsure whether this is even theoretically possible for current FEM software, let alone whether it has been implemented. Information on this issue would be most welcome!

This wiki page is maintained by Rich Morin, an independent consultant specializing in software design, development, and documentation. Please feel free to email comments, inquiries, suggestions, etc!

Topic revision: r10 - 06 Jan 2014, RichMorin
This site is powered by Foswiki Copyright © by the contributing authors. All material on this wiki is the property of the contributing authors.
Foswiki version v2.1.6, Release Foswiki-2.1.6, Plugin API version 2.4
Ideas, requests, problems regarding CFCL Wiki? Send us email