Chapter 9 The Finite Element Method for 2D elliptic PDEs

1 IiiiiiChapter 9 The Finite ElementMethod for 2D ellipticPDEsThe procedure of the Finite Element Method to solve 2D problems is the same as that for1D problems, as the flow chart below Integration by parts weak form inV:a(u,v) =L(v)or minv VF(v) Vh( Finite dimensional space and basis functions) a(uh,vh) =L(vh) uhand error The second Green s theorem and integration byparts in 2 DLet us first recall the 2D version of the well known divergencetheorem in Cartesian H1( ) H1( )is a vector in 2D, thenZ Z Fdxdy=Z F nds,( )wherenis the unit normal direction pointing outward at the boundary with line elementds, and is the gradient operator once again is = [ x, y] second Green s theorem is a corollary of the divergence theorem if we setF=v u= v u x,v u y since F= x v u x + y v u y = u x v x+v 2u x2+ u y v y+v 2u y2= u v+v u,219iiiiii220 Chapter 9.

2 The Finite Element Method for 2D elliptic PDEs nnFigure diagram of a two dimensional domain , its boundary and its unit normal u= u=uxx+uyy, we obtainZ Z Fdxdy=Z Z ( u v+v u)dxdy=Z F nds=Z v u nds=Z v u nds,wheren= (nx,ny) (n2x+n2y= 1) is the unit normal direction, and u n= u n=nx u x+ny u y, the normal derivative derivative ofu, see Fig. for an illustration. Thisresult immediately yields the formula for integration by parts in (x,y) H2( )andv(x,y) H1( )where is a bounded domain,thenZ Z v u dxdy=Z v u nds Z Z u v dxdy.( )Note: the normal derivative u/ nis sometimes written more concisely important elliptic PDEs in 2D Cartesian coordinates are:uxx+uyy= 0,Laplace equation, uxx uyy=f(x,y),Poisson equation, uxx uyy+ u=f,generalized Helmholtz equation,uxxxx+ 2uxxyy+uyyyy= 0,Bi-harmonic >0, the generalized Helmholtz equation is easier to solve than when < , the expressions involved in these PDEs may also be abbreviated using thegradient operator , ,uxx+uyy= u= uas mentioned before.

3 We also recallthat a general linear second order elliptic PDE has the forma(x,y)uxx+ 2b(x,y)uxy+c(x,y)uyy+d(x,y)ux+e(x,y)uy+g (x,y)u=f(x,y) Weak form of second order self-adjoint elliptic PDEs221with discriminantb2 ac <0. A second order self-adjoint elliptic partial differentialequation has the form (p(x,y) u) +q(x,y)u=f(x,y).( ) Boundary conditionsIn 2D, the domain boundary is one or several curves. We consider the following variouslinear boundary conditions. Dirichlet boundary condition on the entire boundary, ,u(x,y)| =u0(x,y) isgiven. Neumann boundary condition on the entire boundary, , u/ n| =g(x,y) this case, the solution to a Poisson equation may not be unique or even exist, de-pending upon whether a compatibility condition is satisfied. Integrating the Poissonequation over the domain, we haveZ Z fdxdy= Z Z udxdy= Z Z udxdy= Z unds= Z g(x,y)ds,( )which is the compatibility condition to be satisfied for the solution to exist.

4 If asolution does exist, it is not unique as it is determined within an arbitrary constant. Mixed boundary condition on the entire boundary, , (x,y)u(x,y) + (x,y) u n= (x,y)is given, where (x,y), (x,y), and (x,y) are known functions. Dirichlet, Neumann, and Mixed boundary conditions on some parts of the Weak form of second order self-adjoint ellipticPDEsNow we derive the weak form of the self-adjoint PDE ( ) witha homogeneous Dirichletboundary condition on part of the boundary D,u| D= 0 and a homogeneous Neumannboundary condition on the rest of boundary N= D, u n| N= 0. Multiplyingthe equation ( ) by a test functionv(x,y) H1( ), we haveZ Z n (p(x,y) u) +q(x,y)uov dxdy=Z Z fv dxdy;and on using the formula for integration by parts the left-hand side becomesZ Z p u v+quv dxdy Z pvunds,iiiiii222 Chapter 9. The Finite Element Method for 2D elliptic PDEsso the weak form isZ Z (p u v+quv)dxdy=Z Z fvdxdy+Z Npg(x,y)v(x,y)ds v(x,y) H1( ).

5 ( )Here Nis the part of boundary where a Neumann boundary condition isapplied; andthe solution space resides inV= v(x,y), v(x,y) = 0,(x,y) D, v(x,y) H1( ) ,( )where Dis the part of boundary where a Dirichlet boundary conditionis Verification of conditions of the Lax-Milgram LemmaThe bilinear form for ( ) isa(u,v) =Z Z (p u v+quv)dxdy,( )and the linear form isL(v) =Z Z fv dxdy( )for a Dirichlet BC on the entire boundary. As before, we assume that0< pmin p(x,y) pmax,0 q(x) qmax, p C( ), q C( ).We need the Poincar einequality to prove the V- elliptic (x,y) H10( ), R2, ,v(x,y) H1( )and vanishes at theboundary (can be relaxed to a point on the boundary), thenZ Z v2dxdy CZZ | v|2dxdy,( )whereCis a we are ready to check the conditions of the Lax-Milgram It is obvious thata(u,v) =a(v,u).2. It is easy to see that|a(u,v)| max{pmax,qmax} Z Z (| u v|+|uv|)dxdy = max{pmax,qmax}|(|u|,|v|)1| max{pmax,qmax}kuk1kvk1,soa(u,v) is a continuous and bounded bilinear Triangulation and basis functions2233.

6 From the Poincar einequality|a(v,v)|= Z Z p | v|2+qv2 dxdy pminZ Z | v|2dxdy=12pminZ Z | v|2dxdy+12pminZ Z | v|2dxdy 12pminZ Z | v|2dxdy+pmin2CZ Z |v|2dxdy 12pminminn1,1 Cokvk21,thereforea(u,v) is Finally, we show thatL(v) is continuous:|L(v)|=|(f,v)0| kfk0kvk0 , the solutions to the weak form and the minimization form are unique andbounded inH10( ). Triangulation and basis functionsThe general procedure of the Finite Element Method is the samefor any dimension, andthe Galerkin Finite Element Method involves the following main steps. Generate a triangulation over the domain. Usually the triangulation is composedof either triangles or rectangles. There are a number of meshgeneration softwarepackages available, , the Matlab PDE toolbox from Mathworks, Triangle fromCarnegie Mellon University,etc. Some are available through the Internet. Construct basis functions over the triangulation.

7 We mainly consider the conformingfinite Element Method in this book. Assemble the stiffness matrix and the load vector Element by Element , using eitherthe Galerkin Finite Method (the weak form) or the Ritz Finite Method (the minimiza-tion form). Solve the system of equations. Do the error Fig. , we show a diagram of simple mesh generation process. The circulardomain is approximated by a polygon with five vertices (selected points on the boundary).We then connect the five vertices to get initial five triangles(solid line) to obtain an initialcoarse mesh. We can refine the mesh using the so called middle point rule by connectingall the middle points of all triangles in the initial mesh to obtain a finer mesh (solid anddashed lines).iiiiii224 Chapter 9. The Finite Element Method for 2D elliptic PDEsFigure diagram of a simple generation process and the middle Triangulation and mesh parametersGiven a general domain, we can approximate the domain by a polygon and then generatea triangulation over the polygon, and we can refine the triangulation if necessary.

8 A simpleapproach is the mid-point rule by connecting all the middle points of three sides of existingtriangles to get a refined triangulation usually has the mesh parameters p: polygonal region =K1 K2 K3 Knelem,Kj: are non-overlapping triangles, j= 1,2, , nelem,Ni: are nodal points, i= 1,2, , nnode,hj: the longest side ofKj, j: the diameter of the circle inscribed inKj(encircle),h: the largest of allhj, h= max{hj}, : the smallest of all j, = min{ j},with1 jhj >0,where the constant is a measurement of the triangulation quality, see Fig. for anillustration of such a s andh s. The larger the , the better the quality of the triangula-tion. Given a triangulation, a node is also the vertex of all adjacent triangles. We do notdiscuss hanging nodes The FE space of piecewise linear functions over atriangulationFor linear second order elliptic PDEs, we know that the solution space is in theH1( ).

9 Unlike the 1D case, an elementv(x,y) inH1( ) may not be continuous under the Sobolevembedding theorem. However, in practice most solutions areindeed continuous, especiallyfor second order PDEs with certain regularities. Thus, we still look for a solution in Triangulation and basis functions225continuous function spaceC0( ). Let us first consider how to construct piecewise linearfunctions over a triangulation with the Dirichlet BCu(x,y)| = a triangulation, we defineVh=nv(x,y) is continuous in and piecewise linear over eachKj,v(x,y)| = 0}.( )We need to determine the dimension of this space and construct a set of basis each triangle, a linear function has the formvh(x,y) = + x+ y ,( )where , and are constants (three free parameters). LetPk={p(x,y),a polynomial of degree ofk}.( )We have the following A linear functionp1(x,y) = + x+ ydefined on a triangle is uniquely determinedby its values at the three Ifp1(x,y) P1andp2(x,y) P1are such thatp1(A) =p2(A)andp1(B) =p2(B),whereAandBare two points in thexy-plane, thenp1(x,y) p2(x,y), (x,y) IAB, whereIABis the line segment betweenAandB.

10 A1a2p1p2a3 Figure diagram of a triangle with three verticesa1,a2, anda3; anadjacent triangle with a common side; and the local coordinate system in whicha2is the origin anda2a3is the 9. The Finite Element Method for 2D elliptic PDEsProof:Assume the vertices of the triangle are (xi,yi),i= 1,2,3. The linear functiontakes the valueviat the vertices, ,p(xi,yi) =vi,so we have the three equations + x1+ y1=v1, + x2+ y2=v2, + x3+ y3= determinant of this linear algebraic system isdet 1x1y11x2y21x3y3 = 2 area of the triangle6= 0 since jhj >0,( )hence the linear system of equations has a unique let us prove the second part of the theorem. Suppose that the equation of theline segment isl1x+l2y+l3= 0, l21+l226= can solve forxor fory:x= l2y+l3l1ifl16= 0,ory= l1x+l3l2ifl26= loss of generality, let us assumel26= 0 such thatp1(x,y) = + x+ y= + x l1x+l3l2 = l3l2 + l1l2 x= 1+ , we havep2(x,y) = 1+ (A) =p2(A) andp1(B) =p2(B), 1+ 1x1=p(A), 1+ 1x1=p(A), 1+ 1x2=p(B), 1+ 1x2=p(B), Triangulation and basis functions227where both of the linear system of algebraic equations have the same coefficient matrix"1x11x2#that is non-singular sincex16=x2(because pointsAandBare distinct).