Galerkin finite element method Boundary value problem   Lu = f      u=g 0

 n · ∇u = g1      n · ∇u + αu = g 2

→

weighted residual formulation

in Ω

partial differential equation

on Γ0

Dirichlet boundary condition

on Γ1

Neumann boundary condition

on Γ2

Robin boundary condition

1. Multiply the residual of the PDE by a weighting function w vanishing on the Dirichlet boundary Γ0 and set the integral over Ω equal to zero 2. Integrate by parts using the Neumann and Robin boundary conditions 3. Represent the approximate solution uh ≈ u as a linear combination of polynomial basis functions ϕi defined on a given mesh (triangulation)

4. Substitute the functions uh and ϕi for u and w in the weak formulation 5. Solve the resulting algebraic system for the vector of nodal values ui

Construction of 1D finite elements 1

1. Linear finite elements

2

Consider the barycentric coordinates x2 − x λ1 (x) = , x2 − x1

x − x1 λ2 (x) = x2 − x1

defined on the element e = [x1 , x2 ] • λi ∈ P1 (e),

=

=

Basis functions

x2

'1

'2

2 − dλ dx

e

i, j = 1, 2

• λ1 (x) + λ2 (x) = 1, 1 − x2 −x 1

x1

i = 1, 2

• λi (xj ) = δij , dλ1 dx

e

∀x ∈ e

x1

x2 u2

constant derivatives

ϕ1 |e = λ1 ,

uh (x) = u1 ϕ1 (x) + u2 ϕ2 (x),

u1

ϕ2 |e = λ2 ∀x ∈ e

e x1

x2

Construction of 1D finite elements 1

2. Quadratic finite elements {λ1 (x), λ2 (x)} x1 = {1, 0}, x12 =

x1 +x2 2

barycentric coordinates e

x2 = {0, 1} = { 21 , 21 }

Basis functions

2

endpoints

midpoint

x1

x2

'1

'12

ϕ1 , ϕ2 , ϕ12 ∈ P2 (e)

ϕ1 (x) = λ1 (x)(2λ1 (x) − 1)

e x12

x1

ϕ2 (x) = λ2 (x)(2λ2 (x) − 1) ϕ12 (x) = 4λ1 (x)λ2 (x)

x2 u2

u1

Shape function uh |e uh (x) = u1 ϕ1 (x) + u2 ϕ2 (x) + u12 ϕ12 (x)

'2

u12

e x1

x12

x2

Construction of 1D finite elements 3. Cubic finite elements {λ1 (x), λ2 (x)}

1

barycentric coordinates

x1 = {1, 0},

x12 =

2x1 +x2 3

= { 32 , 13 }

x2 = {0, 1},

x21 =

x1 +2x2 3

= { 31 , 23 }

Basis functions

e x1

x2

ϕ1 , ϕ2 , ϕ12 , ϕ21 ∈ P3 (e)

ϕ1 (x) = 21 λ1 (x)(3λ1 (x) − 2)(3λ1 (x) − 1) ϕ2 (x) = 12 λ2 (x)(3λ2 (x) − 2)(3λ2 (x) − 1) ϕ12 (x) = 92 λ1 (x)λ2 (x)(3λ1 (x) − 1) ϕ21 (x) = 92 λ1 (x)λ2 (x)(3λ2 (x) − 1) Shape function

2

u2 u12 u1

u21

e x1

x12 x21 x2

uh (x) = u1 ϕ1 (x) + u2 ϕ2 (x) + u12 ϕ12 (x) + u21 ϕ21 (x)

Construction of triangular finite elements x3

Let x = (x, y) = {λ1 (x), λ2 (x), λ3 (x)}

e

Construct 2D barycentric coordinates λi ∈ P1 (e),

λi (xj ) = δij ,

x

i, j = 1, 2, 3

x1 Polynomial fitting: λi (x) = ci1 + ci2 x + ci3 y       λ (x ) = c + c x + c y = δ i 1 i1 i2 1 i3 1 i1  c 1 x y  i1 1 1  1 x2 y2   ci2  =  λi (x2 ) = ci1 + ci2 x2 + ci3 y2 = δi2    ci3 1 x3 y3 λi (x3 ) = ci1 + ci2 x3 + ci3 y3 = δi3 {z } |

x2  δi1 δi2  δi3

A

We have 3 systems of 3 equations for 9 unknowns. They can be solved for the unknown coefficients cij by resorting to Cramer’s rule. det A = x2 y3 + x1 y2 + x3 y1 − x2 y1 − x3 y2 − x1 y3 Area of the triangle

|e| = 12 | det A|

(also needed for quadrature rules)

Construction of triangular finite elements Connect the point x to the vertices to construct the barycentric splitting

i = 1, 2, 3 3 S ei e= xi ,

x3 e e2 e1

i=1

Areas of the triangles   1 x y A1 =  1 x2 y2  , 1 x3 y3

1 2 | det Ai (x)|

|ei (x)| =   1 x1 y1 A2 =  1 x y  , 1 x3 y3

Solution of the linear systems:

λi (x) =

|ei (x)| |e|

x1

∀x ∈ e

x2



 1 x1 y1 A3 =  1 x2 y2  1 x y det Ai (x) = det A ,

It is obvious that the barycentric coordinates satisfy |e1 (x)| + |e2 (x)| + |e3 (x)| = |e|,

e3

⇒

i = 1, 2, 3

λi (xj ) = δij

λ1 (x) + λ2 (x) + λ3 (x) ≡ 1 e2

A similar interpretation is possible in one dimension: x1

e1 e

x2

Construction of triangular finite elements Nodal barycentric coordinates x1 = {1, 0, 0}, x12 = { 12 , 21 , 0},

x2 = {0, 1, 0},

ϕ1 = λ 1 ,

midpoint-oriented

ϕ12 = 1 − 2λ3 ,

2. Quadratic elements

ϕ3 = λ3 (2λ3 − 1),

x1

x12

x2

ϕ2 = λ2 ,

ϕ3 = λ3

(standard)

ϕ13 = 1 − 2λ2 ,

ϕ23 = 1 − 2λ1

uh (x) = c1 + c2 x + c3 y + c4 x2 + c5 xy + c6 y 2 ∈ P2 (e)

Standard P2 basis (6 nodes)

ϕ2 = λ2 (2λ2 − 1),

x23 = {0, 21 , 21 }

x123 x23

uh (x) = c1 + c2 x + c3 y ∈ P1 (e)

vertex-oriented

ϕ1 = λ1 (2λ1 − 1),

x13

x3 = {0, 0, 1}

x13 = { 12 , 0, 12 },

1. Linear elements

x3

x123 = { 13 , 31 , 31 }

ϕ12 = 4λ1 λ2 ϕ13 = 4λ1 λ3 ϕ23 = 4λ2 λ3

Extended P2+ basis (7 nodes) ϕi = λi (2λi − 1) + 3λ1 λ2 λ3

ϕij = 4λi λj − 12λ1 λ2 λ3

ϕ123 = 27λ1 λ2 λ3 ,

i, j = 1, 2, 3

Construction of tetrahedral finite elements x4

Let x = (x, y, z) = {λ1 (x), λ2 (x), λ3 (x), λ4 (x)} Construct

λi ∈ P1 (e) : λi (xj ) = δij ,

Polynomial fitting: 

1  1   1 1

x1 x2 x3 x4

y1 y2 y3 y4

x124

i, j = 1, . . . , 4









x134 x234

x3

x123

x1

λi (x) = ci1 + ci2 x + ci3 y + ci4 z

δi1 ci1 z1     z2    ci2  =  δi2  z3   ci3   δi3  δi4 ci4 z4

1. Linear elements

e

x2

Barycentric splitting

e=

4 S

ei

i=1

λi (x) =

|ei (x)| , |e|

i = 1, . . . , 4

uh (x) = c1 + c2 x + c3 y + c4 z ∈ P1 (e)

vertex-oriented

ϕ1 = λ 1 ,

ϕ2 = λ2 ,

face-oriented

ϕ123 = 1 − 3λ4 ,

ϕ134 = 1 − 3λ2 ,

ϕ3 = λ3 ,

ϕ124 = 1 − 3λ3

ϕ234 = 1 − 3λ1

ϕ4 = λ4

(standard)

xijk =

xi +xj +xk 3

2. Higher-order approximations are possible but rather expensive in 3D

Coordinate and element transformations Idea: define the basis functions on a geometrically simple reference element n = 1 n = 2n = 1 y^ y^

n=1

x^ 3 0

e^

1

x^1

0

x^2 x^1

x^ 3

1

x^ 3 x^14

1

e^0 e^10 e^1 0 1 e^ x^ 1 x^x^21 x^ x^ 2 x^x^12 x^

01 0

y^ Fe e x1

Fe

x1

x3Fe

e

x1

x2 e

e x2 x x1

Fe y

x2 x1

Linear mapping in Rn

z^

x^ 4

1

1 1 0

e^

x^

Fey z e

e

y^

x3

x1

e

x2x2

Fe : y eˆ −→ e

x3 x1

xx

x y

1

e^

1

x^ 2 x^

x^x^2 1 x^

1

y^ Fe 3 Fe x 4 z

x3 x1

01

x^ 2 x^x2^ x^x^ 1 1

x^ 3

x2

xx2 x

e^ e^

n=3

n=3 z^

x^ 4

1

x^x^21 x^x^1

1

y

n= 2 3 n=

n=2 y^ z^

x^ 3

Fe

Fe x4 z e x2

x4 x3 x1

x

y

e

x3

x2

x

Coordinate and element transformations Mapping of the reference element eˆ onto an element e with vertices xi Fe :

e = Fe (ˆ e)

x = Fe (ˆ x) =

n+1 X

xi λi (ˆ x),

∀x ∈ e

i=1

where λi are the barycentric coordinates. This mapping is of the form ˆ + be , Fe (ˆ x) = Be x

Be ∈ Rn×n ,

be ∈ Rn

It is applicable to arbitrary simplex elements with straight sides n=1 n=2

n=3

Be = x2 − x1 ,  x2 − x1 Be = y2 − y1  x2 − x1 Be =  y2 − y1 z2 − z1

be = x1  x3 − x1 , y3 − y1

be =

 x3 − x1 x4 − x1 y3 − y1 y4 − y1  , z3 − z1 z4 − z1



x1 y1

  x1 be =  y1  z1 

Coordinate and element transformations Properties of the linear mapping • vertices are mapped onto vertices

ˆ + be Fe (ˆ x) = Be x xi = Fe (ˆ xi )

• midpoints of sides are mapped onto midpoints of sides   ˆ i +ˆ xi +xj x xj xij = 2 = Fe = Fe (ˆ xij ) 2 • barycenters are mapped onto barycenters   ˆ i +ˆ xi +xj +xk x xj +ˆ xk = Fe (ˆ xijk ) = Fe xijk = 3 3 The values of ϕi on the physical element e are defined by the formula ϕi (x) = ϕˆi (Fe−1 (x)),

∀x ∈ e

ϕi (x) = ϕˆi (ˆ x),

x = Fe (ˆ x)

Note that ϕi (xj ) = ϕˆi (ˆ xj ) = δij and the degree of basis functions (linear, ˆ quadratic, cubic etc.) is preserved since x depends linearly on x

Coordinate and element transformations Derivative transformations

ϕi (x) = ϕˆi (ˆ x),

ˆ ϕˆi = J∇ϕi ∇

Chain rule

x ∈ e,

ˆ ∈ eˆ x

where J is the Jacobian of the (inverse)

mapping as introduced before in the context of the finite difference method ∂ϕi ∂x ∂ϕi ∂y

= =

1 det J 1 det J

h h

∂ϕ ˆi ∂y ∂x ˆ ∂ yˆ ∂ϕ ˆi ∂x ∂ yˆ ∂ x ˆ

−

∂ϕ ˆi ∂y ∂ yˆ ∂ x ˆ

−

∂ϕ ˆi ∂x ∂x ˆ ∂ yˆ

i

i



J =

∂x ∂x ˆ ∂x ∂ yˆ

∂y ∂x ˆ ∂y ∂ yˆ

 

must be nonsingular for Fe to be invertible

Isoparametric mappings: it is possible to define curved elements e using a mapping Fe of the same degree as the basis functions on the reference element eˆ Example. Extended quadratic element P2+ e = Fe (ˆ e) x = Fe (ˆ x) =

ˆ i → xi , x 3 P

i=1

ˆ ij → xij , x

xi ϕ ˆi (ˆ x) +

P ij

ˆ 123 → x123 x

xij ϕ ˆij (ˆ x) + x123 ϕ ˆ123 (ˆ x)

y^

Fe

1

e

e^

1 2

0

y

1 2

1

x^

0

x

Construction of quadrilateral finite elements Idea: construct 2D basis functions as a tensor product of 1D ones defined on eˆ 1. Bilinear finite elements Let

λ1 (t) = 1 − t,

y^

λ2 (t) = t,

t ∈ [0, 1]

ϕˆ1 (ˆ x) = λ1 (ˆ x)λ1 (ˆ y ),

ϕˆ3 (ˆ x) = λ2 (ˆ x)λ2 (ˆ y)

ϕˆ2 (ˆ x) = λ2 (ˆ x)λ1 (ˆ y ),

ϕˆ4 (ˆ x) = λ1 (ˆ x)λ2 (ˆ y)

1

y x^ 4 x^ 1

x^ 3

e^

0

x^ 2 1

e

Fe x1

x^

x3

x4

0

x2

x

The space Q1 (ˆ e) spanned by ϕˆi consists of functions which are P1 for each variable In general

Qk (ˆ e) = span{xk11 xk22 . . . xknn },

Isoparametric mapping

x = Fe (ˆ x) =

4 P

i=1

0 ≤ ki ≤ k,

xi ϕˆi (ˆ x),

i = 1, . . . , n

ϕi (x) = ϕˆi (Fe−1 (x))

The physical element e = Fe (ˆ e) is a quadrilateral with straight sides which must be convex for Fe to be invertible. It is easy to verify that Fe (ˆ xi ) = xi , i = 1, . . . , 4

Construction of quadrilateral finite elements 2. Nonconforming rotated bilinear elements Let

˜ 1 (ˆ Q e) = span{1, x ˆ, yˆ, x ˆ2 − yˆ2 } 2

yˆ

x = Fe (ˆ x) =

αi xi

ˆ4 x

bilinear mapping 1 |Si |

R

Edge-oriented basis functions:

Si

0

ˆ1 x

uh (x) =

ϕˆ = [ϕˆ1 , ϕˆ2 , ϕˆ3 , ϕˆ4 ]T ,

c = [c1 , c2 , c3 , c4 ]T ,

ψˆ = [1, x ˆ, yˆ, x ˆ2 − yˆ2 ]T , ⇒

ψˆT c = ψˆT A−1 u = ϕˆT u

x4

x2

e x1

xˆ

0

edge mean values

uj ϕˆj (ˆ x) =

u = [u1 , u2 , u3 , u4 ]T ,

Ac = u

1

uh (x(s)) ds ≈ uh (xi ) 4 P

Fe

ˆ2 x

eˆ

j=1

Coefficients:

x3

1

i=1

Degrees of freedom: ui =

y

ˆ3 x

2

uh (x) = c1 + c2 x ˆ + c3 yˆ + c4 (ˆ x − yˆ ) 4 P

(Rannacher and Turek, 1992)

4 P

j=1

ui =

cj ψˆj (ˆ x), ∀x ∈ e P

aij cj R x)ds aij = |S1i | Si ψˆj (ˆ ⇒

j

ϕˆT = ψˆT A−1

x

Construction of quadrilateral finite elements Midpoint rule:

aij ≈ ψˆj (ˆ xi ),

ui ≈ uh (xi )

˜a Q 1

{edge mean values}   1 1 1 2 0 3   2   1 1 1  2 3  A=   1 1 1 −2  2 3   1 0 12 − 31

ϕˆ1 (ˆ x) =

3 4

+ 23 x ˆ − 52 yˆ − 32 (ˆ x2 − yˆ2 )

ϕˆ2 (ˆ x) = − 41 − 12 x ˆ + 32 yˆ + 23 (ˆ x2 − yˆ2 )

ˆ − 12 yˆ − 23 (ˆ x2 − yˆ2 ) ϕˆ3 (ˆ x) = − 41 + 32 x ϕˆ4 (ˆ x) =

3 4

− 25 x ˆ + 32 yˆ + 32 (ˆ x2 − yˆ2 )

exact for linear functions

˜b Q 1

{edge midpoint values}   1 1 1 2 0 4   3   1 1 1  2 4  A=   1 1 1 −3  2 4   1 0 12 − 41

ϕˆ1 (ˆ x) =

3 4

+x ˆ − 2ˆ y − (ˆ x2 − yˆ2 )

ϕˆ2 (ˆ x) = − 41 + yˆ + (ˆ x2 − yˆ2 )

ˆ − (ˆ x2 − yˆ2 ) ϕˆ3 (ˆ x) = − 41 + x

ϕˆ4 (ˆ x) =

3 4

− 2ˆ x + yˆ + (ˆ x2 − yˆ2 )

Nonparametric version: construct the basis functions directly using a local coordinate system rather than the transformation to a reference element

Construction of quadrilateral finite elements 3. Biquadratic finite elements θ1 (t) = (1 − t)(1 − 2t),

θ2 (t) = t(1 − 2t),

θ3 (t) = 4t(1 − t),

t ∈ [0, 1]

Products of 1D quadratic basis functions spanning the space Q2 (ˆ e) ϕˆ1 (ˆ x) = θ1 (ˆ x)θ1 (ˆ y ),

ϕˆ4 (ˆ x) = θ1 (ˆ x)θ2 (ˆ y ),

ϕˆ7 (ˆ x) = θ3 (ˆ x)θ2 (ˆ y)

ϕˆ2 (ˆ x) = θ2 (ˆ x)θ1 (ˆ y ),

ϕˆ5 (ˆ x) = θ3 (ˆ x)θ1 (ˆ y ),

ϕˆ8 (ˆ x) = θ1 (ˆ x)θ3 (ˆ y)

ϕˆ3 (ˆ x) = θ2 (ˆ x)θ2 (ˆ y ),

ϕˆ6 (ˆ x) = θ2 (ˆ x)θ3 (ˆ y ),

ϕˆ9 (ˆ x) = θ3 (ˆ x)θ3 (ˆ y)

ϕi (x) = ϕˆi (Fe−1 (x)),

Basis functions on the physical element

∀x ∈ e

Mapping: subparametric (bilinear) or isoparametric e

y

x7

x4 x8 x1 0

x3 x9

x6

Fe1

y^ 1

x5 x2

x

0

e^

Fe2

x^ 4 x^ 7 x^ 3 x^ 8 x^ x^ 6 9 x^ 1 x^ 5 x^ 2 1

x^

e

y

x7 x3 x6 x8 x9 x2 x1 x5

x4

0

x

x = Fe (ˆ x) =

9 P

xi ϕˆi (ˆ x)

i=1

e = Fe (ˆ e) is curved

Construction of hexahedral finite elements 1. Trilinear finite elements

λ1 (t) = 1 − t,

λ2 (t) = t,

t ∈ [0, 1]

Products of 1D linear basis functions spanning the space Q1 (ˆ e) ϕˆ1 (ˆ x) = λ1 (ˆ x)λ1 (ˆ y )λ1 (ˆ z ),

ϕˆ5 (ˆ x) = λ1 (ˆ x)λ1 (ˆ y )λ2 (ˆ z)

ϕˆ2 (ˆ x) = λ2 (ˆ x)λ1 (ˆ y )λ1 (ˆ z ),

ϕˆ6 (ˆ x) = λ2 (ˆ x)λ1 (ˆ y )λ2 (ˆ z)

ϕˆ3 (ˆ x) = λ2 (ˆ x)λ2 (ˆ y )λ1 (ˆ z ),

ϕˆ7 (ˆ x) = λ2 (ˆ x)λ2 (ˆ y )λ2 (ˆ z)

ϕˆ4 (ˆ x) = λ1 (ˆ x)λ2 (ˆ y )λ1 (ˆ z ),

ϕˆ8 (ˆ x) = λ1 (ˆ x)λ2 (ˆ y )λ2 (ˆ z)

Basis functions on the physical element z^

x^ 4

1

x^ 8 x^ 5

y^

x^ 1 0 1

e^

z

x^ 3

Fe

x^ 7 x^ 6

x5

x^ 2 1

x8

x^

0

x4 x1

e

ϕi (x) = ϕˆi (Fe−1 (x)),

∀x ∈ e

x3 x7

x2

Isoparametric mapping

x6

x

x = Fe (ˆ x) =

8 P

xi ϕˆi (ˆ x)

i=1

y

2. Rotated trilinear elements (6 nodes, face-oriented degrees of freedom)

Finite element matrix assembly Example: 1D Poisson equation   − d2 u = f in (0, 1) dx2  u(0) = 0, du (1) = 0

'k

ek

dx

Galerkin discretization:

'k

1

xk

uh =

N P

u j ϕj

1

xk

(linear finite elements)

j=1 N X

u0 = 0,

uj

Z

1

0

j=1

dϕi dϕj dx = dx dx

Z

1

f ϕi dx,

0

Decomposition of integrals into element contributions ak ij

Au = F,

N X j=1

uj

z N Z X

k=1

|

ek

ek = [xk−1 , xk ]

Fik

}| { z }| { N Z X dϕi dϕj f ϕi dx, dx = dx dx k=1 ek {z } | {z } aij

∀i = 1, . . . , N

Fi

∀i = 1, . . . , N

Example: 1D Poisson equation / linear elements Idea: evaluate element contributions and insert them into the global matrix aij = Fi =

N P

k=1 N P

k=1

akij =

Fik =

R1

dϕi dϕj 0 dx dx

R1 0

f ϕi dx

dx

akij 6= 0 only for i, j ∈ {k − 1, k} Fik 6= 0 only for i ∈ {k − 1, k}

Element stiffness matrix and load vector ek = [xk−1 , xk ]  R  R   R dϕk−1 dϕk dϕk−1 dϕk−1 f ϕk−1 dx dx dx dx dx dx dx  ek ek k k ek    a = , f = R R dϕk dϕk R dϕk dϕk−1 f ϕk dx dx ek dx dx ek dx dx ek dx Coefficients of the global system Au = F  · · ·  · ak−1 k−1 ak−1 k A=  · ak k−1 ak k · · ·

· · · ·

which are to be augmented    ·  Fk−1   ,  F =  Fk   ·

Example: 1D Poisson equation / linear elements Special case: ϕk−1 (x) = ϕk (x) =

3 elements, xk −x xk −xk−1

x−xk−1 xk −xk−1

=

=

∆x = 31 ,

k∆x−x ∆x ,

x−(k−1)∆x , ∆x

∀x ∈ ek = [xk−1 , xk ]

Hence,

dϕk−1 dx dϕk dx

f ≡1

'0

'1

1 = − ∆x

=

ak =

0

1 ∆x 1 ∆x



e1 x0

1 −1 −1 1

e2 x1



,

'3

'2

fk =

e3

1

x3

x2

∆x 2



1 1



Assembly of the global stiffness matrix and load vector       1 −1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 −1 0  1  −1 1  0 0 1  0 A = ∆x  0 0 0 0  + ∆x  0 −1 1 0  + ∆x  0 0 1 −1  ⇒ 0 0 −1 1 0 0 0 0 0 0 0 0           1 −1 0 0 1 0 0 1 2 −1 0   1  ∆x  1  ∆x  0  ∆x  2  1  −1 A = ∆x  0 −1 2 −1  , F = ∆x 2  0 + 2  1 + 2  1 = 2  2  0 0 −1 1 1 0 0 1

Example: 1D Poisson equation / linear elements Recall that u0 = 0 so the first equation drops out and the system shrinks to      1 2 −1 0 u 1 u1 = 25 (∆x)2 , u2 = 4(∆x)2 1  −1 2 −1   u2  =  1  ⇒ 2 (∆x) 1 u3 = 92 (∆x)2 , ∆x = 31 0 −1 1 u3 2

Implementation of Dirichlet boundary conditions 1. Row/column elimination: u0 = g0 whereas the second one turns into

⇒ the first equation is superfluous a11 u1 + a12 u2 + a13 u3 = F1 − a10 g0

2. Row modification (replacement by a row of the identity matrix) a00 := 1,

a0j := 0,

j = 1, 2, 3,

F0 := g0

3. Penalty method / addition of a large number α to the diagonal a00 := a00 + α,

F0 := F0 + αg0

symmetry is preserved

Implementation of Neumann boundary conditions R du ⇒ FN = eN f ϕN dx + g1 a surface integral is added dx (1) = g1

Example: 1D Poisson equation / quadratic elements N P

Galerkin FEM:

uj

k=1

j=1/2

ek = [xk−1 , xk ],

ek

dϕi dϕj dx dx

dx =

N R P

k=1

x = {λ1 (x), λ2 (x)}

ϕk−1 = λ1 (2λ1 − 1),

ϕk = λ2 (2λ2 − 1), ϕk−1/2 = 4λ1 λ2 ,

N R P

ek

'k

dϕk−1 4λ1 −1 = − dx ∆x dϕk−1/2 2 −1 = 4λ∆x dx dϕk λ1 −λ2 dx = 4 ∆x

∀i = 1, . . . , N

f ϕi dx, 1

1=2

'k

1=2

xk

'k

ek xk

1

xk

Element stiffness matrix and load vector

ak =

dϕk−1 dϕk−1 dx dx dx ek 6 R dϕk−1/2 dϕk−1 6 dx dx dx 4 ek R dϕk dϕk−1 dx dx ek dx

2 R

dϕk−1 dϕk−1/2 dx dx dx ek R dϕk−1/2 dϕk−1/2 dx dx dx ek R dϕk dϕk−1/2 dx dx ek dx

 R  f ϕ dx k−1 R ek k f =  ekR f ϕk−1/2 dx  = f ϕk dx ek

R

∆x 6





1 4 1

3 dϕk−1 dϕk dx dx dx ek 7 R dϕk−1/2 dϕk dx 7 dx dx 5 ek R dϕk dϕk dx ek dx dx R

=

1 3∆x

2

3 7 −8 1 4 −8 16 −8 5 1 −8 7

Global system: Au = F, where u = [u1/2 u1 u3/2 . . . uN −1/2 uN ]T

Numerical integration for finite elements R

Change of variables theorem

f (x) dx =

e

ϕi (x) = ϕˆi (Fe−1 (x)),

R

fˆ(ˆ x)| det J| dˆ x

eˆ

ˆ ϕˆi = J∇ϕi ∇

∀x ∈ e

ˆ ϕˆi ∇ϕi = J −1 ∇

⇒

For instance, the entries of the element stiffness matrix are given by Z Z ˆ ϕˆj )| det J| dˆ ˆ ϕˆi ) · (J −1 ∇ x aij = ∇ϕi · ∇ϕj dx = (J −1 ∇ e

Numerical integration

eˆ

R

eˆ

gˆ(ˆ x) dˆ x≈

n P

w ˆi gˆ(ˆ xi ),

gˆ(ˆ x) = fˆ(ˆ x)| det J|

i=0

Newton-Cotes formulae can be used but Gaussian quadrature is preferable: R

gˆ(ˆ x) dˆ x≈ eˆ

g ˆ(ˆ x1 )+ˆ g (ˆ x2 ) , 2

x ˆ1 =

1 2

−

1 6

√

3,

x ˆ2 =

1 2

+

1 6

√

3,

eˆ = [0, 1]

exact for gˆ ∈ P3 (ˆ e) as compared to P1 (ˆ e) for the trapezoidal rule

Storage of sparse matrices Banded matrices: store the nonzero diagonals as 1D arrays Arbitrary matrices: store the nonzero elements as a 1D array 1. Coordinate storage

(inconvenient access)

A(NNZ)

nonzero elements in arbitrary order

IROW(NNZ)

auxiliary array of row numbers

ICOL(NNZ)

auxiliary array of column numbers

2. Compact storage

(convenient access)

A(NNZ)

nonzero elements stored row-by-row

ILD(N+1)

pointers to the beginning of each row

ICOL(NNZ)

auxiliary array of column numbers

Example



1 2 A= 0 7

2 4 3 0

0 3 6 5

 7 0 5 8

A = (1, 2, 7, 4, 2, 3, 6, 3, 5, 8, 7, 5) ICOL = (1, 2, 4, 2, 1, 3, 3, 2, 4, 4, 1, 3) IROW = (1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4) NNZ = 12, ILD = (1, 4, 7, 10, 13)