Method of Finite Elements I PhD Candidate - Charilaos Mylonas HIL - PowerPoint PPT Presentation

Method of Finite Elements I PhD Candidate - Charilaos Mylonas HIL H33.1 Quadrature and Boundary Conditions , 26 March, 2018 Institute of Structural Engineering Method of Finite Elements I 1

Quadrature Boundary conditions Outline Quadrature 1 Integration in 1D Integration in 2D and 3D Boundary conditions 2 Penalty method Lagrange multiplier method Institute of Structural Engineering Method of Finite Elements I 2

Quadrature Integration in 1D Boundary conditions Integration in 2D and 3D Outline Quadrature 1 Integration in 1D Integration in 2D and 3D Boundary conditions 2 Penalty method Lagrange multiplier method Institute of Structural Engineering Method of Finite Elements I 3

Quadrature Integration in 1D Boundary conditions Integration in 2D and 3D Integration Basic numerical integration Riemann sum Comments: Both simple but inefficient! useful techniques in other contexts we use many more points than needed Trapezoidal rule in order to achieve good accuracy Institute of Structural Engineering Method of Finite Elements I 4

Quadrature Integration in 1D Boundary conditions Integration in 2D and 3D Numerical Integration in 1D example: Euler-Bernoulli beam We are going to apply the Galerkin method to the Euler Bernoulli beam, but with the basis N 1 ( x ) = x 2 N 2 ( x ) = x 3 N 3 ( x ) = x 4 · · · Instead of the Hermite basis functions, and discretize only displacements! Strong form: d 2 EI ( x ) d 2 w � � = f ( x ) dx 2 dx 2 multiply with test function u ( x ) and integrate by parts twice a � 1 � 1 EI ( x ) d 2 w ( x ) d 2 u ( x ) dx = f ( x ) u ( x ) dx dx 2 dx 2 0 0 a boundary terms [ · ] 1 0 zero by basis selection Institute of Structural Engineering Method of Finite Elements I 5

Quadrature Integration in 1D Boundary conditions Integration in 2D and 3D Numerical integration - Galerkin method code The full commented code of the following demo will be made available to you in case you would like to play around. This code is a comparisson between analytical and numerical evaluation of integrals needed in the implementation of the Galerkin procedure. In the following slides the most important parts of the code are highlighted. Namely: The definition of the load the assembly of the system with analytical integration (symbolic integration) the assembly of the system with numerical integration (Gauss-Legendre) Institute of Structural Engineering Method of Finite Elements I 6

Quadrature Integration in 1D Boundary conditions Integration in 2D and 3D Numerical Integration in 1D results syms x % Coordinate a x i s 1 syms q ( x ) 2 % The l o a d i n g of the beam : 3 q ( x ) = x .ˆ2 − 0.3 ∗ x .ˆ4 4 NBASIS = 4 ; 5 f o r i = 1 : NBASIS 6 Phi ( i ) = x . ˆ ( i +1) ; 7 end 8 % Compute the i n t e g r a l s needed f o r the s t i f f n e s s matrix : 9 f o r i =1:NBASIS 10 f o r j = 1: NBASIS 11 K( i , j )= . . . 12 i n t ( EI ∗ d i f f ( Phi ( i ) ,2) ∗ d i f f ( Phi ( j ) ,2) ,0 , L) ; 13 end 14 % A n a l y t i c a l computation of the load v e c t o r : 15 f ( i ) = i n t ( q ( x ) ∗ Phi ( i ) ,0 , L) ; 16 end 17 Institute of Structural Engineering Method of Finite Elements I 7

Quadrature Integration in 1D Boundary conditions Integration in 2D and 3D Numerical integration for LHS q u a d r o r d e r = QUAD ORD; 1 2 %Loading the Gauss − Legendre quadrature nodes 3 % and weights f o r a f i l e : 4 nodes weights = W { quadr order − 1 } ; 5 nodes = ( nodes weights . nodes+1) ∗ L /2; 6 weights = W { quadr order − 1 } . weights ; 7 8 f o r i =1: N b a s i s 9 f o r j =1: N b a s i s 10 V = e v a l ( subs ( EI ∗ d i f f ( Phi ( i ) ,2) ∗ d i f f ( Phi ( j ) ,2) , nodes ) ) ; 11 K approx ( i , j ) = V’ ∗ weights ; 12 % Note that we don ’ t i n t e g r a t e a n a l y t i c a l l y ! 13 % we can get EXACT r e s u l t s without doing so 14 end 15 16 f a p p r o x ( i ) = e v a l ( subs ( Phi ( i ) ∗ q ( x ) , nodes ) ) ’ ∗ weights ; 17 end 18 Institute of Structural Engineering Method of Finite Elements I 8

Quadrature Integration in 1D Boundary conditions Integration in 2D and 3D Numerical Results Institute of Structural Engineering Method of Finite Elements I 9

Quadrature Integration in 1D Boundary conditions Integration in 2D and 3D Conclussions Notice the errors in the approximation ( 1 . 7 e − 13 for quadrature of order 4!). This is close to machine precission (the ”eps” command in matlab). 1D Quadrature order ”N” means we evaluate the function inside the integral ”N” times, and we sum the results according to some weight (see next slide). Instead of integrating analytically a function, we can get EXACT results by carefully selecting points of evaluation (instead of regular intervals as in trapezoidal and Riemann sums). The accuracy of the approximation depends on the Gaussian quadrature order Institute of Structural Engineering Method of Finite Elements I 10

Quadrature Integration in 1D Boundary conditions Integration in 2D and 3D Optimal points and weights integration for polynomials: Quadrature - Optimal weights given the nodes For a polynomial f ( x ) = a 0 + a 1 x + · · · + a n x n , the integral is approximated as � 1 N � f ( x ) dx ≃ w i f ( x i ) − 1 i =1 Assume x i are given. Also, for simplicity, f ( x ) = a 0 + a 1 x + a 2 x 2 + a 3 x 3 . Then, � 1 � 1 � 1 � 1 x 2 dx + a 3 x 3 dx = a 0 dx + a 1 xdx + a 2 − 1 − 1 − 1 − 1 4 � w i ( a 0 + a 1 x i + a 2 x 2 i + a 3 x 3 i ) i =1 We equate all the terms, for every a i term. we arrive at the following linear system: � 1       − 1 dx 1 1 1 1 w 1      � 1          − 1 xdx x 1 x 2 x 3 x 4 w 2   =   � 1 x 2 x 2 x 2 x 2 − 1 x 2 dx w 3   1 2 3 4         � 1 x 3 x 3 x 3 x 3     w 4 − 1 x 3 dx   1 2 3 4 � �� � Vandermonde matrix Institute of Structural Engineering Method of Finite Elements I 11

Quadrature Integration in 1D Boundary conditions Integration in 2D and 3D Non-Gaussian Quadrature Numerical computation of weights N-point quadrature, extact for polynomials up to degree N-1 Weights can be negative, in practice they have values of very different scales and solution of the system for high order quadrature is numerically unstable! The technique will approximate well integrals of functions that are approximated well with the corresponding polynomials. Lower order quadrature works, but as we will see it is far from the best we can do! (half as good). Institute of Structural Engineering Method of Finite Elements I 12

Quadrature Integration in 1D Boundary conditions Integration in 2D and 3D Legendre polynomials In the following, we will see how the Legendre polynomials are important for quadrature in 1D. We are to use x i such that this integral is exact for every polynomial f ( x ) with degree n < 2 N − 1 . For a given N , solution x i are roots of the Legendre orthogonal polynomial of degree N Institute of Structural Engineering Method of Finite Elements I 13

Quadrature Integration in 1D Boundary conditions Integration in 2D and 3D Orthogonal polynomials Legendre orthogonal polynomials L n satisfy (by construction) a : � 1 2 � L k , L m � = L k ( x ) L m ( x ) dx = 2 m + 1 δ mk − 1 From now on this integral is going to be refered to as L 2 inner product . Formulas for the first few polynomials are b , L 0 ( x ) = 1 L 1 ( x ) = x L 2 ( x ) = 1 2 (3 x 2 − 1) L 3 ( x ) = 1 8 (5 x 3 − 3 x ) L 4 ( x ) = 1 8 (35 x 4 − 30 x 2 + 3) a There are more general definitions of orthogonality and Gaussian quadrature where � f, g � w = � Ω f ( x ) g ( x ) w ( x ) dx . The resulting polynomials are useful in other important contexts ( stochastic finite elements!) b They are computed with a process known as Gram-Schmidt orthogonalization Institute of Structural Engineering Method of Finite Elements I 14

Quadrature Integration in 1D Boundary conditions Integration in 2D and 3D Legendre orthogonal polynomials and quadrature The roots of Legendre polynomials of degree N , used as integration nodes, yield the exact integral for f polynomial of maximum degree 2 N − 1 . Outline of proof: consider a polynomial f 2 N − 1 ( x ) = a 0 + a 1 x + · · · + a 2 N − 1 x 2 N − 1 1 Apply polynomial divission with Legendre polynomial of degree N : 2 f 2 N − 1 ( x ) = g ( x ) L N ( x ) + r ( x ) ���� ���� degree ≤ N − 1 degree ≤ N − 1 { x ( N ) , · · · , x ( N ) } are the N roots of the Legendre polynomial of degree N . For brevity, we 3 1 N ommit the superscript x ( N ) = x i . Chose the quadrature nodes to be these exact points: i ✿ 0 ✿ 0 � 1 � 1 N N r ( x ) dx = ✘✘✘✘✘✘✘ ✘ ✘ � � ✘✘✘✘✘✘ g ( x ) L n ( x ) dx + w i g ( x i ) L n ( x i ) + w i r ( x i ) − 1 − 1 i =1 i =1 Since L n ( x i ) = 0 , the first term of the RHS vanishes. Moreover, since the degree of g ( x ) is 4 lower than N , the first integral term on the left is also zero! That is, through orthogonality and the abillity to express polynomials of order N − 1 with a polynomial basis that contains all the monomials (such as the Legendre basis order N )! Institute of Structural Engineering Method of Finite Elements I 15