You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert
Here we will present a very brief introduction to partial differential equations (PDEs) and
to the finite element method (FEM). Perhaps the simplest PDE of all is the (steady-state, linear)
heat equation, also known as the Poisson equation. We will use this equation as a demonstrative
example of the method, and demonstrate how we go from the strong form of the equation, to
the weak form, and then finally to the discrete FE problem.
Strong form
The strong form of the heat equation may be written as:
where $u$ is the unknown temperature field, $\mathbf{q}$ is the heat flux, $f$ is an
internal heat source, and $\Omega$ is the domain on which the equation is defined. To
complete the problem we need to specify what happens at the domain boundary $\Gamma$.
This set of specifications is called boundary conditions. There are different types of
boundary conditions, where the most common ones are Dirichlet -- which means that the solution
$u$ is known at some part of the boundary, and Neumann -- which means that the gradient
of the solution, $\nabla u$ is known. Formally we write for our example
i.e. the temperature is prescribed to a known function $u^\mathrm{p}$ at the Dirichlet part
of the boundary, $\Gamma_\mathrm{D}$, and the heat flux is prescribed to $q^\mathrm{p}$ at
the Neumann part of the boundary, $\Gamma_\mathrm{N}$, where $\mathbf{n}$ describes the outward
pointing normal vector at the boundary.
We also need a constitutive equation which links the temperature field, $u$, to the heat
flux, $\mathbf{q}$. The simplest case is to use Fourier's law
$$\mathbf{q}(u) = -k \nabla u$$
where $k$ is the conductivity of the material. In general the conductivity can vary
throughout the domain as a function of the coordinate, i.e. $k = k(\mathbf{x})$, but for
simplicity we will consider only constant conductivity $k$.
Weak form
The solution to the equation above is usually calculated from the corresponding weak
form. By multiplying the equation with an arbitrary test function$\delta u$, integrating
over the domain and using partial integration we obtain the weak form. Now our problem
can be stated as:
Find $u \in \mathbb{U}$ s.t.
$$\int_\Omega \nabla \delta u \cdot (k \nabla u) \, \mathrm{d}\Omega =
\int_{\Gamma_\mathrm{N}} \delta u \, q^\mathrm{p} \, \mathrm{d}\Gamma +
\int_\Omega \delta u \, b \, \mathrm{d}\Omega \quad \forall \, \delta u \in \mathbb{T}$$
where $\mathbb{U}, \mathbb{T}$ are suitable function spaces with sufficiently regular
functions. Under very general assumptions it can be shown that the solution to the weak
form is identical to the solution to the strong form.
Finite Element approximation
Using the finite element method to solve partial differential equations is usually
preceded with the construction of a discretization of the domain $\Omega$ into a finite
set of elements or cells. We call this geometric discretization grid (or mesh)
and denote it with $\Omega_h$. In this example the corners of the triangles are called
nodes.
Next we introduce the finite element approximation $u_\mathrm{h} \approx u$ as a sum of N nodal
shape functions, where we denote each of these function by $\phi_i$ and the corresponding
nodal values$\hat{u}_i$. Note that shape functions are sometimes referred to as
base functions or trial functions, and instead of $\phi_i$ they are sometimes denoted $N_i$.
In this example we choose to approximate the test function in the same way. This approach is known
as the Galerkin finite element method. Formally we write the evaluation of our approximations
at a specific point $\mathbf{x}$ in our domain $\Omega$ as:
Since test and trial functions are usually chosen in such a way, that they build the basis of
some function space (basis as in basis of a vector space), sometimes are they are also called basis
functions. In the following the argument $\mathbf{x}$ is dropped to keep the notation compact.
We may now insert these approximations in the weak form, which results in
Since this equation must hold for arbitrary $\delta u_\mathrm{h}$, the equation must especially
hold for the specific choice that only one of the nodal values $\delta \hat{u}_i$ is fixed to 1 while
an all other coefficients are fixed to 0. Repeating this argument for all $i$ from 1 to N we obtain
N linear equations. This way the discrete problem can be written as a system of linear equations
where we call $\underline{\underline{K}}$ the (tangent) stiffness matrix, $\underline{\hat{u}}$ the solution
vector with the nodal values and $\underline{\hat{f}}$ the force vector. The specific naming is for
historical reasons, because the finite element method has its origins in mechanics. The elements
of $\underline{\underline{K}}$ and $\underline{\hat{f}}$ are given by
Finally we also need to take care of the Dirichlet boundary conditions. These are enforced by
setting the corresponding $\hat{u}_i$ to the prescribed values and eliminating the associated equations
from the system. Now, solving this equation system yields the nodal values and thus an
approximation to the true solution.
Notes on the implementation
In practice, the shape functions $\phi_i$ are only non-zero on parts of the domain $\Omega_\mathrm{h}$.
Thus, the integrals are evaluated on sub-domains, called elements or cells.
Each cell gives a contribution to the global stiffness matrix and force vector. The process
of constructing the system of equations is also called assembly. For clarification,
let us rewrite the formula for the stiffness matrix entries as follows:
This formulation underlines the element-centric perspective of finite element methods and
reflects how it is usually implemented in software.
Computing the element integrals by hand can become a tedious task. To avoid this issue we
approximate the element integrals with a technique called numerical integration. Skipping any
of the mathematical details, the basic idea is to evaluate the function under the integral at
specific points and weighting the evaluations accordingly, such that their sum approximates the
volume properly. A very nice feature of these techniques is, that under quite general
circumstances the formula is not just an approximation, but the exact evaluation of the integral.
To avoid the recomputation of the just mentioned evaluation positions of the integral for each
individual element, we perform a coordinate transformation onto a so-called reference element.
Formally we write
where J is the Jacobian of the coordinate transformation function. The computation of the
transformation, weights, positions and of the Jacobi determinant is handled by Ferrite.
On an intuitive level, and to explain the notation used in the implementation, we think of
$$\mathrm{d}\Omega \approx \, w \, \textrm{det}(J)$$
being the chosen approximation when changing from the integral to the finite summation.
For an example of the implementation to solve a heat problem with Ferrite check out [this
thoroughly commented example](@ref tutorial-heat-equation).
More details
We finally want to note that this quick introduction barely scratches the surface of the finite element
method. Also, we presented some things in a simplified way for the sake of keeping this article short
and concise. There is a large corpus of literature and online tutorials containing more details about
the finite element method. To give a few recommendations there is:
Introduction to the Finite Element Method by Niels Ottosen and Hans Petersson
The Finite Element Method for Elliptic Problems by Philippe Ciarlet
Finite Elements: Theory, Fast Solvers, and Applications in Elasticity Theory by Dietrich Braess
An Analysis of the Finite Element Method by Gilbert Strang and George Fix
Finite Element Procedures by Klaus-Jürgen Bathe
The Finite Element Method: Its Basis and Fundamentals by Olgierd Cecil Zienkiewicz, Robert Taylor and J.Z. Zhu
Higher-Order Finite Element Methods by Pavel Šolín, Karel Segeth and Ivo Doležel
This list is neither meant to be exhaustive, nor does the absence of a work mean that it is in any way
bad or not recommendable. The ordering of the articles also has no particular meaning.