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

1

Quadrature Integration in 1D Integration in 2D and 3D

2

Boundary conditions Penalty method Lagrange multiplier method

Institute of Structural Engineering Method of Finite Elements I 2

Quadrature Boundary conditions Integration in 1D Integration in 2D and 3D

Outline

1

Quadrature Integration in 1D Integration in 2D and 3D

2

Boundary conditions Penalty method Lagrange multiplier method

Institute of Structural Engineering Method of Finite Elements I 3

Quadrature Boundary conditions Integration in 1D Integration in 2D and 3D

Integration

Basic numerical integration Riemann sum Trapezoidal rule Comments: Both simple but inefficient! useful techniques in other contexts we use many more points than needed in order to achieve good accuracy

Institute of Structural Engineering Method of Finite Elements I 4

Quadrature Boundary conditions Integration in 1D 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 N1(x) = x2 N2(x) = x3 N3(x) = x4 · · · Instead of the Hermite basis functions, and discretize only displacements! Strong form: d2 dx2

• EI(x) d2w

dx2

• = f(x)

multiply with test function u(x) and integrate by parts twice a 1 EI(x) d2w(x) dx2 d2u(x) dx2 dx = 1 f(x)u(x)dx

aboundary terms [·]1 0 zero by basis selection Institute of Structural Engineering Method of Finite Elements I 5

Quadrature Boundary conditions Integration in 1D 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 Boundary conditions Integration in 1D Integration in 2D and 3D

Numerical Integration in 1D results

1

syms x % Coordinate a x i s

2

syms q ( x )

3

% The l o a d i n g

• f

the beam :

4

q ( x ) = x .ˆ2 − 0.3∗ x .ˆ4

5

NBASIS = 4 ;

6

f o r i = 1 : NBASIS

7

Phi ( i ) = x . ˆ ( i +1) ;

8

end

9

% 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 :

10

f o r i =1:NBASIS

11

f o r j = 1: NBASIS

12

K( i , j )= . . .

13

i n t ( EI∗ d i f f ( Phi ( i ) ,2) ∗ d i f f ( Phi ( j ) ,2) ,0 , L) ;

14

end

15

% A n a l y t i c a l computation

• f

the load v e c t o r :

16

f ( i ) = i n t ( q ( x ) ∗ Phi ( i ) ,0 , L) ;

17

end

Institute of Structural Engineering Method of Finite Elements I 7

Quadrature Boundary conditions Integration in 1D Integration in 2D and 3D

Numerical integration for LHS

1

q u a d r o r d e r = QUAD ORD;

2 3

%Loading the Gauss−Legendre quadrature nodes

4

% and weights f o r a f i l e :

5

nodes weights = W { quadr order −1};

6

nodes = ( nodes weights . nodes+1)∗L /2;

7

weights = W { quadr order −1}. weights ;

8 9

f o r i =1: N b a s i s

10

f o r j =1: N b a s i s

11

V = e v a l ( subs ( EI ∗ d i f f ( Phi ( i ) ,2) ∗ d i f f ( Phi ( j ) ,2) , nodes ) ) ;

12

K approx ( i , j ) = V’ ∗ weights ;

13

% 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 !

14

% we can get EXACT r e s u l t s without doing so

15

end

16 17

f a p p r o x ( i ) = e v a l ( subs ( Phi ( i )∗q ( x ) , nodes ) ) ’∗ weights ;

18

end

Institute of Structural Engineering Method of Finite Elements I 8

Quadrature Boundary conditions Integration in 1D Integration in 2D and 3D

Numerical Results

Institute of Structural Engineering Method of Finite Elements I 9

Quadrature Boundary conditions Integration in 1D Integration in 2D and 3D

Conclussions

Notice the errors in the approximation (1.7e−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 Boundary conditions Integration in 1D Integration in 2D and 3D

Optimal points and weights integration for polynomials:

Quadrature - Optimal weights given the nodes For a polynomial f(x) = a0 + a1x + · · · + anxn, the integral is approximated as 1

−1

f(x)dx ≃

N

• i=1

wif(xi) Assume xi are given. Also, for simplicity, f(x) = a0 + a1x + a2x2 + a3x3. Then, a0 1

−1

dx + a1 1

−1

xdx + a2 1

−1

x2dx + a3 1

−1

x3dx =

4

• i=1

wi(a0 + a1xi + a2x2

i + a3x3 i )

We equate all the terms, for every ai term. we arrive at the following linear system:     1 1 1 1 x1 x2 x3 x4 x2

1

x2

2

x2

3

x2

4

x3

1

x3

2

x3

3

x3

4

   

• Vandermonde matrix

       w1 w2 w3 w4        =          1

−1 dx

1

−1 xdx

1

−1 x2dx

1

−1 x3dx

        

Institute of Structural Engineering Method of Finite Elements I 11

Quadrature Boundary conditions Integration in 1D 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 Boundary conditions Integration in 1D 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 xi such that this integral is exact for every polynomial f(x) with degree n < 2N − 1. For a given N, solution xi are roots of the Legendre orthogonal polynomial of degree N

Institute of Structural Engineering Method of Finite Elements I 13

Quadrature Boundary conditions Integration in 1D Integration in 2D and 3D

Orthogonal polynomials Legendre orthogonal polynomials Lnsatisfy (by construction)a: Lk, Lm = 1

−1

Lk(x)Lm(x)dx = 2 2m + 1 δmk From now on this integral is going to be refered to as L2 inner product. Formulas for the first few polynomials areb, L0(x) = 1 L1(x) = x L2(x) = 1 2 (3x2 − 1) L3(x) = 1 8 (5x3 − 3x) L4(x) = 1 8 (35x4 − 30x2 + 3)

aThere are more general definitions of orthogonality and Gaussian quadrature where

f, gw =

Ω f(x)g(x)w(x)dx. The resulting polynomials are useful in other important contexts (stochastic

finite elements!)

bThey are computed with a process known as Gram-Schmidt orthogonalization Institute of Structural Engineering Method of Finite Elements I 14

Quadrature Boundary conditions Integration in 1D 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 2N − 1. Outline of proof:

1

consider a polynomial f2N−1(x) = a0 + a1x + · · · + a2N−1x2N−1

2

Apply polynomial divission with Legendre polynomial of degree N: f2N−1(x) = g(x)

• degree≤N−1

LN(x) + r(x)

• degree≤N−1

3

{x(N)

1

, · · · , x(N)

N

} are the N roots of the Legendre polynomial of degree N. For brevity, we

• mmit the superscript x(N)

i

= xi. Chose the quadrature nodes to be these exact points:

✘✘✘✘✘✘ ✘ ✿ 0

1

−1

g(x)Ln(x)dx + 1

−1

r(x)dx =✘✘✘✘✘✘✘

✘ ✿ 0

N

• i=1

wig(xi)Ln(xi) +

N

• i=1

wir(xi)

4

Since Ln(xi) = 0, the first term of the RHS vanishes. Moreover, since the degree of g(x) is 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

Quadrature Boundary conditions Integration in 1D Integration in 2D and 3D

Considering again that we are evaluating the quadrature rule on the roots of Legendre polynomial of degree N, the term LN becomes zero f2N−1(xi) = g(xi)✘✘✘

✘ ✿ 0

LN(xi) + r(xi) = r(xi) Finally we are left with 1

−1

f2N−1(x)dx = 1

−1

r(x) =

N

• i=1

wif2N−1(xi) And we are left with the task of determining the weights as in the non-gaussian case. Therefore, surprizingly, a N-point gauss quadrature rule is exact for polynomials up to degree 2N − 1. In practice, the computational determination of weights wi and nodes xi of 1D Gaussian quadrature rules happens through the so-called Golub-Welsch algorithm.

Institute of Structural Engineering Method of Finite Elements I 16

Quadrature Boundary conditions Integration in 1D Integration in 2D and 3D

Quadrilaterals

Tensor-product Quadrature Take tensor product of the 1D rules. Example: Quadrature in x: w(ξ)

i

= {w(ξ)

1

, w(ξ)

2

, w(ξ)

3

, w(ξ)

4

} ξi = {ξ1, ξ2, ξ3, ξ4} Quadrature in y w(η)

i

= {w(η)

1

, w(η)

2

, w(η)

3

, x(η)

4

} ηi = {η1, η2, η3, η4} Quadrature rule for quadrilateral: 1

−1

1

−1

f(ξ, η)dxdy =

4

• i=1

4

• j=1

w(η)

i

w(ξ)

j

f(ξi, ηj)

Institute of Structural Engineering Method of Finite Elements I 17

Quadrature Boundary conditions Penalty method Lagrange multiplier method

Outline

1

Quadrature Integration in 1D Integration in 2D and 3D

2

Boundary conditions Penalty method Lagrange multiplier method

Institute of Structural Engineering Method of Finite Elements I 18

Quadrature Boundary conditions Penalty method Lagrange multiplier method

Homogeneous constraints

Linear system from FE discretization Consider a linear system from a FE discretization: [K]{u} = {f}     K11 . . . K1N . . . ... . . . KN1 . . . KNN            u1 . . . uN        =        f1 . . . fN        Homogeneous constraints the constraint un = 0, in a matrix form consistent with our FE linear system reads,

N

• i=0

viui = 0

• . . .

1 . . .                  u1 . . . un . . . uN                  = 0

Institute of Structural Engineering Method of Finite Elements I 19

Quadrature Boundary conditions Penalty method Lagrange multiplier method

Imposing homogeneous constr. with the Penalty method

Imposing the constraint The constraint vector should be in a form consistent with the matrix (a bilinear form that corresponds to a matrix). In order to impose the constraint, we simply multiply with a large number p >> tr([K]) and add that to our linear system. [V] = Vij = vivj [V] = {v}T {v} [V] =            . . . ... . . . 1 . . . ... . . .            This forces the degrees of freedom in the system to satisfy the constraint numerically (not explicitly).

Institute of Structural Engineering Method of Finite Elements I 20

Quadrature Boundary conditions Penalty method Lagrange multiplier method

Adding to the linear system

Adding the constraint Simply set, [K] ← [K]pen with [K]pen = [K] + [V] · p with p >> tr([K]) and solve {u} = [K]−1

pen{f}

Inhomogeneous constraints and MFC: What if, we need to set a displacement to a given value (inhomogeneous constraint), or a displacement the same as another displacement (multi-DOF constraints)?

Institute of Structural Engineering Method of Finite Elements I 21

Quadrature Boundary conditions Penalty method Lagrange multiplier method

Multi- degree of freedom constraints

Equation for a multi-DOF constraint Consider we have a constraint of the form 2u5 + 3u7 − 2u10 = 0. Again, written as an equation,

• . . .

2 3 −2 . . .                                      u1 . . . u5 u6 u7 u8 u9 u10 . . . uN                                      = 0

Institute of Structural Engineering Method of Finite Elements I 22

Quadrature Boundary conditions Penalty method Lagrange multiplier method

Penalty, multi-DOF constraints

Application of the constraints Not much change when being systematic about the application! [V] ={v}T {v} [V] =                          . . . ... 4 6 −4 6 9 −6 . . . . . . −4 −6 4 ... . . .                          The rest are the same! (multiply by a large number and add to [K] before solution).

Institute of Structural Engineering Method of Finite Elements I 23

Quadrature Boundary conditions Penalty method Lagrange multiplier method

Non-homogeneous constraints with the Penalty method and multiple constraints

Non-zero RHS Say, we want to constrain u3 − 2u5 = 3 and at the same time u2 − u3 = 0. In matrix form, the constraints read:

• 1

−2 . . .

• {v1}

                 u1 u2 u3 u4 u5 . . .                  = 3

• a1
• 1

−1 . . .

• {v2}

                 u1 u2 u3 u4 u5 . . .                  =

• a2

Institute of Structural Engineering Method of Finite Elements I 24

Quadrature Boundary conditions Penalty method Lagrange multiplier method

Application of BC and solution

Application of BC on system In order to impose the constraint on the system, we proceed as such: Multiplying both equations with the corresponding {vi}T and summing up the constraints for the left-hand side we get, [V ] =

2

• i=1

{vi}T {vi} and for the right-hand-side: {A} =

2

• i=1

{vi}T ai [Kpen] = [K] + [V ] · p {fpen} = {f} + {A} · p What changes, is that we also have a right-hand-side due to the non-homogeneous constraints! Solution Simply compute {u} = [Kpen]−1{fpen}

Institute of Structural Engineering Method of Finite Elements I 25

Quadrature Boundary conditions Penalty method Lagrange multiplier method

Lagrange multipliers

constraints for LHS Another method, just as general and simple to implement, is the Lagrange multiplier method: given the constraint matrix [v] from the previous slide, set [KLagr] =

• [K]

[v]T [v]

• [KLagr] =

       K11 . . . K1N v11 v21 . . . ... . . . . . . . . . KN1 . . . KNN v1N v2N v11 . . . v1N v21 . . . v2N        Constraints at RHS Instead of adding, we are again appending. {fLagr} =

• {f}

{P }

• Where, {P } is a vector containing the values of the non-homogeneous constraints. From the

previous example, {P } =

• 3
• .

Institute of Structural Engineering Method of Finite Elements I 26

Quadrature Boundary conditions Penalty method Lagrange multiplier method

Comparisson

The clear advantage of Lagr. Multiplier method over penalty is accuracy. Penalty has the advantage that does not introduce degrees of freedom to the linear system. However, this rarely becomes a real problem with modern hardware. Some linear system solvers rely on the matrix being positive definite for fast convergence - Lagrange multipliers may not work with some linear system solvers Lagrange multiplier method provides directly the constraint forces (as the values corresponding to the rows and columns of the constraints). When constraints are linearly dependent (a constraint vector can be expressed as linear combination of other constraints in the same system) iterative solvers (i.e. GMRES) have to be used. Linear independence of constraints is hard to judge in some cases (however, there are linear algebra techniques one can perform to remedy that). The multi-freedom homogeneous case for the Master-slave elimination method, is described partially in the project description. This method is also exact (as Lagr. Multipliers).

Institute of Structural Engineering Method of Finite Elements I 27