Finite Element Method: One Dimension - GitHub Pages

University of Notre Dame Dept Aerospace & Mechanical Engrng M. J. Zahr 1) the functions are nodal, i.e., each node I “ 1, . . . , Nv has an associated basis function the Lagrangian property I px̂J q “ IJ I that possesses (4.19) for J “ 1, . . . , Nv , 2) the functions are linear when restricted to an element of Eh , i.e., for any e P Eh I e P P 1 p e q (4.20) for I “ 1, . . . , Nv , and 3) the basis function associated with node I is only non-zero on elements connected to node I, i.e., the basis functions have local support supp 1 “ 1 , supp Nv “ Nv 1 , supp I “ I 1 Y I , I “ 2, . . . , Nv 1. (4.21) The first property implies A is linearly independent (Proposition 4.1). It also provides a strong connection between values of a function in span A and its expansion in the basis A. To see this take f P span A, which we write as Nv ÿ f“ fˆI I , (4.22) I“1 where fˆ1 , . . . , fˆNv P R are the coefficients defining the expansion of f in the basis A. This is an equality between functions and must hold for any x P . Evaluating at a node x “ x̂J , we have f px̂J q “ Nv ÿ fˆI I“1 I px̂J q “ Nv ÿ fˆI IJ I“1 “ fˆJ , (4.23) which shows that the value of the function at node J is equal to the coefficient associated with the basis function J , i.e., the coefficients are the nodal values of the function (intuition behind the nodal terminology). The second property implies that Vh0 Ä span A. Finally, the last property leads to sparsity of the finite element sti ness matrix, which is of extreme practical importance as we will see in later chapters. Proposition 4.1 (Linear independence of nodal functions). Consider the domain Ä Rd with corresponding triangulation Eh and node set Nh (Nv “ Nh ). Let A “ t 1 , . . . , Nv u Ä F ÑR be a collection of functions with the nodal property I px̂J q “ IJ for I “ 1, . . . , Nv and x̂J P Nh . Then the set A is linearly independent. Proof. Take 1 , . . . , Nv P R such that Nv ÿ I I I“1 “ 0. (4.24) Since this is an equation involving functions, it must hold for any x P . Evaluating this equation at a node x “ x̂J , we have Nv Nv ÿ ÿ I I px̂J q “ I IJ “ J “ 0, (4.25) I“1 I“1 which holds for all nodes J “ 1, . . . , Nv . Therefore the functions in A are linearly independent. The collection of functions over the one-dimensional nodal hat functions, defined as x x̂I 1 ’ ’ ’ ’ ’ & hI 1 x̂I 1 x :“ pxq I ’ ’ hI ’ ’ ’ % 0 domain p0, Lq satisfying these properties are the x P I 1 x P I (4.26) x R I 1 Y I Page 53 of 75

University of Notre Dame Dept Aerospace & Mechanical Engrng M. J. Zahr 1 Nv I 1 x̂1 x̂2 x̂I 1 x̂I x̂I 1 x̂Nv 1 x̂Nv Figure 4.4: Hat function for interior node I ( ) and for boundary nodes 1 ( ) and Nv ( ). The functions are only non-zero in elements connected to the particular node and take the value of 1 at their own node and 0 at all other nodes. 1 2 3 4 5 6 x̂1 x̂2 x̂3 x̂4 x̂5 x̂6 1 Figure 4.5: Hat functions for a complete triangulation (Ne “ 5) of a one-dimensional domain Ä R. for interior nodes I “ 2, . . . , Nv 1 and & x̂2 x x P 1 h1 :“ pxq 1 % 0 x R 1 , & x x̂Nv 1 hNv 1 :“ pxq Nv % 0 x P Nv 1 (4.27) x R Nv 1 , for boundary nodes (Figure 4.4-4.5). It is a simple exercise to verify the hat functions satisfy our requirements (piecewise linear nodal functions with local support). In addition, the nodal hat functions are continuous, i.e., I P C 0 pp0, Lqq, which implies any element of span A are continuous (linear combinations of continuous functions are continuous functions). We use these functions to define the FE affine o set 'h and linear subspace Vh0 as 'h :“ ū 1, Vh0 :“ spant 2, . . . , Nv u (4.28) (Figure 4.3). For these to be valid choices according to (4.18), the affine o set must be continuous and satisfy 'h p0q “ ū and any function vh P Vh0 must be continuous and satisfy vh p0q “ 0. Futhermore the Nv 1 vectors defining Vh0 must be linearly independent to span the pNv 1q-dimensional space defined in (4.18). Continuity of 'h and elements of Vh0 follow directly from continuity of the nodal hat functions 1 , . . . , Nv as mentioned previously. Boundary condition enforcement follows directly from the Lagrangian property 'h p0q “ ū 1 p0q “ ū 1 px̂1 q “ ū (4.29) and vh p0q “ vh px̂1 q “ v̂2 2 px̂1 q v̂Nv Nv px̂1 q “ 0, (4.30) where vh P is expanded in the basis t 2 , . . . , Nv u with corresponding coefficients v̂2 , . . . , v̂Nv P R. The Nv 1 vectors defining Vh0 are linearly independent owning to Proposition 4.1. Therefore, the choices in (4.28) are valid. Vh0 4.4.2 Construction from element level Direct construction of the finite element subspace becomes cumbersome in higher dimensions and for higher polynomial degrees. Therefore we consider an alternate approach that, instead of considering the entire domain “ p0, Lq and building up basis functions for Vh0 in a global sense (Section 4.4.1), constructs the finite element space restricted to an arbitrary element and then fits elements together using the continuity requirement. This simplifies the task of constructing a basis for a piecewise polynomial space over the triangulation Eh to constructing one for a polynomial space over e P Eh . Unlike the direct approach, this Page 54 of 75

University of Notre Dame Dept Aerospace & Mechanical Engrng M. J. Zahr e 1 pxq e 2 pxq xe1 xe2 e Figure 4.6: A finite element with local node numbering and P 1 p q nodal basis. construction will lead to a systematic procedure, amenable to computer implementation, that can naturally be generalized to higher dimensions, higher polynomial spaces, and other PDEs. Because Vh0 consists of piecewise linear functions, for any vh P Vh0 we have vh e P P 1 p e q (4.31) for any e P Eh , where vh e : e Ñ R is the restriction of vh to e . In words this says any function of Vh0 is a polynomial of degree one when restricted to an element of Eh . This implies that the restriction of vh to e can be written vh e “ v̂1e e1 v̂2e e2 , (4.32) where t e1 , e2 u is a basis of P 1 p e q (recall dim P 1 p e q “ 2) and v̂1e , v̂2e P R are the corresponding coefficients. For convenience we choose a nodal basis that satisfies e e i px̂j q “ (4.33) ij for i, j “ 1, 2. Repeating the development from the previous section, the nodal property guarantees the functions t e1 , e2 u are linearly independent and, for any function f P P 1 p e q with coefficients fˆ1e , fˆ2e P R in the basis t e1 , e2 u, i.e., f “ fˆ1e e1 fˆ2e e2 , f px̂ei q “ 2 ÿ fˆje j“1 “ fˆie . e e j px̂i q (4.34) In this sense, ei is the basis function associated with node x̂ei and the corresponding coefficient fˆie is equal to the value of the function f at that node. The unique nodal basis of P 1 p e q (Figure 4.6) is e 1 pxq Since basis I :“ x̂e2 x , he e 2 pxq :“ x x̂e1 . he (4.35) P Vh0 for I “ 1, . . . , Nv , its restriction to an element e P Eh can be expanded in the element I e “ 2 ÿ ˆ ie e i. (4.36) i“1 Furthermore since I is a nodal function, it will take the value 1 at node I and zero at all other nodes, which leads to the condition e (4.37) I px̂i q “ I px̂ ie q “ I ie where the first equality used the global-to-local node relationship in (4.14) and the second used the nodal property of I . Substituting this into equation (4.36), we have I ie “ e I e px̂i q “ 2 ÿ j“1 ˆ je e e j px̂i q “ ˆ ie , (4.38) Page 55 of 75

University of Notre Dame Dept Aerospace & Mechanical Engrng M. J. Zahr e 1 e e 1 e 2 x̂e “ x̂e1 e x̂e 1 “ x̂e2 Figure 4.7: The relationship between the global finite element basis functions I ( ) and the element basis functions ei ( ): the element basis functions are the restriction of the global basis functions to a single element. Shaded region: restriction to element e . which leads to the following relationship between the element basis t t 1 , . . . , Nv u 2 ÿ e “ I e I ie i . e 1, e 2u and the global functions (4.39) i“1 This procedure to construct the global basis functions from elementwise basis functions will naturally generalize to higher dimensions and polynomial degrees. In the special case of a scalar PDE on a onedimensional mesh (numbering given in Figure 4.2), we use the fact that 1e “ e and 2e “ e 1 to reduce this to e e (4.40) I e “ Ie 1 Ie1 2 , where e1 “ e 1. For concreteness consider an interior node I “ 2, . . . , Nv 1, then x x̂I 1 ’ I ’ x P I 1 ’ 2 “ ’ hI 1 ’ & x̂I 1 x I I pxq “ x P I ’ 1 “ ’ hI ’ ’ ’ % 0 x R I 1 Y I , (4.41) which agrees with the expression from direct construction of I in (4.26). Figure 4.7 illustrates this relationship between the local and global basis functions. This completes our construction of the Vh0 from the element level. In an implementation, we never need to explicitly form the piecewise polynoimal basis functions t 1 , . . . , Nv u; however, this relationship between the element and global bases will be helpful in deriving the finite element system in a way that will suggest an e ective, modular computer implementation. 4.5. Finite element method Given the construction of the finite element subspace in the previous section, we turn to deriving the finite element system. That is, introducing our subspace approximation into the finite element weak formulation to obtain an algebraic system of equations. Following the precedent set in the previous section, we first consider direct construction of the global sti ness matrix and load vector. While instructive this is not amenable to computer implementation in more complex scenarios. We repeat the derivation using the elementwise construction in Section 4.4.2 to obtain the global algebraic system from elementwise terms. 4.5.1 Direct construction of global system From the construction of the nodal basis of hat functions (4.28) of the FE subspace Vh0 in Section 4.4.2, any element wh P Vh0 and uh P 'h Vh0 can be written as wh “ Nv ÿ I“2 ŵI I, uh “ ū 1 Nv ÿ ûI I, (4.42) I“2 Page 56 of 75

University of Notre Dame Dept Aerospace & Mechanical Engrng M. J. Zahr where ŵ2 , . . . , ŵNv P R and û2 , . . . , ûNv P R are the coefficients of wh and uh , respectively, in the basis t 2 , . . . , Nv u. We introduce ŵ1 , û1 P R so wh and uh can be conveniently written as wh “ Nv ÿ ŵI uh “ I, I“1 Nv ÿ ûI I, (4.43) I“1 provided w1 “ 0 and u1 “ ū; these boundary conditions will be imposed later. Substituting the expansions (4.43) into the FE formulation (4.8) and using linearity of the functionals yields « ff Nv Nv ÿ ÿ ˆ ŵI K̂IJ ûJ fI “ 0 (4.44) I“1 J“1 where K̂ P MNv ,Nv pRq is the global sti ness matrix and fˆ P RNv is the global load vector without considering essential boundary conditions K̂IJ :“ Bp I , J q, fˆI :“ p I q. (4.45) Now we impose the essential boundary conditions, i.e., require ŵ1 “ 0 and û1 “ ū, to yield « ff Nv Nv ÿ ÿ ˆ ŵI K̂IJ ûJ pfI ūK̂I1 q “ 0. I“2 (4.46) J“2 This equation must hold for arbitrary ŵ2 , . . . , ŵNv P R to be equivalent to the FE formulation (4.8), which holds for an arbitrary wh P Vh0 . Equation (4.46) can only be true for arbitrary values of ŵ2 , . . . , ŵNv if each term in the summation over I is zero, i.e., Nv ÿ J“2 K̂IJ ûJ “ fˆI ūK̂I1 (4.47) for I “ 2, . . . , Nv , which is a (square) linear system of equations (of size Nv 1) that can be solved to compute the unknown coefficients û2 , . . . , ûNv . Then the FE solution is reconstructed as uh “ ū 1 û2 2 ûNv Nv . (4.48) To provide a strong connection to the direct sti ness method, we introduce terms that partition quantities based on whether or not an essential boundary condition is prescribed at the corresponding node. Let I c denote the set of node indices constrained by an essential BC and let I u be the unconstrained nodes (without essential BC). Then we use a superscript v̂ u to denote the restriction of a vector over all nodes (v̂ P RNv ) to only the nodes in I u nodes; v̂ c is defined similarly. For matrices (MNv ,Nv pRq) two superscripts are required to specify the restriction of the rows and columns, e.g., for  P MNv ,Nv pRq, Âcu is the restriction of rows of  to those in I c and the restriction to columns in I u . In the present context, I c “ t1u and I u “ t2, . . . , Nv u. With this notation set, define K̂ uu P MNv 1,Nv 1 pRq, K̂ uc P MNv 1,1 pRq, ûu P RNv 1 , ûc P R, and fˆu P RNv 1 where uu uc :“ K̂I 1,J 1 , :“ K̂I 1,1 , K̂IJ K̂I1 ûuI :“ ûI 1 , ûc :“ ū, fˆIu :“ fˆI 1 (4.49) for I, J “ 1, . . . , Nv 1. With these definitions (4.47) reduces to K̂ uu ûu “ fˆu K̂ uc ûc , (4.50) which provides a close parallel to the final system obtained using the direct sti ness method (trusses) in (1.38). Remark 4.2. As mentioned in Section 4.3, the finite element method is a special case of the Ritz method developed in Section 3.6 using the specific subspace Vh0 and affine o set 'h constructed in Section 4.4. To see this, introduce I :“ I 1 for I “ 1, . . . , Nv 1 and recall that t 1 , . . . , Nv 1 u was established to be a basis of Vh0 in Section 4.4.1. From (3.57), applying the Ritz method to the bilinear form in (4.4) with this basis leads to the sti ness matrix ritz K̂IJ “ Bp I, Jq “ Bp I 1 , J 1 q uu “ K̂I 1,J 1 “ K̂IJ (4.51) Page 57 of 75

University of Notre Dame Dept Aerospace & Mechanical Engrng M. J. Zahr for I, J “ 1, . . . , Nv 1 and load vector fˆIritz “ p Iq Bp I , 'h q “ p I 1 q ūBp 1q I 1 , uc c “ fˆI 1 ūK̂I 1,1 “ fˆIu K̂I1 û1 (4.52) for I “ 1, . . . , Nv 1, which is precisely the FE system in (4.50). We close this section by explicitly constructing the global sti ness matrix and force vector corresponding to the PDE in (4.1) based on the basis of nodal hat functions. From the definition of the sti ness matrix in (4.45) and the bilinear terms in (4.5), the entries of the sti ness matrix are K̂IJ :“ Bp I, Jq “ ªL 0 d I d J a dx dx dx (4.53) for I, J “ 1, . . . , Nv . From the expression for the hat functions we observe the integrand will only be non-zero if I J § 1, otherwise the support of the basis functions in the integrand will not overlap and the product will be zero. This reduces the sti ness matrix to K̂IJ ª x̂2 d 1 d 1 ’ ’ a dx ’ ’ ’ dx x̂1 dx ’ ’ ’ ª ’ ’ x̂I 1 d I d I ’ ’ a dx ’ ’ ’ dx x̂I 1 dx ’ ’ ’ª x̂ ’ Nv ’ d Nv d Nv ’ & a dx dx “ x̂Nv 1 dx ’ ª x̂I ’ ’ ’ d I d I 1 ’ ’ a dx ’ ’ dx dx ’ x̂ I 1 ’ ’ ’ ª x̂I 1 ’ ’ d I d I 1 ’ ’ a dx ’ ’ ’ dx dx x̂I ’ ’ % 0 if I “ J “ 1 if 1 † I “ J † Nv if I “ J “ Nv (4.54) if J “ I 1, I 1 if J “ I 1, I † Nv otherwise for I, J “ 1, . . . , Nv . Since the nodal hat functions are piecewise linear, their derivatives are piecewise constant 1 ’ ’ ’ ’ & hI 1 d I pxq :“ 1 ’ dx ’ hI ’ ’ % 0 x P I 1 (4.55) x P I x R I 1 Y I for interior nodes I “ 2, . . . , Nv 1 and & 1 d 1 h1 :“ pxq % dx 0 x P 1 x R 1 , d Nv dx pxq :“ & 1 hNv 1 % 0 x P Nv 1 (4.56) x R Nv 1 , Page 58 of 75

University of Notre Dame Dept Aerospace & Mechanical Engrng M. J. Zahr boundary nodes. From these equations, the sti ness matrix reduces to K̂IJ ª x̂2 1 ’ ’ a dx ’ ’ ’ h21 x̂1 ’ ’ ’ ª x̂I ª ’ ’ 1 1 x̂I 1 ’ ’ a dx a dx ’ ’ ’ h2I 1 x̂I 1 h2I x̂I ’ ’ ’ ª x̂Nv ’ ’ ’ & 1 a dx “ h2Nv 1 x̂Nv 1 ’ ª x̂I ’ ’ ’ ’ 1 ’ a dx ’ ’ ’ h2I 1 x̂I 1 ’ ’ ’ ª ’ ’ 1 x̂I 1 ’ ’ a dx ’ 2 ’ ’ h ’ ’ I x̂I % 0 if I “ J “ 1 if 1 † I “ J † Nv if I “ J “ Nv (4.57) if J “ I 1, I 1 if J “ I 1, I † Nv otherwise. This expression highlights one critical reason for choosing basis functions with local support: the sti ness matrix is sparse, i.e., many of its entries are zero. Furthermore, in the special case of a scalar PDE in one dimension (with the numbering chosen according to Figure 4.2), the matrix is tridiagonal, i.e., of the form » ˆ —ˆ — — K̂ “ — — – ˆ ˆ . . ˆ . . ˆ fi . . ˆ ˆ ffi ffi ffi ffi ffi ˆfl ˆ (4.58) where the ˆ symbol indicates a nonzero value. This special structure means that, in addition to being very sparse, specialized algorithms exist for solving linear systems defined by K̂ extremely efficiently using direct methods (OpNv q operations, whereas direct solvers usually require OpNv3 q operations). To complete the Ritz method, we compute the load vector fˆ P RNv as fˆI :“ p Iq “ I pLqQ̄ ªL If dx (4.59) 0 for I “ 1, . . . , Nv , which can be reduced to ª x̂2 ˆ x̂2 x ’ ’ f dx ’ ’ h1 ’ x̂1 ’ ’ ’ ª x̂I 1 ˆ &ª x̂I ˆ x x̂ x̂I 1 x I 1 ˆ f dx f dx fI “ hI 1 hI ’ x̂I 1 x̂I ’ ’ ’ ª x̂Nv ˆ ’ ’ x x̂Nv 1 ’ ’ f dx %Q̄ hNv 1 x̂Nv 1 if I “ 1 if 1 † I † Nv (4.60) if I “ Nv using the expressions for the hat functions (4.26)-(4.27) and their derivatives (4.55)-(4.56). The approach to directly construct the sti ness matrix from its definition in (4.53) is striaghtforward; however, it is not practical for a number of reasons. The main reason is it relies on explicit expressions for the basis functions 1 , . . . , Nv , which are messy to define for higher degree polynomial spaces and higher dimensional problems. In addition, it does not appear amenable to computer implementation for more complex problems due to the presence of the integrals, each of which must consider two elements (in general) to carry out the computation. In higher dimensions on unstructured meshes, the number of elements connected to a given node varies throughout the mesh, which makes this procedure more complicated. Page 59 of 75

University of Notre Dame Dept Aerospace & Mechanical Engrng 4.5.2 M. J. Zahr Construction of global system from element level assembly Now we turn to an alternative procedure to form the finite element system from element basis functions that that can be easily extended to higher dimensional space and higher degrees polynomials. To begin, we observe that the bilinear form in (4.4) can be written equivalently, owing to the additive property of integration, as a summation over element contributions Bpw, uq “ Ne ÿ e“1 Be pw, uq, pwq “ Ne ÿ e“1 e pwq, (4.61) where Be : W ˆ V Ñ R is the restriction of the bilinear functional B to element e P Eh and e : W Ñ R is the restriction of to e . In the case of the PDE in (4.5), these terms are # ª ª wpLqQ̄ if e “ Ne dw du Be pw, uq :“ a dx, e pwq :“ wf dx , (4.62) 0 otherwise. e dx dx e The boundary term is only included in the element touching the boundary x “ L to avoid counting it multiple times. From the definition of the elementwise functionals as integrals over a single element e P Eh , we have Be pw, uq “ Be p w e , u e q, e pwq “ e p w e q (4.63) for any w P W and u P V. With these definitions, the sti ness matrix K̂ in (4.45) can be expanded as K̂IJ :“ Bp I, Jq “ Ne ÿ e“1 Be p I, Jq “ Ne ÿ e“1 Be p I e , J e q “ Ne ÿ 2 ÿ 2 ÿ e“1 i“1 j“1 e i, Be p e j q I ie J je , (4.64) for I, J “ 1, . . . , Nv and we used (4.61) to write the bilinear term as a summation over element contributions, (4.63) to restrict the basis functions to e , and (4.39) to relate the global basis functions to the element basis functions. Next we expand the load vector fˆ in (4.45) as fˆI :“ p Iq “ Ne ÿ e“1 e p Iq “ Ne ÿ e“1 e p I e q “ Ne ÿ e“1 e I ie e p i q (4.65) for I “ 1, . . . , Nv and we used (4.61) to write the linear functional as a summation over element contributions, (4.63) to restrict the basis functions to e , and (4.39) to relate th

4.3. Finite element method: formulation The finite element method is a Ritz method in that it approximates the weak formulation of the PDE in a finite-dimensional trial and test (Galerkin) space of the form V h:" ' h V0, W h:" V0, (4.6) where ' h is a ane o set satisfying the essential BC of (4.1) and V0 h is a finite-dimensional .