Notes on solving PDEs by the Finite Element Method¶

The Finite Element Method (FEM) is the numerical method for solving Partial Differential Equations (PDEs). FEM was developed in the middle of XX. century and now it is widely used in different areas of science and engineering, including mechanical and structural design, biomedicine, electrical and power design, fluid dynamics and other. FEM is based on a very elegant mathematical theory of weak solution of PDEs. In this section we will briefly discuss basic ideas underlying FEM.

Strong form of Poisson’s equation and its integration¶

Let us start our discussion about FEM with the strong form of Poisson’s equation

(1)¶ $\Delta T = f(x), \quad x \in \Omega,$

(2)¶ $T = u(x), \quad x \in \Gamma_D,$

(3)¶ $\nabla T \cdot \mathbf{n} = g(x), \quad x \in \Gamma_N,$

where $\Omega \subset \mathbb{R}^n$ is the solution domain with the boundary $\partial \Omega$ , $\Gamma_D$ is the part of the boundary where Dirichlet boundary conditions are given, $\Gamma_N$ is the part of the boundary where Neumann boundary conditions are given, $T(x)$ is the unknown function to be found, $f(x), u(x), g(x)$ are known functions.

FEM is based on a weak formulation. The weak form of the equation (1) is

$\int\limits_{\Omega} (\Delta T - f) \cdot s \, \mathrm{d}\Omega = 0,$

where $s$ is a test function. Integrating this equation by parts

$0 = \int\limits_{\Omega} (\Delta T - f) \cdot s \, \mathrm{d}\Omega = \int\limits_{\Omega} \nabla \cdot (\nabla T) \cdot s \, \mathrm{d}\Omega - \int_{\Omega} f \cdot s \, \mathrm{d}\Omega =$

$= - \int\limits_{\Omega} \nabla T \cdot \nabla s \, \mathrm{d}\Omega + \int\limits_{\Omega} \nabla \cdot (\nabla T \cdot s) \, \mathrm{d}\Omega - \int\limits_{\Omega} f \cdot s \, \mathrm{d}\Omega$

and applying Gauss theorem we obtain:

$0 = - \int\limits_{\Omega} \nabla T \cdot \nabla s \, \mathrm{d}\Omega + \int\limits_{\Gamma_D \cup \Gamma_N} \!\!\!\! s \cdot (\nabla T \cdot \mathbf{n}) \, \mathrm{d}\Gamma - \int\limits_{\Omega} f \cdot s \, \mathrm{d}\Omega$

or

$\int\limits_{\Omega} \nabla T \cdot \nabla s \, \mathrm{d}\Omega = \int\limits_{\Gamma_D \cup \Gamma_N} \!\!\!\! s \cdot (\nabla T \cdot \mathbf{n}) \, \mathrm{d}\Gamma - \int\limits_{\Omega} f \cdot s \, \mathrm{d}\Omega.$

The surface integral term can be split into two integrals, one over the Dirichlet part of the surface and second over the Neumann part

(4)¶ $\int\limits_{\Omega} \nabla T \cdot \nabla s \, \mathrm{d}\Omega = \int\limits_{\Gamma_D} s \cdot (\nabla T \cdot \mathbf{n}) \, \mathrm{d}\Gamma + \int\limits_{\Gamma_N} s \cdot (\nabla T \cdot \mathbf{n}) \, \mathrm{d}\Gamma - \int\limits_{\Omega} f \cdot s \, \mathrm{d}\Omega.$

The equation (4) is the initial weak form of the Poisson’s problem (1)–(3). But we can not work with it without applying the boundary conditions. So it is time to talk about the boundary conditions.

Dirichlet Boundary Conditions¶

On the Dirichlet part of the surface we have two restrictions. One is the Dirichlet boundary conditions $T(x) = u(x)$ as they are, and the second is the integral term over $\Gamma_D$ in equation (4). To be consistent we have to use only the Dirichlet conditions and avoid the integral term. To implement this we can take the function $T \in V(\Omega)$ and the test function $s \in V_0(\Omega)$ , where

$V(\Omega) = \{v(x) \in H^1(\Omega)\},$

$V_0(\Omega) = \{v(x) \in H^1(\Omega); v(x) = 0, x \in \Gamma_D\}.$

In other words the unknown function $T$ must be continuous together with its gradient in the domain. In contrast the test function $s$ must be also continuous together with its gradient in the domain but it should be zero on the surface $\Gamma_D$ .

With this requirement the integral term over Dirichlet part of the surface is vanishing and the weak form of the Poisson equation for $T \in V(\Omega)$ and $s \in V_0(\Omega)$ becomes

$\int\limits_{\Omega} \nabla T \cdot \nabla s \, \mathrm{d}\Omega = \int\limits_{\Gamma_N} s \cdot (\nabla T \cdot \mathbf{n}) \, \mathrm{d}\Gamma - \int\limits_{\Omega} f \cdot s \, \mathrm{d}\Omega, T(x) = u(x), \quad x \in \Gamma_D.$

That is why Dirichlet conditions in FEM terminology are called Essential Boundary Conditions. These conditions are not a part of the weak form and they are used as they are.

Neumann Boundary Conditions¶

The Neumann boundary conditions correspond to the known flux $g(x) = \nabla T \cdot \mathbf{n}$ . The integral term over the Neumann surface in the equation (4) contains exactly the same flux. So we can use the known function $g(x)$ in the integral term:

$\int\limits_{\Omega} \nabla T \cdot \nabla s \, \mathrm{d}\Omega = \int\limits_{\Gamma_N} g \cdot s \, \mathrm{d}\Gamma - \int\limits_{\Omega} f \cdot s \, \mathrm{d}\Omega,$

where test function $s$ also belongs to the space $V_0$ .

That is why Neumann conditions in FEM terminology are called Natural Boundary Conditions. These conditions are a part of weak form terms.

The weak form of the Poisson’s equation¶

Now we can write the resulting weak form for the Poisson’s problem (1)–(3). For any test function $s \in V_0(\Omega)$ find $T \in V(\Omega)$ such that

(5)¶ $\boxed { \begin{split} \int\limits_{\Omega} \nabla T \cdot \nabla s \, \mathrm{d}\Omega & = \int\limits_{\Gamma_N} g \cdot s \, \mathrm{d}\Gamma - \int\limits_{\Omega} f \cdot s \, \mathrm{d}\Omega, \quad \mbox{and}\\ T(x) & = u(x), \quad x \in \Gamma_D. \end{split} }$

Discussion of discretization and meshing¶

It is planned to have an example of the discretization based on the Poisson’s equation weak form (5). For now, please refer to the wikipedia page Finite Element Method for a basic description of the disretization and meshing.

Numerical solution of the problem¶

To solve numerically given problem based on the weak form (5) we have to go through 5 steps:

Define geometry of the domain $\Omega$ and surfaces $\Gamma_D$ and $\Gamma_N$ .
Define the known functions $f$ , $u$ and $g$ .
Define the unknown function $T$ and the test functions $s$ .
Define essential boundary conditions (Dirichlet conditions) $T(x) = u(x), x \in \Gamma_D.$
Define equation and natural boundary conditions (Neumann conditions) as the set of all integral terms $\int\limits_{\Omega} \nabla T \cdot \nabla s \, \mathrm{d}\Omega$ , $\int\limits_{\Gamma_N} g \cdot s \, \mathrm{d}\Gamma$ , $\int\limits_{\Omega} f \cdot s \, \mathrm{d}\Omega$ .