FEM (Finite Element Method)

FEM (Finite Element Method) is a special kind of numerical analysis technique developed mainly to solve various Partial Differential Equations. It splits a complicated region into many small and simple pieces, and it turns the equation into a matrix equation that a computer can solve.

I will put more description later.. for now refer to the following examples and try to get some intuitive understanding.

The sections below give the background that those examples assume. First, we'll look at the problem FEM solves and the steps every FEM analysis follows. Then we'll put numbers into the Simple Spring Model. Finally, we'll apply the same steps to a small differential equation, and you'll see that it ends in the same kind of matrix equation.

What problem does FEM solve ?

Many engineering questions are written as a partial differential equation over a region. Examples are the temperature inside an engine block, the stress in a bridge and the electric field around an antenna. For a simple shape, you can sometimes solve such an equation by hand. For a real shape with real materials and real boundaries, a closed form solution almost never exists.

FEM handles this by splitting the region into many small pieces called elements. The elements meet at points called nodes. Inside each element, a simple function approximates the unknown quantity. This function is usually a low order polynomial, and its values at the element's nodes fix it completely. So the unknown function over the whole region is replaced by a finite list of nodal values.

This is where the name comes from. The exact problem has an infinite number of unknowns, one value at every point. The FEM problem has a finite number, one value per node. The differential equation then becomes a system of linear algebraic equations for those nodal values. In structural problems it is written as K u = F. Here K is the stiffness matrix, u is the vector of nodal displacements, and F is the vector of nodal forces.

FEM replaces a function with a list of numbers : the unknowns are the values at the nodes, and a simple function inside each element fills in the values between them.
The result is always a matrix equation : whatever the physics, the analysis ends in K u = F. So solving an FEM problem means building and solving a large linear system.
Smaller elements give a better approximation : refining the mesh adds nodes, and the approximate solution comes closer to the exact one. The cost is a larger matrix.

What steps does every FEM analysis follow ?

Every FEM analysis follows the same sequence, from a two spring toy model to a full crash simulation. Commercial software hides most of these steps, but it runs all of them. If you know the sequence, you can follow what a software package asks for and what it reports. The table below lists the steps, and what each one means for the spring examples linked above.

Step	What happens	In the spring examples
1. Discretize	Split the region into elements that meet at nodes. The result is called the mesh.	Each spring is one element. Each end point or cart is one node.
2. Element equations	For one element, write the relation between the forces at its nodes and the displacements of its nodes.	A spring with stiffness k between nodes i and j gives the 2 x 2 element matrix k [1 -1; -1 1].
3. Assemble	Add every element matrix into the global matrix, at the rows and columns of that element's nodes.	A diagonal term is the sum of the stiffness of the springs attached to that node. An off-diagonal term is minus the stiffness of the springs joining the two nodes.
4. Boundary conditions	Set the nodal values that are known, and remove them from the unknowns.	A node fixed to the wall has zero displacement.
5. Solve	Solve K u = F for the unknown nodal values.	This gives the displacement of each free node.
6. Post-process	Compute derived quantities from the nodal values.	Spring elongation, spring force and the reaction force at the wall.

Step 3 is the rule that the More Complicated Spring Model example derives entry by entry. The same assembly rule works for any mesh, so a computer can apply it without knowing anything about the physics. Step 4 is not optional. Before the boundary conditions are applied, K is singular, and the next section shows why.

Assembly is addition : each element matrix goes into the rows and columns of its own nodes. Where two elements share a node, their entries are summed.
The stiffness matrix is symmetric and sparse : K_ij = K_ji, and K_ij is zero unless nodes i and j belong to the same element. A large mesh therefore gives a large matrix with mostly zero entries.
The truss adds directions, not new steps : in the Truss example each node moves in two directions, so each node has two rows in K. The six steps stay the same.

How does the two spring model become K u = F ?

Let's put numbers into the Simple Spring Model, so that you can follow every step with a calculator. The model and its node numbering are the same as on that page. Node 1 is fixed to a wall, and a force pulls Node 3 away from the wall.

The diagram below shows the model with the numbers used in this section. Spring 1 with stiffness k₁ joins Node 1 and Node 2. Spring 2 with stiffness k₂ joins Node 2 and Node 3. The displacements u₁, u₂ and u₃ are measured along the axis of the springs, and F₃ is the force applied at Node 3.

Figure 1. Two springs in series with Node 1 fixed. Node 2 belongs to both elements, so its row in K collects both stiffness values.

Each spring gives a 2 x 2 element matrix. After assembly, the 3 x 3 equation K u = F reads as follows. F₁ is the unknown reaction force at the wall. F₂ = 0, because nothing pushes Node 2 from outside, and F₃ = 10 N.

K	u₁	u₂	u₃	F
Row 1	k₁ = 100	-k₁ = -100	0	F₁
Row 2	-k₁ = -100	k₁ + k₂ = 300	-k₂ = -200	0
Row 3	0	-k₂ = -200	k₂ = 200	10

This matrix has determinant 0, so the equation cannot be solved as it stands. The reason is physical. Without a fixed node, the whole chain can slide along the axis without stretching any spring. That rigid motion, u₁ = u₂ = u₃, creates no spring force, so K times the vector [1, 1, 1] is zero. You can see the same fact in the table, because every row of K sums to 0.

The wall removes that freedom. Because u₁ = 0, you delete row 1 and the u₁ column. The remaining two equations are 300 u₂ - 200 u₃ = 0 and -200 u₂ + 200 u₃ = 10, and this 2 x 2 system has a unique solution.

The solution is u₂ = 0.1 m and u₃ = 0.15 m. So Spring 1 stretches by 0.1 m and carries 100 x 0.1 = 10 N. Spring 2 stretches by 0.15 - 0.1 = 0.05 m and carries 200 x 0.05 = 10 N. Both springs carry the full 10 N, as they must in a series chain. The deleted row 1 now gives the reaction at the wall: F₁ = -100 x u₂ = -10 N. It is equal and opposite to the applied force.

Rigid body motion makes K singular : a structure that can move without deforming gives K a zero eigenvalue. The boundary conditions remove that motion.
Deleted rows are not wasted : after the solve, the rows of the fixed nodes give the reaction forces.
Check the answer with physics : in a series chain every spring carries the same force, and here both springs carry 10 N.

How does a differential equation give the same matrix ?

Springs are an easy first example, but FEM was built for differential equations. Let's take the simplest one and follow the same steps. You'll see that it produces a matrix with exactly the same pattern as the spring chain.

The equation is -u''(x) = 1 on the interval from 0 to 1, with u(0) = 0 and u(1) = 0. One physical example is the sag of a taut string under a uniform load. Its exact solution is u(x) = x(1 - x)/2, so you can check the FEM answer against it.

Split the interval into 4 elements of length h = 1/4. The nodes are at x = 0, 0.25, 0.5, 0.75 and 1. Inside each element, u is a straight line between its two nodal values. The Galerkin method is the standard way to derive element equations from a differential equation. For this equation it gives each element the matrix (1/h) [1 -1; -1 1], which is the spring matrix with k = 1/h = 4. The uniform load puts h/2 on each end node of each element, so each interior node receives h = 0.25.

After assembly, the two boundary nodes are removed, exactly like the wall node above. Three unknowns remain, and the system reads 8 u(0.25) - 4 u(0.5) = 0.25, -4 u(0.25) + 8 u(0.5) - 4 u(0.75) = 0.25 and -4 u(0.5) + 8 u(0.75) = 0.25.

The solution is u(0.25) = 0.09375, u(0.5) = 0.125 and u(0.75) = 0.09375. The exact solution gives the same three values, for example 0.5 x 0.5 / 2 = 0.125. For this equation in one dimension, linear elements give exact values at the nodes. Between the nodes, the FEM answer is a straight line while the exact answer is a parabola. So the error sits inside the elements, and its largest value is h²/8 = 0.0078 at the middle of each element.

Different physics can give the same matrix : a chain of equal springs and this differential equation both give the tridiagonal pattern of 2 on the diagonal and -1 beside it. Heat conduction along a rod gives it too.
Boundary conditions play the role of the wall : u(0) = 0 and u(1) = 0 remove the rows of the end nodes, exactly as u₁ = 0 did for the springs.
Mesh refinement controls the error : the error inside an element is proportional to h². Halving the element length cuts it by a factor of about 4.