Discontinuous Galerkin Methods: An Overview and Some Applications - - PowerPoint PPT Presentation

Discontinuous Galerkin Methods: An Overview and Some Applications Daya Reddy

UNIVERSITY OF CAPE TOWN

Joint work with Beverley Grieshaber and Andrew McBride

SANUM, Stellenbosch, 22 – 24 March 2016

§ A model elliptic problem: weak or variational formulations § The Galerkin finite element method: analysis and approximations § Discontinuous Galerkin (DG) formulations § Near-incompressibility in elasticity

Structure of talk

Minimization of an “energy” The solution satisfies the weak problem for all functions which satisfy on the boundary Sufficiently smooth satisfies the Poisson equation and boundary condition

min

v

J(v) J(v) = 1

2

Z

Ω

|rv|2 dx Z

Ω

fv dx Z

Ω

ru · rv dx = Z

Ω

fv dx −∆u = f Ω u = Γ

• n ¡
• n ¡

Ω Γ

f

u v v = 0

u

Model problem: deformation of a membrane

∆u = ∂2u ∂x2 + ∂2u ∂y2

Define the bilinear form and linear functional Thus the above problem is

• r equivalently

The membrane problem, continued

a(·, ·) `(·) a : V ⇥ V ! R , a(u, v) = Z

Ω

ru · rv dx ` : V → R , `(v) = Z

Ω

fv dx Z

Ω

ru · rv dx = Z

Ω

fv dx min

v∈V J(v) := 1 2

Z

Ω

|rv|2 dx Z

Ω

fv dx min

v∈V 1 2a(v, v) − `(v)

a(u, v) = `(v) ∀v ∈ V

Interlude: the Sobolev spaces

Built from the Lebesgue space of square-integrable functions: Define, for integer ,

Seminorm for problem in

Hilbert space with induced norm We will also need

L2(Ω) = ⇢ v : Z

Ω

v2 dx := kvk2

0 < 1

• m ≥ 0

Hm(Ω) =

• v : Dαv ∈ L2(Ω) , |α| ≤ m

Dαv =

∂|α|v ∂xα1

1 ∂xα2 2

|α| = α1 + α2

R2

kvk2

m =

X

|α|≤m

|v|2

m

H1

0(Ω) = {v ∈ H1(Ω) : v = 0 on Γ}

|v|2

m =

X

α=m

Z

Ω

|Dαv|2 dx

Hm(Ω)

e.g. kvk2

1 =

Z

Ω

h |v|2 + ⇣ ∂v ∂x1 ⌘2 + ⇣ ∂v ∂x2 ⌘2i dx

Well-posedness of the variational problem

This problem has a unique solution if: § W is a closed subspace of a Hilbert space H § a is coercive or W-elliptic: § a is continuous: § is continuous: ¡ a(w, w) αkwk2

H

` |a(w, z)|  MkwkHkzkH |`(z)|  CkzkH min

w∈W 1 2a(w, w) − `(w)

The model problem has a unique solution in H1

0(Ω)

• 1. Partition the domain into subdomains or finite elements
• 2. Construct a basis for comprising continuous functions that are

polynomials on each element

Finite element approximations

{ϕi}N

i=1

Aim: to pose the variational problem on a finite-dimensional subspace V h ⊂ V V h

portion of hip replacement: physical object and finite element model

• 3. The piecewise-polynomial approximations can

be written

• 4. Substitute in the weak formulation to obtain

ϕi uh = X

i

ϕi(x)ui ≡ ϕu vh = X

i

ϕi(x)vi ≡ ϕv X

i

a('i, 'j) | {z }

Kji

ui = `('j) | {z }

Fj

Ku = F

X

i,j

vi[a('i, 'j)ui − `('j)] = 0 a(uh, vh) = `(vh)

The Galerkin finite element method

Construct and seek such that for all mesh size Define the error by : under what conditions do we have convergence in the sense that Orthogonality of the error: a(uh, vh) = `(vh) for all vh ∈ Vh lim

h→0 uh = u?

Vh ⊂ V vh ∈ Vh uh ∈ Vh u − uh a(u − uh, vh) = a(u, vh) − a(uh, vh) = `(vh) − `(vh) =

Ku = F

hT = diameter of T h = max

T∈T hT

hT

Convergence of finite element approximations

ku uhk1  ¯ C inf

vh∈Vh ku vhk1

(V-­‑ellipticity) αku uhk2

V

 a(u uh, u uh) = a(u uh, u vh) + a(u uh, vh uh)  a(u uh, u vh)  Mku uhkV ku vhkV (orthogonality of error) (continuity) Céa’s lemma Strategy for obtaining error bound: a) choose to be the interpolate of in b) use interpolation error estimate to bound actual error vh u Vh

An a priori estimate

Ciarlet and Raviart Arch. Rat. Mech. Anal. 1972

Let be a triangulation of a bounded domain with polygonal boundary: Define the mesh size

A family of triangulations is regular as if there exists such that hT = diameter of T ρT = sup{diameter of B; B a ball contained in T} T Ω ¯ Ω = ∪T∈T T h = max

T∈T hT

h → 0 σT ≤ σ for all T ∈ Th σ > 0 σT = hT /ρT

hT ρT

Finite element interpolation theory

(Global estimate) Let be a uniformly regular triangulation of a polygonal

• domain. Define

(global interpolator) T

Vh = {vh ∈ C(¯ Ω) : vh|T ∈ Pk} ∩ V

h = max

T

hT (Local estimate) For a regular triangulation with and the interpolation operator which maps functions to polynomials of degree , v ∈ Hk+1(T) , k + 1 ≥ m |v − πv|m,T ≤ Chk+1−m|v|k+1,T Πh : H2(Ω) → Vh kv Πhvkm,Ω  Chk+1−m|v|k+1,Ω m = 0, 1

v Πhv hT

π ≤ k

Finite element interpolation theory

ku uhkV  C infvh∈Vh ku vhkV  Cku ΠhukV  Chmin(k,r−1)|u|r ¡ ¡ ¡ ¡ ¡ So for the simplest approximation, by piecewise-linear simplices, for the second-order elliptic equations Could use piecewise-quadratic simplices in which case

• - provided that the solution is smooth enough to belong to !

¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ku uhk1  Ch|u|2 (k = 2) ku uhk1 = O(h2) H3(Ω)

Convergence of finite element approximations

Discontinuous Galerkin (DG) methods

• Drop the condition of continuity
• So
• A member

is now a polynomial on each element, but not continuous across element boundaries

• An immediate consequence: much greater number of degrees of freedom

Vh 62 V

K1 K2 u e12 { {u} } N1 N2

vh ∈ Vh

Discontinuous Galerkin (DG) methods

0.02 0.04 0.06 0.08 0.1 0.12 0.14

(a) Conforming approximation

0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16

(b) Discontinuous Galerkin approximation

(c) Discontinuous Galerkin approximation with a low penalty parameter

What do solutions look like?

Can accommodate hanging nodes Useful in adaptive mesh refinement

Why bother with DG?

0.1 1 10 1 10 100 1000 10000 Displacement Error (m) Total DOF Conforming DG

Figure 6. Plot of the L2-norm over B0 of the error in the displacement field as a function of the total number of degrees of freedom for the test in Section 6.1. Curves are shown for the discontinuous Galerkin method as well as for the conforming one. Remarkably, the two curves

• verlap, indicating that both methods provide the same accuracy for the same computational

cost, i.e. for this example they are equally efficient.

ten ¡Eyck ¡and ¡Lew ¡2010 ¡ ¡ (nonlinear ¡elas7city) ¡

Efficiency and accuracy

• For the problem of deformations of elastic bodies DG methods show

good behaviour for near-incompressibility with the use of low-order triangles in two dimensions, and tetrahedra in three

• DG with quadrilaterals and hexahedra less straightforward:

– poor behaviour, including locking, for low-order approximations

• We show why this is so, and propose some remedies

Our objective

First, review derivation of heat equation

heat flux vector ¡

Heat equation: Balance of energy + Fourier’s law give the (steady) heat equation

div q = s q = krϑ −k ∆ϑ = s

heat source ¡

q = s =

Elasticity: Equilibrium + Hooke’s law give Navier’s equation

σ = λ div u + (ru + [ru]T ) −div σ = f ¯ ∆u = [∆u + (1 + λ)rdiv u] = f

stress tensor or matrix σ ¡ displacement vector u

Governing equations for elasticity

Γ Ω Γ u(x) x

Equivalent to minimization problem Define Then the minimization problem is

• r equivalently

which has a unique solution in Furthermore (Brenner and Sung 1992) the solution satisfies min u J(u) = 1

2

Z

Ω

|rsu|2 + (1 + λ)(div u)2 dV Z

Ω

f · u dV

`(v) = Z

Ω

f · v dx

a(u, v) = Z

Ω

[rsu : rsv + (1 + λ)div u div v] dx min

v∈V 1 2a(v, v) − `(v)

rsu = 1

2[ru + (ru)T ]

(rsu)ij = 1 2 ✓ ∂ui ∂xj + ∂uj ∂xi ◆

a(u, v) = `(v) V := [H1

0(Ω)]d

kukH2 + λk div ukH1  CkfkL2 u ∈ [H2(Ω)]2

The use of low-order elements leads to a non-physical solution in the incompressible limit We explore the use of DG methods as a remedy for locking

Q1

• ­‑0,2 ¡
• ­‑0,1 ¡

0 ¡ 0,1 ¡ 0,2 ¡ 0,3 ¡ 0,4 ¡ 0,5 ¡ 0,6 ¡ 0,7 ¡ 0,8 ¡

• ­‑0,6 ¡
• ­‑0,4 ¡
• ­‑0,2 ¡

0 ¡ 0,2 ¡ 0,4 ¡ 0,6 ¡ Case ¡III ¡ Case ¡V ¡(EAS) ¡ MHW ¡Case ¡1 ¡ HW ¡Case ¡I ¡or ¡Q1 ¡

Q1

Locking

a(u, v) = Z

Ω

[rsu : rsv + (1 + λ)div u div v] dx λ → ∞ div u → 0

Jumps and averages:

• n an interior edge of element ,

For example On an edge that forms part of the boundary,

T + T −

n+ n−

E

[v] = v+ · n+ + v− · n−, [[v]] = v+ ⊗ n+ + v− ⊗ n−, {v} = 1

2(v+ + v−)

[[τ]] = τ +n+ + τ −n−, {τ} = 1

2(τ + + τ −)

[[v]] = v ⊗ n {τ} = τn

[[τ]] = (τ + − τ −)n+

T

E

DG formulation: some definitions

Z

Ω

div σ · v dV = X

T

Z

T

div σ · v dV

• take the inner product with a test function
• integrate, and integrate by parts

Consider the left-hand side:

− Z

Ω

div σ · v dV = Z

V

f · v dV

Setting it up

¯ ∆u = [∆u + (1 + λ)rdiv u] | {z }

−div σ

= f σ = λ div u + (ru + [ru]T ) Z

Ω

div σ · v dV = X

T

Z

∂T

σn · v ds X

T

Z

T

σ : rsv dV

X

T

Z

∂T

v · σn ds = X

E

Z

E

[[v]] : {σ} ds + X

Eint

Z

Eint

{v} · [[σ]] ds

Thus weak equation is now Next, we need the “magic formula”

T

E ∂T

• X

T

Z

∂T

σn · v ds + X

T

Z

T

σ : rsv dV = Z

Ω

f · v dV

This gives Since the exact solution is smooth ( ) and the stress is continuous we can assume that which leaves ¡

u ∈ [H2(Ω)]d

• X

E

Z

E

[[v]] : {σ} ds X

Eint

Z

Eint

{v} · [[σ]] ds + Z

Ω

σ : rsv dV = Z

Ω

f · v dV [[σ]] = 0

• X

E

Z

E

[[v]] : {σ} ds + Z

Ω

σ : rsv dV = Z

Ω

f · v dV

Exploiting again the smoothness of the exact solution we can add a term to symmetrize the problem: Finally, we may need to stabilize the problem: this gives us the symmetric interior penalty (SIP) formulation (Douglas and Dupont 1976) ¡

[[u]] = 0

• X

E

Z

E

[[v]] : {σ(u)} ds X

E

Z

E

[[u]] : {σ(v)} ds + X

E

Z

Ω

σ : rsv dV = Z

Ω

f · v dV

• X

E

Z

E

[[v]] : {σ(u)} ds X

E

Z

E

[[u]] : {σ(v)} ds + k X

E

Z

E

1 hE [[u]][[v]] + X

E

Z

Ω

σ : rsv dV = Z

Ω

f · v dV

The DG formulation

P1(T) : vh = a0 + a1x + a2y Vh =

• vh : vh ∈ [L2(Ω)]d , vh|T ∈ P1(T) or Q1(T)

Q1(T) : vh = a0 + a1x + a2y + a3xy

ah(u, v)

DOUGLAS & DUPONT 1976, ARNOLD 1982, HANSBO & LARSON 2002 RIVIÈRE & WHEELER 1999, 2000; ODEN, BABUSKA & BAUMANN 1998 DAWSON, SUN AND WHEELER 2004

Nonsymmetric Interior Penalty Galerkin (NIPG)

θ =              +1 −1

Symmetric Interior Penalty Galerkin (SIPG) Incomplete Interior Penalty Galerkin (IIPG) X

T

Z

T

σ(u) : rsv dV + θ X

E

Z

E

[[u]] : {σ(v)} ds X

E

Z

E

[[v]] : {σ(u)} ds +k X

E

Z

E

1 hE [[u]][[v]] = Z

Ω

f · v dV Vh =

• vh : vh ∈ [L2(Ω)]d , vh|T ∈ P1(T) or Q1(T)

ah(uh, vh) = `(vh)

WIHLER 2002

• x

y O

10

10

10

10

10

10

10

number of degrees of freedom absolute L2 error DGFEM on graded mesh λ=1 λ=100 λ=500 λ=1000 λ=5000 0.99 1

−1 −0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1

DG with triangles is uniformly convergent

• Fig. 8. Tensile stress contours ðry > 0:15 MPaÞ of DG solutions of beam problem with Poisson’s ratio 0.499. (a) OBB; (b) SIPG; dp ¼ 30; (c) NIPG; dp ¼ 10; and (d) IIPG; dp ¼ 50.

LIU, WHEELER AND DAWSON 2009 These authors report good results using all methods OBB ¡ IIPG ¡ NIPG ¡ SIPG ¡

Some computational results with quads appear to show locking-free behaviour …

ALL ¡ ¡ ALL ¡ ¡ Grieshaber, McBride and R. (2015) ¡

… but not so in other cases

• We know that DG works for simplicial elements
• Some authors report good results for quads
• We find locking for a simple example
• What’s going on?

(I) ¡ (II) ¡ (III) ¡ (IV) ¡

(note the use of two stabilization terms)

Let’s see why DG with triangles works

ah(u, v) := X

T ∈Th

Z

T

σ(u) : rsv + θ X

E∈Γ

Z

E

[[u]] : {σ(v)} X

E∈Γ

Z

E

{σ(u)} : [[v]] +kµµ X

E∈Γ

Z

E

1 hE [[u]] : [[v]] + kλλ X

E∈Γ

1 hE Z

E

[[u]] : [[v]] := 2µ " X

T ∈Th

Z

T

rsu : rsv dx + θ X

E∈Γ

Z

E

[[u]] : {rsv} X

E∈Γ

Z

E

{ru} : [[v]] + kµ X

E∈Γ

Z

E

1 2hE [[u]] : [[v]] # +λ " + X

T ∈Th

Z

T

(div u)(div v) + θ X

E∈Γ

Z

E

[[u]] : {(div v)I} X

E∈Γ

Z

E

{(div u)I} : [[v]] ds +kλ X

E∈Γ

1 hE Z

E

[[u]] : [[v]] #

Error analysis for linear triangles: NIPG

Use Crouzeix-Raviart interpolant: linear on triangles, continuous at midpoints of edges Properties:

interpolation error Approximation error

e = u − uh Z

E

(u − Πu) · n ds = 0 Z

T

div(u − Πu) dV = 0 = u − Πu | {z } η + Πu − uh | {z } ξ

Πu = interpolant of u

λ-terms in error bound can be handled as follows: for example, to give eventually Use estimate to get uniform convergence ¡

¡

e = η + ξ η = u − Πu ξ = Πu − uh ξ ∈ Vh ⇒ div ξ, σ(ξ) const.

kek2

DG

 Ch2 ⇣ kuk2

H2(Ω) + λ2k div uk2 H1(Ω)

⌘

kukH2 + λk div ukH1  CkfkL2

(I) ¡= ¡ λ

X

T

Z

T

(div η)(div ξ) ds = λ X

T

div ξ Z

T

div η ds = 0

Convergence analysis for quadrilaterals

First need to construct a suitable interpolant As before and must satisfy

e = u − uh Z

E

(u − Πu) · n ds = 0 Z

E

[[u − Πu]] · n ds = 0 Z

∂T

(u − Πu) · n ds = 0 Πu ∈ Vh Πu = u − Πu | {z } η + Πu − uh | {z } ξ

The interpolant

Inspired by DOUGLAS, SANTOS, SHEEN & YE 1999; CIA, DOUGLAS, YE 1999: Construct interpolant with span

⇢ 1, x, y, θ(x) + κ(y) 1, x, y, κ(x) + θ(y)

• θ(x) = 3x2 − 10x4 + 7x6

κ(x) = −3x2 + 5x4

Define as orthogonal projection on element ¡

Πu

Error bound

Eventually get so not possible to bound λ-dependent term k˜ u uhk2

DG  C

X

T

h2

T

⇣ kuk2

H2(T ) + λ2kuk2 H2(T ) + λ2k div uk2 H1(T )

⌘

Example: square plate

u1(x, y) = sin 2πy

• −1 + cos 2πx
• +

1 1 + λ sin πx sin πy u2(x, y) = sin 2πx

• 1 − cos 2πy
• +

1 1 + λ sin πx sin πy

BRENNER ¡SINUM ¡1993 ¡ Ω = (0, 1)2 µ = 1

NIPG ¡ ¡ IIPG ¡ ¡ ALL ¡ ¡ SIPG ¡ ¡

Recall classical SRI for elasticity problem: to overcome volumetric locking one replaces by

• r

A remedy: selective reduced integration (SRI)

Π0 div v = 1 |T| Z

T

div v dV

Z

Ω

⇥ rsuh : rsvh + λ(div u)(div v) | {z }

volumetric term

⇤ dV Z

Ω

⇥ rsuh : rsvh dV + λ X

T

wT div u(xT ) div v(xT ) | {z }

underintegration of volumetric term

Z

Ω

⇥ rsuh : rsvh dV + λ X

T

Z

T

(Π0 div u)(Π0 div v) dV

Underintegration applied to (II) gives Replace (IV) with But we now have to check coercivity and consistency of the modified bilinear form! ¡ θλ Z

E

[[ξ]]{div η} ds ' θλ Z

E

Π0[[ξ]] Π0{div η} ds = θλ Π0{div η} Z

E

Π0[[ξ]] ds = θλ Π0{div η} Z

E

[[ξ]] ds = 0 kλλθ X

E

Z

E

Π0[[ξ]] Π0[[η]] = kλλθ X

E

Z

E

[[Π0ξ]] [[η]] = kλλθ X

E

Π0[[ξ]] Z

E

[[η]] =

Underintegration applied to (II) gives Eventually get But we now have to check coercivity and consistency of the modified bilinear form! ¡ θλ Z

E

[[ξ]]{div η} ds ' θλ Z

E

Π0[[ξ]] Π0{div η} ds = θλ Π0{div η} Z

E

Π0[[ξ]] ds = θλ Π0{div η} Z

E

[[ξ]] ds = 0 ah(η, ξ)  CkξkDG ⇣ P

T h2 T

⇣ kuk2

H2 + λ2kdiv uk2 H1

| {z }

Ckfk2

⌘⌘1/2

Coercivity

NIPG: OK for SRI on terms (II) and (III) SIPG: OK for SRI on terms (II), (III), (IV) IIPG: OK for SRI on terms (III) and (IV)

Term (III) when underintegrated leads to a consistency error which can be controlled

|ERI(u, ξ)| = λ

• X

E

Z

E

• Π0{div u} tr Π0[[ξ]] − {div u} tr[[ξ]]
• ds

Locking-free (uniformly convergent) behaviour

SIP ✔ ✖ ✔ NIP ✔ ✖ ✔ IIP ✔ ✖ ✔

Q1 P1

Q1 + SRI ✕ ¡ ✕ ¡ ✕ ¡ ✕ ¡

Square plate: DG with triangles quads with and without SRI

Standard finite elements / DG quads without SRI ¡ DG triangles / quads with SRI ¡

NIPG ¡ IIPG ¡ SIPG ¡ DG+SRI ¡

DG ¡ DG ¡ DG ¡

DG+SRI ¡ DG+SRI ¡

u1 = (cos 2πx − 1)(sin 2πy sin πz − sin πy sin 2πz) + 1 1 + λ sin πx sin πy sin πz, u2 = (cos 2πy − 1)(sin 2πz sin πx − sin πz sin 2πx) + 1 1 + λ sin πx sin πy sin πz, u3 = (cos 2πz − 1)(sin 2πx sin πy − sin πx sin 2πy) + 1 1 + λ sin πx sin πy sin πz.

Cube with prescribed body force

10

−1

10 10

−1

10 10

1

10

2

log(h) log(H1 error) Q1 full Q1 UI P1 full SG

1

NIPG ¡

10

−1

10 10

−1

10 10

1

10

2

log(h) log(H1 error) Q1 full Q1 UI P1 full P1 UI SG

1

IIPG ¡

10

−1

10 10

−1

10 10

1

10

2

log(H1 error) Q1 full Q1 UI P1 full P1 UI SG

1

SIPG ¡

Q1 ¡SRI ¡ ¡ P1 ¡full ¡ ¡ Q1 ¡full ¡ ¡ SG ¡ ¡ Q1 ¡full ¡ ¡ Q1 ¡SRI ¡ ¡ SG, ¡P1 ¡full ¡ ¡ P1 ¡SRI ¡ ¡ SG, ¡P1 ¡full, ¡Q1 ¡full ¡ ¡ Q1 ¡SRI ¡ ¡ P1 ¡SRI ¡ ¡

Back to the T-bar example

Without SRI With SRI ¡

• LIU. WHEELER, DAWSON 2009 appeared to show good behaviour for IIPG

3D problem but effectively plane stress (Neumann or flux-free boundary condition) in out-of-plane direction, along ABCD Here treated as plane strain (Dirichlet or constrained displacement) along ABCD A ¡ B ¡ C ¡ C ¡ ¡

Concluding remarks

• For near-incompressibility DG with quadrilateral elements is not

straightforward

• A remedy, viz. selective under-integration λ-dependent edge terms, has been

proposed, analysed, and shown to converge uniformly at the optimal rate

• Arbitrary quadrilaterals: numerical experiments indicate behaviour similar to

that for rectangles. Analysis would be quite complex

• Current work: nonlinear problems; other parameter-dependent problems

References

Grieshaber BJ, McBride AT and Reddy BD, Uniformly convergent interior penalty approximations using multilinear approximations for problems in elasticity. SIAM J. Numer. Anal. 53 (2015) 2255–

• 2278. ¡

Grieshaber BJ, McBride AT and Reddy BD, Computational aspects of Discontinuous Galerkin approximations using quadrilateral and hexahedral elements. In review.

Alternative approach: P1 approximation on quads

Much of the manipulations carry over, all the bounds are as before, but this time Problematic terms: given properties of interpolant But so remains a problem for SIPG and IIPG, but is absent for NIPG

kλλ X

E

Z

E

[[γ]] [[w]] 6=

(II) ¡= ¡ (IV) ¡= ¡

uh ∈ P1 ξ = Πu − uh ∈ Vh ∼ P1 θλ Z

E

[[ξ]]{div η} ds = θλ X

E

]{div η} Z

E

[[ξ]] ds = kλλ 1 hE Z

E

[[ξ]] [[η]] ds 6= 0

Square plate, DG with quads and P1

NIPG ¡ IIPG ¡ SIPG ¡