simpleFEM is a tiny finite element code created for educational purposes by Roger Donaldson. The C++ source for simpleFEM can be downloaded here.

simpleFEM solves the following problem:

\begin{align} \left\{ \begin{array}{rl} -\Delta u = f,& x \in \Omega, \\ u = 0,& x \in \partial \Omega, \end{array} \right. \end{align}

where $\Omega \in \mathbb{R}^2$ is the pie-shaped domain

\begin{align} \Omega = \{(r,\theta) \, | \, 0 \leq r \leq 1, \, 0 \leq \theta \leq 3\pi/2 \}. \end{align}

simpleFEM calls the Triangle program to create a triangular mesh, assembles the finite element problem using piecewise-linear triangular elements, and solves the problem using a diagonally-preconditioned conjugate gradient solver.

The solution procedure for this specific case using piecewise-linear triangular elements is outlined in this talk.

The mesh and the solution for $f = \sin(x+y)$ are shown.

What simpleFEM does

simpleFEM illustrates each of the components finite element codes tend to have. In particular, it

  • reads in user-set parameters
  • subdivides the domain using a polyhedral mesh
  • assembles the problem as a linear system
  • solves the linear system
  • prints the solution

In addition, simpleFEM contains structures for

  • managing calls to each program component
  • storing and navigating through the mesh
  • solving the linear system and storing the solution

Of the tasks and structures required, the only substantial original code provided by simpleFEM is that which stores the mesh, assembles the linear system, and performs the linear solve. Here are some of the other components:

Parameter handling

Richard Wagner's ConfigFile program is used to read input parameters. ConfigFile supports an intuitive human-readable file format.

Mesh generation

Jonathan Shewchuk's Triangle program is used to triangulate the domain and generate the mesh. This fantastic code is written in C (compiles nearly everywhere) and has a remarkable array of features for producing high quality triangulations. Triangle can also refine existing meshes, generate convex hulls of points, and generate Voronoi diagrams of triangulations. Triangle runs quickly, and interacts with the outside world through several intuitive file formats. Be sure and send Jonathan a quick email if you like his program.

Linear solver

Strictly speaking, the linear solver is hand-written, but the algorithm is taken from Shewchuk's An Introduction to the Conjugate Gradient Method Without the Agonizing Pain [1]. If you don't know how the conjugate gradient method works, and can be patient with Shewchuk's conversational style, this paper is for you.

Numerical support

In addition to the cmath library, simpleFEM draws upon the Gnu Scientific Library. This is the only library not supplied with the simpleFEM code, but is available precompiled for many linux distributions, and is available through macports to OS X users.

What simpleFEM doesn't do

Most importantly, simpleFEM does not assemble the stiffness matrix using a change of coordinates, but computes the stiffness sub-matrices using the cotangent formula. This simplification is possible because simpleFEM uses piecewise-linear triangular finite elements, and because the diffusion coefficient is constant over the entire domain (ie, we solve the Poisson equation). Thus, the stiffness matrix elements on triangles, namely

\begin{align} K_{ij} = \int_T \nabla \phi_i^T \nabla \phi_j, \end{align}

can be computed exactly, without the need for a fancy quadrature rule, since its integrand is constant over each triangle.

This point deserves further elaboration. More will be said as this page is developed.

1. Jonathan Richard Shewchuk.An Introduction to the Conjugate Gradient Method Without the Agonizing Pain.
Unless otherwise stated, the content of this page is licensed under Creative Commons Attribution-ShareAlike 3.0 License