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.
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.
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) = \{f(x) \in H^1(\Omega)\},
V_0(\Omega) = \{f(x) \in H^1(\Omega); f(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.
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.
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} }
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.
To solve numerically given problem based on the weak form (5) we have to go through 5 steps: