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

Discontinuous Galerkin Methods: An Overview and Some Applications Daya Reddy U NIVERSITY OF C APE T OWN Joint work with Beverley Grieshaber and Andrew McBride SANUM, Stellenbosch, 22 – 24 March 2016

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

Model problem: deformation of a membrane f Minimization of an “energy” Z Z | r v | 2 dx � J ( v ) = 1 min J ( v ) fv dx 2 v Ω Ω The solution satisfies the weak problem u Ω Z Z Γ r u · r v dx = fv dx Ω Ω for all functions which satisfy on the boundary v = 0 v Sufficiently smooth satisfies the Poisson equation and boundary condition u on ¡ − ∆ u f Ω = ∆ u = ∂ 2 u ∂ x 2 + ∂ 2 u ∂ y 2 on ¡ u Γ = 0

The membrane problem, continued Z Z | r v | 2 dx � v ∈ V J ( v ) := 1 min fv dx 2 Ω Ω Z Z r u · r v dx = fv dx Ω Ω Define the bilinear form and linear functional a ( · , · ) ` ( · ) Z a : V ⇥ V ! R , a ( u, v ) = r u · r v dx Ω Z ` : V → R , ` ( v ) = fv dx Ω Thus the above problem is 1 min 2 a ( v, v ) − ` ( v ) v ∈ V or equivalently a ( u, v ) = ` ( v ) ∀ v ∈ V

Interlude: the Sobolev spaces H m ( Ω ) Built from the Lebesgue space of square-integrable functions: ⇢ � Z v 2 dx := k v k 2 L 2 ( Ω ) = v : 0 < 1 Ω Define, for integer , m ≥ 0 � v : D α v ∈ L 2 ( Ω ) , | α | ≤ m H m ( Ω ) = ∂ | α | v D α v = ∂ x α 1 1 ∂ x α 2 2 | α | = α 1 + α 2 Z | D α v | 2 dx Seminorm X | v | 2 m = R 2 for problem in Ω α = m k v k 2 X | v | 2 Hilbert space with induced norm m = m | α | ≤ m ⇣ ∂ v ⇣ ∂ v Z ⌘ 2 ⌘ 2 i h | v | 2 + e.g. k v k 2 dx 1 = + ∂ x 1 ∂ x 2 Ω We will also need H 1 0 ( Ω ) = { v ∈ H 1 ( Ω ) : v = 0 on Γ }

Well-posedness of the variational problem 1 min 2 a ( w, w ) − ` ( w ) w ∈ W This problem has a unique solution if: § W is a closed subspace of a Hilbert space H a ( w, w ) � α k w k 2 § a is coercive or W -elliptic: H § a is continuous: | a ( w, z ) |  M k w k H k z k H § is continuous: ¡ | ` ( z ) |  C k z k H ` The model problem has a unique solution in H 1 0 ( Ω )

Finite element approximations Aim: to pose the variational problem on a finite-dimensional subspace V h ⊂ V 1. Partition the domain into subdomains or finite elements { ϕ i } N 2. Construct a basis for comprising continuous functions that are V h i =1 polynomials on each element portion of hip replacement: physical object and finite element model

The Galerkin finite element method ϕ i 3. The piecewise-polynomial approximations can be written X u h = ϕ i ( x ) u i ≡ ϕ u i X v h = ϕ i ( x ) v i ≡ ϕ v i 4. Substitute in the weak formulation to obtain a ( u h , v h ) = ` ( v h ) X X v i [ a ( ' i , ' j ) u i − ` ( ' j )] = 0 a ( ' i , ' j ) u i = ` ( ' j ) Ku = F | {z } | {z } i,j i K ji F j

Convergence of finite element approximations Construct and seek such that for all V h ⊂ V u h ∈ V h v h ∈ V h a ( u h , v h ) = ` ( v h ) for all v h ∈ V h Ku = F mesh size h T = diameter of T h = max T ∈ T h T h T Define the error by : under what conditions do we have convergence in the u − u h sense that h → 0 u h = u ? lim Orthogonality of the error: a ( u − u h , v h ) = a ( u, v h ) − a ( u h , v h ) = ` ( v h ) − ` ( v h ) = 0

An a priori estimate α k u � u h k 2  a ( u � u h , u � u h ) ( V -­‑ ellipticity) V = a ( u � u h , u � v h ) + a ( u � u h , v h � u h )  a ( u � u h , u � v h ) (orthogonality of error)  M k u � u h k V k u � v h k V (continuity) k u � u h k 1  ¯ v h ∈ V h k u � v h k 1 C inf Céa’s lemma Strategy for obtaining error bound: a) choose to be the interpolate of in u V h v h b) use interpolation error estimate to bound actual error

Finite element interpolation theory Ciarlet and Raviart Arch. Rat. Mech. Anal. 1972 h T = diameter of T h T ρ T = sup { diameter of B ; B a ball contained in T } ρ T σ T = h T / ρ T Let be a triangulation of a bounded domain with polygonal boundary: T Ω ¯ Ω = ∪ T ∈ T T Define the mesh size h = max T ∈ T h T A family of triangulations is regular as if there exists such that h → 0 σ > 0 for all T ∈ T h σ T ≤ σ

Finite element interpolation theory (Local estimate) For a regular triangulation with v ∈ H k +1 ( T ) , k + 1 ≥ m and the interpolation operator which maps functions to π polynomials of degree , ≤ k | v − π v | m,T ≤ Ch k +1 − m | v | k +1 ,T (Global estimate) Let be a uniformly regular triangulation of a polygonal T domain. Define V h = { v h ∈ C (¯ Ω ) : v h | T ∈ P k } ∩ V Π h v Π h : H 2 ( Ω ) → V h (global interpolator) v h = max h T T h T k v � Π h v k m, Ω  Ch k +1 − m | v | k +1 , Ω m = 0 , 1

Convergence of finite element approximations k u � u h k V  C inf v h ∈ V h k u � v h k V  C k u � Π h u k V ¡ ¡ Ch min( k,r − 1) | u | r  ¡ ¡ ¡ So for the simplest approximation, by piecewise-linear simplices, k u � u h k 1  Ch | u | 2 for the second-order elliptic equations Could use piecewise-quadratic simplices in which case ( k = 2) k u � u h k 1 = O ( h 2 ) -- provided that the solution is smooth enough to belong to ! H 3 ( Ω ) ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡

Discontinuous Galerkin (DG) methods

Discontinuous Galerkin (DG) methods • Drop the condition of continuity • So V h 62 V • A member is now a polynomial on each element, but not continuous v h ∈ V h across element boundaries • An immediate consequence: much greater number of degrees of freedom � u � { { u } } N 1 K 1 K 2 N 2 e 12

What do solutions look like? 0.16 0.14 0.14 0.12 0.12 0.1 0.1 0.08 0.08 0.06 0.06 0.04 0.04 0.02 0.02 (a) Conforming approximation (b) Discontinuous Galerkin approximation 0.5 0.6 0.5 0.4 0.4 0.3 0.3 0.2 0.1 0 0.2 − 0.1 0.5 0.4 0.3 1 0.1 0.2 0.8 0.1 0.6 0.4 0 0.2 − 0.1 0 − 0.2 − 0.2 − 0.4 − 0.3 − 0.6 − 0.4 − 0.8 0 − 0.5 − 1 (c) Discontinuous Galerkin approximation with a low penalty parameter

Why bother with DG? Can accommodate hanging nodes Useful in adaptive mesh refinement

Efficiency and accuracy 10 Conforming DG Displacement Error (m) ten ¡Eyck ¡and ¡Lew ¡2010 ¡ ¡ (nonlinear ¡elas7city) ¡ 1 0.1 1 10 100 1000 10000 Total DOF Figure 6. Plot of the L 2 -norm over B 0 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 overlap, indicating that both methods provide the same accuracy for the same computational cost, i.e. for this example they are equally efficient.

Our objective • 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

First, review derivation of heat equation heat flux vector ¡ q = Heat equation : Balance of energy div q = s heat source ¡ s = + Fourier’s law q = � k r ϑ give the (steady) heat equation − k ∆ ϑ = s

Governing equations for elasticity displacement vector u Γ Ω x stress tensor or matrix σ ¡ u ( x ) Γ Elasticity: Equilibrium − div σ = f + σ = λ div u + ( r u + [ r u ] T ) Hooke’s law � ¯ give Navier’s equation ∆ u = � [ ∆ u + (1 + λ ) r div u ] = f

Equivalent to minimization problem Z Z | r s u | 2 + (1 + λ )(div u ) 2 dV � u J ( u ) = 1 min f · u dV 2 Ω Ω 2 [ r u + ( r u ) T ] r s u = 1 Z ✓ ∂ u i Define a ( u , v ) = [ r s u : r s v + (1 + λ )div u div v ] dx ( r s u ) ij = 1 ◆ + ∂ u j Ω 2 ∂ x j ∂ x i Z ` ( v ) = f · v dx Ω 1 Then the minimization problem is min 2 a ( v , v ) − ` ( v ) v ∈ V or equivalently a ( u , v ) = ` ( v ) which has a unique solution in 0 ( Ω )] d V := [ H 1 Furthermore (Brenner and Sung 1992) the solution satisfies k u k H 2 + λ k div u k H 1  C k f k L 2 u ∈ [ H 2 ( Ω )] 2

Locking The use of low-order elements leads to a non-physical solution in the incompressible limit Z a ( u , v ) = [ r s u : r s v + (1 + λ )div u div v ] dx Ω λ → ∞ div u → 0 Q 1 0,8 ¡ 0,7 ¡ Q 1 0,6 ¡ 0,5 ¡ Case ¡III ¡ 0,4 ¡ Case ¡V ¡(EAS) ¡ 0,3 ¡ MHW ¡Case ¡1 ¡ 0,2 ¡ HW ¡Case ¡I ¡or ¡Q1 ¡ 0,1 ¡ 0 ¡ -­‑0,6 ¡ -­‑0,4 ¡ -­‑0,2 ¡ 0 ¡ 0,2 ¡ 0,4 ¡ 0,6 ¡ -­‑0,1 ¡ -­‑0,2 ¡ We explore the use of DG methods as a remedy for locking