Finite element method

The Finite Element Method (FEM) is arguably the most powerful method known for the numerical solution of boundary- and initial-value problems characterized by partial differential equations. Consequently, it has had a monumental impact on virtually all areas of engineering and applied science.

There are two fundamental attributes of the method that are at the heart of its great utility and success. Firstly, it is based on the idea of partitioning bounded domains $\mathrm{\Omega }$ in ${\mathbb{R}}^n$ into a number $N$ of small, non-overlapping subdomains, the finite elements, over which functions are approximated by local functions, generally polynomials. Secondly, the boundary- and initial-value problems, to which the method is applied, are formulated in a so-called weak, or integral form, so that the contributions of each subdomain to the global integrals sum up to produce an integral characterizing the problem over the whole domain.

Some point to the first attribute of FEMs as its most important: the systemic representation of approximations of function spaces. From the viewpoint of approximation theory, the method provides a systematic approach to piecewise approximation over subdomains that produces sequences of functions that can approximate arbitrary members of, say, Sobolev spaces, arbitrarily closely in appropriate norms (e.g. Sobolev norms). This attribute is referred to as a property of "interpolation" for FEMs. It basically lifts the idea of approximate function representation from the classical notion of interpolation across values of functions or values of derivatives to a general approach for approximating functions with generalized derivatives in spaces such as $L^p\left(\mathrm{\Omega }\right)$ or $W^{m,p}\left(\mathrm{\Omega }\right),$ $(1\le p\le \mathrm{\infty },m\ge 0).$

But, were it only its ability to represent functions on quite general domains in ${\mathbb{R}}^n,$ the FEM would not have attained its popularity in science and engineering. It is the powerful notion of weak formulations of boundary- and initial-value problems that enables the summation of element contributions to integrals over the whole domain, a property that allows one to conceptualize a formulation of the problem at hand over a simple domain, the finite element, so as to approximate it locally. Then, these local subdomain approximations, as noted earlier, are summed up to determine the contributions from the full domain. This is the essence of the FEM.

The basic idea can be demonstrated by a simple model problem. Consider the Poisson problem on a two-dimensional domain which involves seeking a real-valued function $u=u(x,y)$ such that

$$\begin{array}{}\text{(1)}& \{\begin{array}{rl}-\mathrm{\Delta }u=f& \text{in}\mathrm{\Omega }\\ u=g& \text{on}\partial \mathrm{\Omega }\end{array}\end{array}$$

where $\mathrm{\Delta }=\text{div}\left(\text{grad}\right)={\mathrm{\nabla }}^2$ is the Laplacian, $f=f(x,y)$ is given source data, and $\mathrm{\Omega }$ is an open bounded domain in ${\mathbb{R}}^2$ with boundary $\partial \mathrm{\Omega }.$ Formally, the differential equation in (1) is multiplied by a test function $v=v(x,y),$ integrated by parts via Green's theorem, to obtain

$${\int }_{\mathrm{\Omega }}(-\mathrm{\Delta }u-f)vdxdy={\int }_{\mathrm{\Omega }}\mathrm{\nabla }u\cdot \mathrm{\nabla }vdxdy-{\int }_{\mathrm{\Omega }}fvdxdy-{\oint }_{\partial \mathrm{\Omega }}\mathbf{n}\cdot \mathrm{\nabla }uvds=0$$

where $\mathbf{n}$ is the outward unit vector normal to the boundary $\partial \mathrm{\Omega }.$ If now one demands that $v=0$ on $\partial \mathrm{\Omega },$ (1) is formally equivalent to the following weak or integral form of the boundary-value problem:

Find $u\in U$ such that

$$\begin{array}{}\text{(2)}& {\int }_{\mathrm{\Omega }}\mathrm{\nabla }u\cdot \mathrm{\nabla }vdxdy={\int }_{\mathrm{\Omega }}fvdxdy,\mathrm{\forall }v\in V.\end{array}$$

Here $U$ is called the space of trial functions and $V$ the space of test functions. Conceptually, the formulation (2) is "weaker" than (1) as it imposes weaker conditions on the smoothness of the solutions $u$ and test functions $v.$ In fact

$$U=\{w=w(x,y):w,\frac{\partial w}{\partial x},\frac{\partial w}{\partial y}\in L^2(\mathrm{\Omega }),w=0\text{on}\partial \mathrm{\Omega }\}$$ and, in this case, $V=U.$ The space $U$ is usually denoted

$$U=H_0^1(\mathrm{\Omega })=\{w\in H^1(\mathrm{\Omega }),w=0\text{on}\partial \mathrm{\Omega }\}$$

where $H^1(\mathrm{\Omega })$ is the Sobolev space of functions on $\mathrm{\Omega }$ with generalized derivatives $L^2(\mathrm{\Omega }),$ a special case of the $W^{m,p}(\mathrm{\Omega })$ of functions with derivatives of orders $m\ge 0$ in $L^p(\mathrm{\Omega }),$ $1\le p\le \mathrm{\infty }.$

Problem (2) is a simple example of a broad class of linear boundary- and initial-value problems of the form

$$\begin{array}{}\text{(3)}& \text{Find}u\in U\text{such that}B(u,v)=F(v),\mathrm{\forall }v\in V,\end{array}$$

where $B(\cdot ,\cdot )$ is a bilinear form mapping (reflexive Banach spaces) $U\times V$ into $\mathbb{R}$ and $F$ is a continuous linear functional defined on $V.$

The standard finite element approximation of (2) amounts to the construction of Galerkin or weighted residual approximations of (1) over finite dimensional subspaces $U^h$ and $V^h$ of $U$ and $V,$ respectively. That is, one seeks

$$\begin{array}{}\text{(4)}& u^h\in U^h\subset U\text{such that}B(u^h,v^h)=F(v^h),\mathrm{\forall }{v}^h\in V^h,\end{array}$$

with $V^h\subset V,$ $h$ being a discretization parameter generally chosen as a measure of the mesh size. The great utility of FEM's arises from its approach to constructing the approximation spaces $U^h$ and $V^h$ for functions defined on quite general domains $\mathrm{\Omega }.$

Returning to the example (1) (or (2)), consider the case in which $\mathrm{\Omega }$ is the domain shown in Figure 1, partitioned with a mesh of triangular elements, as shown. This mesh is viewed as one generated by a sequence of maps $F_k$ of a master triangle into the triangles so that, when connected along their edges and at their vertices, cover $\mathrm{\Omega }$ or produce a close approximation ${\mathrm{\Omega }}^h$ of $\mathrm{\Omega }.$ Suppose there are $N$ elements in the mesh. Over a typical isolated triangle ${\mathrm{\Omega }}_k,$ assume, for example, that the restriction $v_k^h$ of the test function $v^h\in V^h$ varies linearly. Then,

$$v_k^h(x,y)=v^h{|}_{{\mathrm{\Omega }}_k}(x,y)=a_1+a_{2}x+a_{3}y,(x,y)\in {\bar{\mathrm{\Omega }}}_k.$$

Figure 1: A triangular finite element partition of a domain $\mathrm{\Omega }$ into elements overwhich test functions are piecewise linear functions of the coordinates $(x,y).$

The constants $a_1,$ $a_2,$ and $a_3$ are determined so that $v^h$ is completely determined by its values at the vertices of ${\mathrm{\Omega }}_k,$ called the element nodes. If $v_1,$ $v_2,$ and $v_3$ denote the values of $V^h$ at nodes numbered 1, 2, and 3, then one easily determines that $v^h$ can be written,

$$v_k^h(x,y)=v_1{\psi }_1^k(x,y)+v_2{\psi }_2^k(x,y)+v_3{\psi }_3^k(x,y),k=1,2,\dots ,N,$$

where the ${\psi }_i^k$ are the element shape functions,

$$\{\begin{array}{l}{\displaystyle {\psi }_1^k(x,y)=\frac{1}{2A_k}[(x_2{y}_3-x_3{y}_2)+(y_2-y_3)x+(x_3-x_2)y]}\\ {\displaystyle {\psi }_2^k(x,y)=\frac{1}{2A_k}[(x_3{y}_1-x_1{y}_3)+(y_3-y_1)x+(x_1-x_3)y]}\\ {\displaystyle {\psi }_3^k(x,y)=\frac{1}{2A_k}[(x_1{y}_2-x_2{y}_1)+(y_1-y_2)x+(x_2-x_1)y]}\end{array}$$

$(x_i,y_i)$ are the coordinates of node $i,$ and $A_k$ is the area of ${\mathrm{\Omega }}_k.$ These linear functions have the property that

$${\psi }_i^k(x_j,y_j)=\{\begin{array}{l}1\text{if}i=j\\ 0\text{if}i\ne j\end{array}1\le i,j\le 3.$$

All the elements are then connected together to form the full computational domain $\mathrm{\Omega }$ (or an approximation ${\mathrm{\Omega }}^h$ of it), and in so doing create global piecewise linear basis functions ${\varphi }_i,$ corresponding to each node $i$ in the connected mesh. These "pyramid" or "hat" functions result from patching the element shape functions together, as indicated in Figure 2. Now the (global) test functions are linear combinations of these basis functions: For a mesh with $M$ interior nodes (boundary nodes are excluded since $u=0$ on the boundary),

$$v^h(x,y)=\sum _{i=1}^M{v}^h(x_i,y_i){\varphi }_i(x,y)$$ $$V^h=\text{span}{\{{\varphi }_i(x,y)\}}_{i=1}^M$$

Figure 2: A global basis function ${\varphi }_i$ for the spaces $U^h$ and $V^h$ made up by piecing together element shape functions. Note that the basis function is defined everywhere in $\mathrm{\Omega }$ but is nonzero only over the patch of elements surrounding node $i.$

In the finite element approximation of (2), one seeks

$$u^h\in U^h=V^h:u^h(x,y)=\sum _{i=1}^M{u}_i{\varphi }_i(x,y)$$

that satisfies the weak problem (4) where

$$\{\begin{array}{ll}& {\displaystyle B(u^h,v^h)={\int }_{\mathrm{\Omega }}\mathrm{\nabla }{u}^h\cdot \mathrm{\nabla }{v}^{h}dxdy=\sum _{i,j=1}^M{v}_i{K}_{ij}{u}_j}\\ & {\displaystyle F(v^h)=\sum _{i=1}^M{f}_i{v}_i}\end{array}$$

for arbitrary coefficients $v_i.$ Here $K_{ij}$ are the elements of the global stiffness matrix:

$$K_{ij}={\int }_{\mathrm{\Omega }}\mathrm{\nabla }{\varphi }_j\cdot \mathrm{\nabla }{\varphi }_{i}dxdy$$

and $f_i$ are the components of the load vector,

$$f_i={\int }_{\mathrm{\Omega }}f{\varphi }_{i}dxdy.$$

This leads to a system of linear equations for the coefficients $u_j,$ called the degrees of freedom:

$$\begin{array}{}\text{(5)}& \sum _{j=1}^M{K}_{ij}{u}_j=f_i,1\le i\le M.\end{array}$$

The final solution is then,

$$u^h(x,y)=\sum _{i=1}^M\left(\sum _{j=1}^M{\left\{K^{-1}\right\}}_{ij}{f}_j\right){\varphi }_i(x,y)$$

A property of fundamental importance is that this entire process can essentially be done at the local element level, and then summed up to apply to the global version of the problem. Thus, on element ${\mathrm{\Omega }}_k,$ the local equalities for element $k$ are

$$\sum _{j=1}^3{K}_{ij}^k{u}_j^k=f_i^k+g_i^k(u),1\le i\le 3.$$

where $K_{ij}^k$ and $f_i^k$ are the elements of the local stiffness matrix and of the local load vector for element $k,$ respectively, and $g_i^k(u)$ are element interface flux terms:

$$\begin{array}{}\text{(6)}& \{\begin{array}{l}{\displaystyle K_{ij}^k={\int }_{{\mathrm{\Omega }}_k}\mathrm{\nabla }{\psi }_j^k\cdot \mathrm{\nabla }{\psi }_i^{k}dxdy}\\ {\displaystyle f_i^k={\int }_{{\mathrm{\Omega }}_k}f{\psi }_i^{k}dxdy}\\ {\displaystyle g_i^k(u)={\oint }_{\partial {\mathrm{\Omega }}_k}\mathbf{n}\cdot \mathrm{\nabla }u{\psi }_i^{k}ds}\end{array}\end{array}$$

In the expression for $g_i^k(u),$ $u$ is the exact solution of the boundary-value problem (2); while it is unknown, it does not enter the discrete global problem because

$$\sum _{k=1}^M{g}_i^k(u)=0.$$

In practice, the local element matrices (6) are integrated using numerical quadrature methods when higher-order shape functions are used. In general, FEM computer programs are written so that element matrices such as those in (6) are calculated from a family of element types, and these are used to generate the global matrices of (5) when the elements are connected together to form the global element mesh. Of course, there are many other possible choices for finite elements, including for instance, quadratic approximations on triangles with 6 nodes, cubic with 10 nodes, quadrilateral elements, etc., and generalizations to tetrahedral and hexahedral elements for three-dimensional elements are commonly used.

The mathematical properties of FEM's and the fact that convergent and accurate approximations of solutions to problems defined on domains of quite general geometry are worth noting. Returning to the general weak boundary-value problem (3), one observes that $B(\cdot ,\cdot )$ is typically assumed to be continuous on $U\times V$ into $\mathbb{R}$ with continuity constant $M,$ and to satisfy an "inf-sup" condition

$$\begin{array}{}\text{(7)}& \underset{u}{inf}\underset{v}{sup}\frac{\left|B(u,v)\right|}{{\Vert u\Vert }_U{\Vert v\Vert }_V}\ge \gamma >0,\end{array}$$

where $\gamma$ is a constant. Under the conditions that $B(\cdot ,\cdot )$ and $F(\cdot )$ are continuous (on $U\times V$ and $V,$ respectively), that $\underset{u}{sup}\left|B(u,v)\right|>0$ for all $v\in V,$ and that (7) holds, it can be shown that (3) possesses a unique solution.

To determine the FE approximation, the domain $\mathrm{\Omega }$ is partitioned into $N$ subdomains ${\mathrm{\Omega }}_k$ such that

$$\bar{\mathrm{\Omega }}=\text{closure of}\mathrm{\Omega }={\cup }_k{\bar{\mathrm{\Omega }}}_k,{\mathrm{\Omega }}_k\cap {\mathrm{\Omega }}_j=\text{empty},k\ne j.$$

As noted earlier, each element ${\mathrm{\Omega }}_k$ is regarded as the image of a master subdomain (a master element) $\hat{\mathrm{\Omega }}$ under an invertible element map $F_k:\hat{\mathrm{\Omega }}\to {\mathrm{\Omega }}_k,$ $1\le k\le N.$ In this way, a sequence of subspaces $\left\{V^h\right\}$ is constructed such that

$$V^h=\{v^h\in C^0(\bar{\mathrm{\Omega }}):{v^h|}_{{\mathrm{\Omega }}_k}\circ F_k^{-1}=\hat{v}\in {\mathbb{P}}_k\},$$

where ${\mathbb{P}}_k$ is the space of local FEM approximations, e.g.

$$\begin{array}{}\text{(8)}& {\mathbb{P}}_k={\mathcal{P}}_p(\hat{\mathrm{\Omega }})=\text{space of polynomials of degree}\le p\text{on}\hat{\mathrm{\Omega }},\end{array}$$

or $${\mathbb{P}}_k={\mathcal{P}}_p(\hat{I})\otimes {\mathcal{P}}_p(\hat{I})=\text{space of tensor products of polynomials on}\hat{\mathrm{\Omega }}=\hat{I}\otimes \hat{I},$$

Then the FEM approximation of (3) reduces to the sequence of discrete problems (4). From (3) and (4), one can derive the standard "orthogonality condition"

$$B(u-u^h,v^h)=0,\mathrm{\forall }{v}^h\in V^h.$$

So, in this sense, the FEM (Galerkin) approximation of (4) is optimal: the error $e^h=u-u^h$ is "orthogonal" to the space $V^h$ of approximations with respect to the bilinear form $B(u,v).$

Under reasonable assumptions on the sequence of meshes producing the subspaces $\{V^h\},$ one can approximate functions in $W^{r,s}\left(\mathrm{\Omega }\right)$ with interpolation error satisfying

$$\Vert u-{\mathrm{\Pi }}_{h}u{\Vert }_{r,s,{\mathrm{\Omega }}_k}\le C{h}_k^{min(p+1-r,t-r)}\Vert u{\Vert }_{W^{t,s}(\mathrm{\Omega })}$$

where $p$ is the highest order of complete polynomials contained in the space of local FEM approximations (as in (8)), ${\mathrm{\Pi }}_h$ is a projection operator of $V^h$ into $W^{r,s}(\mathrm{\Omega })$ that preserves polynomials of degree $\le p$ for triangular elements, $h_k=\text{diam}({\mathrm{\Omega }}_k),$ $1\le s\le \mathrm{\infty },$ $t\ge 0,$ and $C$ is a constant.

When interpolation properties such as this hold, one can typically obtain a priori error estimates for elliptic boundary-value problems of the form $$\Vert u-u_h{\Vert }_{H^m(\mathrm{\Omega })}\le C{h}^{min(p+1-m,q-m)}\Vert u{\Vert }_{H^q(\mathrm{\Omega })}$$

where now $s=2,$ $r=m,$ $t=q,$ and $h=\underset{1\le h\le N}{max}\left(h_k\right).$ Thus, the rate of convergence of such methods in the $H^m(\mathrm{\Omega })$-norm is $min(p+1-m,q-m)>0.$ The constant $C$ is independent of $u$ and $h$ but depends upon the continuity constant $M$ of $B(\cdot ,\cdot )$ and the constant $\gamma$ of (7).

Many variations of such results can be obtained. Among them are cases in which the local spectral order $p=p_k$ is an approximation parameter. These characterize the "$p$-versions" of the FEM. Likewise, both $h_k$ and $p_k$ can be varied, resulting in $h$-$p$ FEMs. The "$h$-versions" of the FEM commonly refers to the case where $p$ is uniform over all elements.

The field of FEM's, both as a tool for solving problems of science and engineering and as an area of research in numerical mathematics, is large and expanding and continues to be a vital area of study and practical application.

[edit] References