Subsections

• 2.3.1 Basic Concepts for a Galerkin Method

2.3 Finite Elements

The finite element method is a systematic approach to approximate the unknown exact solution of a partial differential equation based on basis functions and the projection of a given domain onto a consistent finite cell complex. The finite elements correspond to the $p$ -cells of the complex.

The origin of the finite element method is in the solid mechanics of rigid bodies [29]. Some constraints should be placed on the selection of the shape functions to guarantee the fact that with an arbitrary number of shape functions the exact solution is approximated best, and in the limit, the exact solution should be obtained.

A good approximation is obtained by a residual formulation where the residuum is formulated as the difference of the unknown exact solution and the calculated approximate solution. This residuum is weighted over the simulation domain and integrated with the requirement that the integral vanishes with a set of linearly independent weighting functions. The other possible mechanism is the variational formulation of the partial differential equation [28,68].

One of the main advantages of the finite element method is the possibility to adapt the basis functions to the eigenfunctions of the differential operators [28]. Thereby a high precision of this method can be obtained with a moderate number of mesh elements. Different areas of materials as well as anisotropic, inhomogeneous, and non-linear quantities can be treated. If no boundary conditions are declared, homogeneous/natural Neumann boundary conditions are implicitly given.

2.3.1 Basic Concepts for a Galerkin Method

In the following, a theoretical part of the finite element method is summarized, which is a special Galerkin method [28,29] based on the following construction of finite dimensional subspaces ${\ensuremath{\mathcal{V}}}_h$ :

• ${\ensuremath{\Omega}}$ is triangulated into a consistent cell complex
• ${\ensuremath{\mathcal{V}}}_h$ consists of piecewise polynomials
• ${\ensuremath{\mathcal{V}}}_h$ has a finite basis ${\ensuremath{\phi}}$ with local support

The main part of this section and the corresponding notation is based on [69,29]. Galerkin's method can be briefly explained as a general technique for the construction of solution approximations not represented by infinite space basis functions ${\ensuremath{\phi}}_i \in {\ensuremath{\mathcal{V}}}$ but instead by a finite dimensional space ${\ensuremath{\mathcal{V}}}_h \subset {\ensuremath{\mathcal{V}}}$ of variational problems. The examples presented in this work are modeled by the following general form:

$$\mathcal{L}(u) = f(\ensuremath{\mathbf{x}})$$ (2.22)

which is defined on an arbitrary dimensional, single connected domain ${\ensuremath{\Omega}}$ with boundary $\partial {\ensuremath{\Omega}}$ . Equation 2.22 is fulfilled by functions of class $C^{2}$

$$u = u(\ensuremath{\mathbf{x}}) \in {\ensuremath{\mathcal{V}}}$$ (2.23)

The notation $\mathcal{L}()$ being a linear spatial differential operator. Additionally, it is assumed that the domain ${\ensuremath{\Omega}}$ has a piecewise smooth boundary ${\ensuremath{\Gamma}}_i \in {\ensuremath{\Gamma}}(\forall i \neq j \rightarrow... ...amma}}_j = \emptyset, {\ensuremath{\Gamma}}= \bigcup_i {\ensuremath{\Gamma}}_i)$ . A general form of boundary conditions can thereby be specified by:

$$\ensuremath{\mathbf{A}}_i  u + \ensuremath{\mathbf{B}}_i \partia... ... = \ensuremath{\mathbf{g}}_i(u) \quad \textnormal{on}  {\ensuremath{\Gamma}}_i$$ (2.24)

where $\ensuremath{\mathbf{A}}_i$ and $\ensuremath{\mathbf{B}}_i$ are matrices consisting of functions sufficiently smooth on ${\ensuremath{\Gamma}}_i$ and $\ensuremath{\mathbf{g}}_i$ is a vector of continuous linear functionals. $\partial_n u$ denotes the outward normal derivative. A weak formulation of Equation 2.22 is obtained by using a test function $v$ :

$$\textnormal{Find}\; u \in {\ensuremath{\mathcal{V}}}: \; \langle ... ...u) \rangle = \langle v,f \rangle \quad \forall v \in C^1({\ensuremath{\Omega}})$$ (2.25)

The infinite-dimensional space ${\ensuremath{\mathcal{V}}}$ is then replaced by a sequence of finite-dimensional spaces ${\ensuremath{\mathcal{V}}}_h$ with $n$ linearly independent functions ${\ensuremath{\phi}}_i \in {\ensuremath{\mathcal{V}}}$ which span the space ${\ensuremath{\mathcal{V}}}_h$ , e.g.,

$${\ensuremath{\mathcal{V}}}_h = \left \{ v: v(\ensuremath{\mathbf{... ...ensuremath{\phi}}_i(\ensuremath{\mathbf{x}})\right \}, \quad s_i \in \mathbb{R}$$ (2.26)

The basis should build a complete function system to guarantee the approxicatiom for $n \rightarrow \infty$ . Any function $v_h \in {\ensuremath{\mathcal{V}}}_h$ is uniquely determined by a finite number of degrees of freedom, e.g., function values. This results in the following discrete variational problem for conforming methods ( ${\ensuremath{\mathcal{V}}}_h \subset {\ensuremath{\mathcal{V}}}$ ):

$$\textnormal{Find} \; u_h \in {\ensuremath{\mathcal{V}}}_h: \; \la... ...ngle = \langle v_h,f \rangle \quad \forall v_h \in {\ensuremath{\mathcal{V}}}_h$$ (2.27)

By selecting ${\ensuremath{\phi}}_i \in {\ensuremath{\mathcal{V}}}_h$ the following system is obtained:

$$\langle {\ensuremath{\phi}}_i,\mathcal{L}(u_h) \rangle = \langle {\ensuremath{\phi}}_i,f \rangle \quad \textnormal{for}  i=1, ..,n$$ (2.28)

The solution variable $u_h$ is then also expanded in terms of the basis $({\ensuremath{\phi}}_i)_{1 \leq i \leq n}$ of ${\ensuremath{\mathcal{V}}}_h$ .

$$u_h = \sum_{i=1}^n u_{i}  {\ensuremath{\phi}}_i$$ (2.29)

A graphical representation of $u_h$ is given in Figure 2.7.

Figure 2.7: Geometrical interpretation of the basis of the expanded solution variable $u_h$ for finite elements for a one-dimensional cell complex.

By the given selection of the test function and the expansion of the solution variable, it can be observed [69] that the scalar product in Equation 2.28 transforms the differential operator $\mathcal{L}()$ into the discrete operator

$$L: \mathbb{R}^n \rightarrow \mathbb{R}^n$$ (2.30)

The residuum is obtained by using the operator of Equation 2.30

$$R_i = L(u_1, ..,u_n) - \langle {\ensuremath{\phi}}_i , f \rangle$$ (2.31)

The next step for the discrete approximation of a given problem is the subdivision of the domain ${\ensuremath{\Omega}}$ according to cell complex properties (see Section 1.3 for details). The elements of the cell complex, the collection of cells, is then used as finite elements for the finite space ${\ensuremath{\mathcal{V}}}_h$ . The basis functions ${\ensuremath{\phi}}_i \in {\ensuremath{\mathcal{V}}}_h$ are defined on the global vertices $\ensuremath{\tau_{0}^{i}}$ , the 0-cells, only and are therefore called nodal basis functions. They can be expressed as:

$${\ensuremath{\phi}}_i(\ensuremath{\tau_{0}^{i}}) = \delta_{ik}, \quad i,k = 1, .., n$$ (2.32)

Based on this subdivision of the finite element space, the function $u_e$ can be expressed in local cell terms by

$$u_c = \sum_{i=1}^{n_p} u_i  {\ensuremath{\phi}}_i$$ (2.33)

where $i$ represents the local index of the element and $n_p$ the number of element vertices. Now the operator from Equation 2.30 can be defined locally for, e.g., a 3-simplex (tetrahedron):

$$L^c: \mathbb{R}^4 \rightarrow \mathbb{R}^4$$ (2.34)

and the local residuum is then defined as:

$$R_i = L^c(u_1, ..,u_n) - \langle {\ensuremath{\phi}}_i , f \rangle$$ (2.35)

To determine the operator given in Equation 2.35 for a particular cell type, the basic nodal functions have to be calculated, e.g., for a 3-simplex cell (tetrahedron):

$${\ensuremath{\phi}}_k(\ensuremath{\mathbf{x}}) = {\ensuremath{\phi}}_k(x_1, x_2, x_3) = \frac{1}{3V} (a_k x_1 + b_k x_2 + c_k x_3 +d_k)$$ (2.36)

where the coefficients $a_k, b_k, c_k, d_k$ are functions of the vertex coordinates and $V$ is the volume of the element. Next, the following integral has to be calculated:

$$M_{pq} = \int_{\ensuremath{\tau_{3}^{i}}} {\ensuremath{\phi}}_p(\... ...phi}}_q(\ensuremath{\mathbf{x}})  d {\ensuremath{\Omega}}, \quad p,q = 1,2,3,4$$ (2.37)

The key for an efficient practical realization is that global finite elements are defined as transformations of a reference element in a normalized coordinate system. Each point of the cell $(x_1, x_2, x_3) \in \ensuremath{\tau_{3}^{i}}$ can be expressed as a bijective function onto a reference point $(\xi, \eta, \zeta) \in \ensuremath{\tau_{\mathrm{r},3}^{i}}$ :

$$x_1 = x_1^1 + (x_1^2 - x_1^1)\xi + (x_1^3 - x_1^1)\eta + (x_1^4 - x_1^1)\zeta,$$ (2.38)

$$x_2 = x_2^1 + (x_2^2 - x_2^1)\xi + (x_2^3 - x_2^1)\eta + (x_2^4 - x_2^1)\zeta,$$ (2.39)

$$x_3 = x_1^1 + (x_3^2 - x_3^1)\xi + (x_3^3 - x_3^1)\eta + (x_3^4 - x_3^1)\zeta$$ (2.40)

The nodal basis functions on the reference element $\ensuremath{\tau_{\mathrm{r},3}^{i}}$ are:

$${\ensuremath{\phi}}_1^0(\xi, \eta, \zeta) = 1 - \xi - \eta - \zeta,$$ (2.41)

$${\ensuremath{\phi}}_2^0(\xi, \eta, \zeta) = \xi,$$ (2.42)

$${\ensuremath{\phi}}_3^0(\xi, \eta, \zeta) = \eta,$$ (2.43)

$${\ensuremath{\phi}}_4^0(\xi, \eta, \zeta) = \zeta.$$ (2.44)

Equation 2.37 is then calculated by:

$$M_{pq}(\ensuremath{\tau_{\mathrm{r},3}^{0}}) = \mathrm{det}(\ensu... ...nsuremath{\phi}}_q^0(\xi, \eta, \zeta)  d\xi d\eta d\zeta, \quad p,q = 1,2,3,4$$ (2.45)

where $\ensuremath{\mathbf{J}}$ is the Jacobian of the projection Equation 2.38 - Equation 2.40. The final assembly process consists of calculating stencil matrices $\Pi(\ensuremath{\tau_{n}^{i}})$ for each element of the cell complex, mapping the local vertex indices to a global index and then building the global matrix $\ensuremath{\mathbf{A}}$ .

For an electrostatic problem formulated by the Laplace equation $\mathrm{div} \left( \mathrm{grad} \left( {\ensuremath{\varphi}} \right) \right) = 0$ the stencil matrix is represented by:

$$\Pi(\ensuremath{\tau_{3}^{i}}) = M_{pq}(\ensuremath{\tau_{3}^{i}}... ...(\ensuremath{\mathbf{x}}) \right)  d{\ensuremath{\Omega}}, \quad p,q = 1,2,3,4$$ (2.46)

This problem is linear and the global matrix $\ensuremath{\mathbf{A}}$ is simply obtained by assembling $\Pi$ with the corresponding index transformation.