Introduction to the Finite Element Method

1 Introduction to the Finite Element MethodMarek K. MisztalNiels Bohr 28, 20161 IntroductionThe Finite Element methods (FEM) is a family of numerical methods designedto find approximate solutions to linear, partial differential equations (PDEs)on a variety of domains. It relies on domain decomposition into simplerele-ments(triangles, tetrahedra, quadrilaterals, etc.) and interpolation of thesolution based on a discrete set of values defined at specified domain applying the Galerkin Method , the original PDE is converted into aweakformulation, which can be written as a linear system of equations.

2 TheGalerkin Method ensures a solution that minimizes (although not strictly)the residual between the actual solution and the functions from the space ofapproximate note is meant as a practical guide to obtaining the linear, Finite elementformulations of Poisson-like PDEs, with emphasis on deriving equations suit-able for numeric implementation. Because of that, this note makes use ofthe matrix notation, rather than tensor Focus problemsProblem A (Poisson s equation)Let be a bounded, compact subsetofR2. Find suchu: R, that 2u= 0(1)u(x) =u0(x)(2)Problem B (Stokes equation)*Compute the laminar, steady state flowthrough a Finite pipe (in 2D), due to applied pressure difference between the1inlet and the R2is the geometric representation of the pipe.

3 We introduce into denote the inlet surface, outto denote the outlet surface, in out= . Further, s= ( in out) denotes the solid wall. Formally, weseek suchu: R2andp: R, that 2u p+f= 0,(Stokes equation)(3) u= 0,(continuity equation)(4)subject to the Dirichlet boundary conditionsu(x) s= 0,(5)p(x) in=pin,(6)p(x) out=pout,(7)wherepinandpoutare constant Domain discretizationIn order to numerically solve the PDE, we first have to approximate thedomain with a Finite number of simpler, geometric objects (referred to aselements).

4 Popular types of elements include simplicial: segments in 1D, triangles in 2D, tetrahedra in 3D, quadrilateral (in 2D), hexahedral (in 3D).From now we focus on simplicial domains in 2D, which are known astrianglemeshes. Most of the following equations significantly simplify in the 1D case,and the overall procedure is very similar in the 3D case. We will denote thevertices(ornodes) of the meshxi,i= 1,2,..,N, and the triangle elements e,e= 1,2,..,M. Typically, the number of elements in the meshMscales linearly with the number of verticesN. elements are only allowed tointersect along the common edge (no T-junctions allowed).

5 We define theentire discrete domain as =M e=1 e.(8)2 Once we have decided on the domain type, we have to decide on how to dis-cretize the solution. In this note, we focus on linear schemes, ones wherethe approximate solution is given as a piece-wise linear function. Higher-order schemes are often more robust, however, they require maintaining thevalues of the solution s derivatives. Having a piecewise linear approximationof the solution u(x) requires us to store only a Finite number of its values,typically at the vertices or elements centres.

6 For now we will focus on vertex-centred schemes, where we specifyui= u(xi), fori= 1,2,.., Solution spaceThe decisions we have made about the discretization determines thesolutionspace, the set of admissible solutions. These are constructed throughinterpolation, based on the nodal values, usinginterpolantfunctions ifori= 1,2,..,N, such that i(xj) ={1 ifi=j0 ifi6=j,N i=1 i(x) = 1.(9)Then the interpolated (approximate) solution is defined as u(x) =N i=1ui i(x) = (x)T u,(10)whereui=u(xi) fori= 1,2,..,Nare the (stored) nodal values of the solu-tion, (x) = [ 1(x), 2(x).]}

7 , N(x)]Tand u= [u1,u2,..,uN]T. Formally,the solution space is defined as:S={f: T; ( TN)f(x) = (x)T }.(11)Here we focus on linear elements with linear shape functions. In this case i(x)| e= const. The gradient of the interpolated function uis well-definedin the interior of each Element e, and is given by u(x) =N i=1ui i(x) = ( (x))T u,(12)whereRN 23 (x) =[ x 1(x) x 2(x).. x N(x) y 1(x) y 2(x).. y N(x)]T.(13) Galerkin methodThe Galerkin Method allows us to find the quasi-best approximation of thereal solution to the PDE among the functions from the solution specific details of the mathematics behind the Galerkin Method are outof scope of this note, and are available in specialized literature.

8 In short, theGalerkin solutionug Sis found by transforming the PDELu= 0,(14)whereLis a linear operator (for exampleL= 2for the Poisson problem)into itsweak formulation ( w T) w,L u = 0,(15)where is the inner product in the space of the solution functions, andTisthe set of admissibletest functions( w),T={v S:v(x)| = 0}. UndertheL2-norm, this becomes ( w T) w(x)TL u(x) d = 0.(16)3 Poisson s Galerkin Method in practiceThe weak formulation of the Problem A reads ( w T) wT 2 ud = 0,(17)where u, w:R2 R. By applying the divergence theorem, we obtain wT 2 ud = wT( u)T ndS ( w)T ud ,(18)where nis the normal vector to.

9 Note that becausew(x)| = 0, the firstterm on the right hand side vanishes. Now recall that u(x) = ( (x))T uand w(x) = ( (x))T w, which allows us to rewrite Eq. (18) into a linearequation w(x)T u(x) d = wT (x) ( (x))T ud (19)= wT{ (x) ( (x))Td } u(20)= wTK u,(21)4whereKis anN Nreal matrix (called thestiffness matrixof the PDE),with coefficients defined asKij= i(x) j(x)d . We have arrived atan equivalent form of Eq. (17) ( w RN) wTK u= 0,(22)which is fulfilled by usuch thatK u=0.(23) Matrix assembly*We can learn more about the structure of matrixKby rewriting the integralover as a sum of integrals over each Element eK= (x) ( (x))Td =M e=1 e e( e)Td =M e=1Ke,(24)where e= (x)| e= const, andKeis theelement stiffness the 2D case with triangle elements , only three shape functions are non-zero over a given Element e(the three shape functions associated with thetriangle s vertices).

10 That means, only three out ofNrows of eare non-zero. That also means thatKehas at most nine non-zero coefficients, forminga 3 3 submatrixKe3=| e| ei ei ei ej ei ek ej ei ej ej ej ek ek ei ek ej ek ek ,(25)where| e|is the volume (area) of the Element e. Note that matrixKhas atmost 9 Mnon-zero coefficients, which is one order of magnitude lower thanits sizeN2. Thissparsityproperty is important from the point of view ofnumeric performance, allowing to solve the linear system using fast, iterativesolvers (such as the conjugate gradient, GMRES, etc.)