Finite Element Method

The Finite Element Method is widely used by Electrical, Mechanical, and Civil Engineers to solve the partial differential equations that describe things like the stress on bridges and airplane wings, the flow of heat, the propagation of electromagnetic fields, and other such problems. If you don't know what a partial differential equation is, but you've had some calculus, then suffice it to say that we are trying to find the value of one or more variables (like stress, strain, temperature, electric field, magnetic field, voltage, current, etc.) at all points in a defined space. The laws of physics give us equations that define the behavior of these quanitites, and our job is to solve the equations for the variables. The equations usually involve multiple derivatives with respect to the three dimensions of space and perhaps time as well.

If you've had basic calculus and can remember what a derivative is, then you can follow this. I'll only use nasty equations to define the problem, and then I'll wave my hands through the rest, but hopefully in a way that lets you figure out what it's all about. The finite element method is a numerical method for generating an approximate solution to partial differential equations. For example, the electromagnetic laws that govern the behavior of an electrostatic field in any region of space, that does not contain current sources, may be expressed in the form of a differential equation called the Laplace equation:

del [ c(x,y,z) del V(x,y,z) ] = 0

where c(x,y,z) is the conductivity over the region of space (called the domain), V(x,y,z) is the voltage, and "del" is the vector gradient operator. Note that I simply wrote c(x,y,z) as c, and V(x,y,z) as V. This is (almost) as complex as the problem definition gets for a three-dimensional domain. For 2-D problems, simply eliminate one of the dimensions in the preceding equation.

Now, in what is called the "forward problem" we know c(x,y,z) and want to solve this equation for V(x,y,z). In general, the solution cannot be expressed as a simple analytic function (that means it can't be written down using a few simple polynomials or exponentials or trig functions, or whatever). In fact, we often cannot even express c(x,y,z) as a simple function, and have to make due with drawing the domain as a collection of irregularly shaped regions and saying what the conductivity is in each region. We will do the same thing with the approximate solution for V.

That is, we will break the domain into small regions and say that V has approximately the shape of some simple function over each subregion. You might expect that these simple shapes could be constant values of V over each subregion, but it turns out (for reasons that I won't bore you with) that they have to be at least as complicated as linear functions of (x,y,z) over each subregion. Don't panic, let me illustrate what that means. Let's consider a 2-D problem, where V is a function of x and y, and we have formed little triangles to define our subregions. The approximate solution for each little triangle then forms a plane, and the little planes all fit together at the edges to form an approximation of the solution for V(x,y), like this:

The little subregions are called "elements," and hence the name "Finite Element Method." So, how do we figure out how each little plane is oriented? Well, if you have ever had to fit a line to approximate a set of data points, then you know the basic idea already. You may even have heard of a procedure called "least squares" that defines a criterion to select the best line through these data points. This procedure chooses the line that minimizes the sum of the squared error at each data point, where the error is the difference between the data point and the line. Fitting our little planes is the 2-D version of fitting a line. Of course, we have to do it once for each element, and we have to ensure that the edges all fit together (again, for reasons that I won't bore you with). We have one major problem though: we don't have any data points to fit to! We just have that differential equation. We can overcome this difficulty by defining the error we want to minimize as how close our little plane-shaped approximation of V(x,y) evaluates to zero when we plug it into the left-hand side of that differential equation. The right hand side says that this should evaluate to zero for the correct solution, so we will assume that something that evaluates close to zero is close to the true solution (this assumption can be proven to be justified, but that would take us into a long and boring mathematical diversion). There are different ways of summing the error over each element, and this whole idea for solving a differential equation goes under the name "Methods of Weighted Residuals." When we use this idea over these little subregions, then we have a "Finite Element Method."

The significance of my work on adaptive non-conforming meshes is that it extends the application of the finite element method into an arena where it was previously difficult to apply. That is where the space over which we are searching for a solution is very complicated, containing highly irregular structures. A perfect example is the human body. There are a lot of different tissues with complex shapes and different properties that show up in the the equations. It can be very difficult to generate a finite element mesh over these domains using the traditional tools that were developed for man-made structures. In a nutshell, my solution is to create a fairly simple (but unorthodox) mesh over the domain that is able to adapt itself as the equation is solved. The unorthodox structure of my mesh creates some mathematical problems with the finite element method. This just makes the algorithm more complicated, but the end result is the ability to model these types of problems much more rapidly than if traditional tools are used. Here is a picture of the elements in a 3-D adapted mesh over a simple cube-shaped domain. The electric field in this example is zero just off the corner of the cube nearest you. The field forms a large peak everywhere on the surface of a sphere that cuts through the middle of this cube. You can see that the elements are smaller near the surface of this sphere, because smaller elements are required to accurately approximate the field wherever it is rapidly changing. This is because each element approximates the electric field as a linearly varying quantity within its region (it actually allows supralinear, but not quadratic, variation, if that means anything to you).

Note that this is a picture of where the elements are, NOT a picture of the electric field (which would look like the shell of a sphere cutting through this cube-shaped volume). The important thing is that the method figured out where to put smaller elements on its own. I didn't have to tell it where it was necessary to put small elements and where it could get away with larger ones. I only told it how accurate I wanted the solution overall. This is why the method is called "adaptive." This is important for bioengineering applications because there are a lot of different tissues in the human body (it is very "heterogeneous") and the boundaries between them are complex. This makes it pretty much impossible to know in advance where big and small elements should go for a given problem.

Here is an example cross-section through the thorax of a pig. The image on the left is colored according to tissue type, and was generated from a CT image (some call this a CAT scan). The black dot in the right ventricle is the cross-section of an internal defibrillation electrode. A model was built for this pig from a collection of such classified CT images taken at different slices. The model was solved adaptively for the electric field produced by a defibrillation pulse. The image on the right is a cross section through the model at the same level and shows the elements that were generated in this slice. Notice that the elements are small near the black electrode, and also towards the upper right, because that is the primary path of current to another electrode placed just under the skin of the pig (out of the plane of this slice).

By the way, notice that the 2-D example in the picture at the top uses triangular elements. The analog of a triangle in 3-D is a tetrahedra (a 3-D volume with 4 vertices and 4 faces), whereas in these pictures I am using hexahedra (8 vertices and 6 faces). These hexahedra are also parallelipipeds (the opposing sides are always parallel). This is because the models are derived from classified medical images, and it is much more convenient to use an element shape that corresponds to a collection of image voxels (each image pixel represents a discrete volume in the model, hence the term "voxels"). Why do I bother mentioning it? Notice that when the size of adjacent brick-shaped elements are different, a vertex of one element will fall on the center of the side of an adjacent element. This creates some mathematical difficulties, as described below.

In the 2-D example above, the vertices of adjacent elements always line up with each other. The presence of non-aligned (or "irregular") vertices introduces some mathematical complications because it becomes more difficult to make the approximate solution continuous across the element boundaries. Why? Remember that the solution over any triangle forms a plane. If the edge of one plane butts up against the edges of two other planes, then the three won't necessarily line up because of the extra vertex. This could create a gap in the approximate solution, as illustrated in red below. It can be proven that if this is allowed under the normal Finite Element procedures, the approximate solution will NOT approach the true solution as elements are made arbitrarily small.

The fix to this problem is to "constrain" the solution at that extra vertex so that it falls on the edge of the adjacent (larger) plane. This adds some complexity to the algorithm because it creates extra dependencies between the element solutions. Although this example uses triangular elements, the problem is the same using rectangular elements. With triangles, you can achieve gradual variations in element size while keeping the vertices aligned. With rectangles, you can't: if the adjacent rectangle is a different size, then one of the vertices won't line up. This is the price paid for achieving a nonuniform mesh density using only rectangular (2-D) or brick-shaped elements (3-D).

On this site you can find an interactive Windows program for my adaptive nonconforming mesh finite element solver. You can also find some code for one type of inverse problem. Looking at a slight modification of the differential equation, called the Poisson equation, that allows for current sources in the volume (rho):

del [ c(x,y,z) del V(x,y,z) ] = -rho(x,y,z)

an inverse problem is when we know c(x,y,z) everywhere, know V(x,y,z) at some locations (like, say, on the surface of the head), and want to know where the current sources are (representing neural activity in the brain).

Another type of inverse problem is where we know where the sources are (typically voltages or current sources on the boundary), and want to know what the conductivites are. That type of inverse problem goes by the name "impedance imaging." You won't find code for that one here. Not at present anyway. I will post a tutorial, including some code, on impedance imaging if and when I obtain funding to produce that tutorial (there is a grant proposal pending that includes that work).