ANOTHER hp 3 D , ORIENTATION EMBEDDED SHAPE FUNCTIONS FOR ELEMENTS - - PowerPoint PPT Presentation

ANOTHER hp3D, ORIENTATION EMBEDDED SHAPE FUNCTIONS FOR ELEMENTS OF ALL SHAPES, AND HOW TO CODE THE PRIMAL DPG METHOD

• L. Demkowicz
• F. Fuentes, B. Keith and S. Nagaraj

ICES, The University of Texas at Austin ICERM Workshop on Robust Discretization and Fast Solvers for Computable Multi-Physics Models Brown University, May 11-15, 2014

Motto: The difference between 2D and 3D computations is analogous to those that separate boys from men.

Brown U, May 11-15, 2014 How to code primal DPG? 2 / 48

Act One One more hp3D code...

Brown U, May 11-15, 2014 How to code primal DPG? 3 / 48

hp3d ∗ overview

The code is designed for solving 3D coupled multi-physics problems: GMP provides MBG description with globally C0-compatible parameterizations. hp-Mesh stores all topological (i.e., vertex, edge, face, and interior) nodes and supports variable order of approximation; an element is defined as an assembly of nodes. Data struct. includes only two object arrays. ELEMS initial mesh elements only. NODES a family of nodal trees. FEs support for exact sequence elements; shape functions with embedded

• rientations.

Refinements support for isotropic and anisotropic refinements for elements of all shapes (i.e., tet, hexa, prism, and pyramid). Constrained Approx. allows one-level hanging nodes.

∗Developed with Paolo Gatto (now at Brown) and Kyongjoo Kim (now at Sandia).

Brown U, May 11-15, 2014 How to code primal DPG? 4 / 48

Mesh topology

hp-mesh structure includes all topological entities

◮ For multi-physics problem, one may need to support different energy spaces following the

exact sequence.

◮ All topological nodes are defined with orientations.

Figure: Different discretizations associate DOF with different topological nodes.

Brown U, May 11-15, 2014 How to code primal DPG? 5 / 48

Mesh topology

hp-mesh structure includes all topological entities

◮ For multi-physics problem, one may need to support different energy spaces following the

exact sequence.

◮ All topological nodes are defined with orientations. ◮ Mesh refinement translates into an orchestrated node breaking procedure.

Element = Vertices + Edges + Faces + Interior

Brown U, May 11-15, 2014 How to code primal DPG? 5 / 48

Nodal trees

◮ Mesh refinements are explicitly recorded in nodal trees. ◮ No need for additional search tables (e.g., hash-table) to keep track of dynamically

changed element-to-node connectivities.

◮ Generic algorithms are developed to reconstruct element-to-node connectivities.

Figure: A family of nodal trees are constructed during mesh refinements.

Brown U, May 11-15, 2014 How to code primal DPG? 6 / 48

Act Two Shape functions for elements of all shapes and the exact sequence spaces

Brown U, May 11-15, 2014 How to code primal DPG? 7 / 48

Elements of “all shapes”

Brown U, May 11-15, 2014 How to code primal DPG? 8 / 48

Element system of coordinates, node enumeration and local orientations

Hexahedral element vertices:

Brown U, May 11-15, 2014 How to code primal DPG? 9 / 48

Element system of coordinates, node enumeration and local orientations

Hexahedral element edges. Enumeration and local orientation (parametrization):

Brown U, May 11-15, 2014 How to code primal DPG? 9 / 48

Element system of coordinates, node enumeration and local orientations

Hexahedral element edge nodes enumeration :

Brown U, May 11-15, 2014 How to code primal DPG? 9 / 48

Element system of coordinates, node enumeration and local orientations

Hexahedral element face nodes enumeration and orientation (parametrization):

Brown U, May 11-15, 2014 How to code primal DPG? 9 / 48

Element system of coordinates, node enumeration and local orientations

Tetrahedral element vertex and edge nodes enumeration and orientation: and so on... (Show module/element data.F)

Brown U, May 11-15, 2014 How to code primal DPG? 9 / 48

Finite element spaces form exact sequences of the first type †

segment: Pp

∂

− → Pp−1 quad: Pp ⊗ Pq ∇ − → (Pp−1 ⊗ Pq) × (Pp ⊗ Pq−1)

curl

− → Pp−1 ⊗ Pq−1 triangle: Pp ∇ − → Pp−1 ⊕ Rp

curl

− → Pp−1 where: Rp : {E ∈ ˜ Pp : x · E(x) = 0 ∀x} hexa:

Pp ⊗ Pq ⊗ Pr ∇ − → (Pp−1 ⊗ Pq ⊗ Pr) × (Pp ⊗ Pq−1 ⊗ Pr) × (Pp ⊗ Pq ⊗ Pr−1) ∇× − → (Pp ⊗ Pq−1 ⊗ Pr−1) × (Pp−1 ⊗ Pq ⊗ Pr−1) × (Pp−1 ⊗ Pq−1 ⊗ Pr) ∇◦ − → Pp−1 ⊗ Pq−1 ⊗ Pr−1

and so on...

†Refers to N´

ed´ elec’s elements of the first type. The drop in polynomial degree is always one.

Brown U, May 11-15, 2014 How to code primal DPG? 10 / 48

For the pyramid the spaces are far from trivial ‡

x =

ξ 1+ζ

y =

η 1+ζ

z =

ζ 1+ζ

ξ =

x 1−z

η =

y 1−z

ζ =

z 1−z

Wp = {u ∈ Qp,p,p p : ∇u ∈ Qp−1,p,p−1 p × Qp,p−1,p−1 p × Qp,p,p−1 p+1 } Qp = {E ∈ Qp−1,p,p−1 p+1 × Qp,p−1,p−1 p+1 × Qp,p,p−2 p+1 : ∇ × E ∈ Qp,p−1,p−1 p+2 × Qp−1,p,p−1 p+2 × Qp−1,p−1,p p+2 } Vp = {V ∈ Qp,p−1,p−1 p+2 × Qp−1,p,p−1 p+2 × Qp−1,p−1,p p+2 : ∇ · V ∈ Qp−1,p−1,p p+3 } Yp = Qp−1,p−1,p p+3

where

Qpqr k := { xiyj zl (1 + z)k : 0 ≤ i ≤ p, 0 ≤ j ≤ q, 0 ≤ l ≤ r }

‡Nigam, N. and Phillips, J., “High-order Conforming Finite Elements on Pyramids”, IMA

Journal of Numerical Analysis, 32(2): 448-483, 2012.

Brown U, May 11-15, 2014 How to code primal DPG? 11 / 48

A bit of history: Ciarlet’s approach.

Degrees of freedom come first: ψj : X(K) ⊃ X(K) → I R, j = 1, . . . , n where X(K) is continuosly embedded in the energy space, X(K) is the space of FE shape functions, d.o.f. functionals are continuous, n = dim = X(K), and ψj|X(K) are linearly independent (unisolvency condition). Shape functions come next (dual basis): φi ∈ X(K), ψj, φi = δij Global conformity (continuity) is enforced by a proper selection of d.o.f.. The approach implies standard assembly procedure. Shape functions are precomputed as linear combinations of monomials.

Brown U, May 11-15, 2014 How to code primal DPG? 12 / 48

A bit of history: Szabo’s approach.

Shape functions come first: φi ∈ X(K), i = 1, . . . , n = dimX(K) Global conformity (continuity) is enforced by classifying the shape functions into vertex, edge, face and interior “modes” (Szabo’s terminology). The assembly procedure is non-standard as it requires accounting for edge and face orientations (Szabo’s sign factors). Shape functions are computed on the fly using recursive formulas. No d.o.f. are introduced.

Brown U, May 11-15, 2014 How to code primal DPG? 13 / 48

Topological character of shape functions

H1 vertex, edge, face, interior functions H(curl) edge, face, interior functions H(div) face, interior functions L2 interior functions

Brown U, May 11-15, 2014 How to code primal DPG? 14 / 48

Orientation embedded shape functions

§. Edge shape

functions are owned by the edge.

Switch from edge coordinate ξedge to local edge coordinate t: ˆ u(t) = u(ξedge), t = 1 − ξedge

§P. Gatto, L. Demkowicz, “Construction of H1-Conforming Hierarchical Shape Functions for

Elements of All Shapes and Transfinite Interpolation”, Finite Elements in Analysis and Design, 46: 474-486, 2010.

Brown U, May 11-15, 2014 How to code primal DPG? 15 / 48

Orientation embedded shape functions . Edge shape functions are owned by the edge.

Project ξ onto the edge: pick up the value at t = Pξ, and blend with a blending function ψ(ξ): u(ξ) := ˆ u(t)ψ(ξ) transform → project → blend (Show the code.)

Brown U, May 11-15, 2014 How to code primal DPG? 15 / 48

In a triangle we need to blend with a higher order polynomial (Dubinger’s construction)

It is convenient to work with edge affine coordinates λ1, λ2: φ(ξ) = ˆ un( λ1 λ1 + λ2

• µ1

, λ2 λ1 + λ2

• µ2

) (λ1 + λ2)n

• blending function

where ˆ u(µ1, µ2) is one of possible representation of the edge shape function in terms of edge affine coordinates µ1, µ2.

Brown U, May 11-15, 2014 How to code primal DPG? 16 / 48

This leads to the concept of scaled polynomials §

Scaled polynomials: un ∈ Pn(−1, 1), ˜ un(x; t) := un(x t )tn Scaling produces a homogeneous polynomials of (total) degree n in (x, t). Scaled Legendre Polynomials are evaluated using the standard recursion formula: n ˜ Pn(x; t) = (2n − 1)x ˜ Pn−1(x; t) − (n − 1)t2 ˜ Pn−2(x; t) . Shifted scaled Legendre polynomials are invariant (by construction) with affine transformations: transform to the symmetric triangle → scale → backtransform (shift)

§J. Schoeberl, J. Zaglmayr, “High order N´

ed´ elec Elements with Local Complete Sequence Property”, International Journal for Computation and Mathematics in Electrical and Electronic Engineering, 24(2): 374-384, 2005.

Brown U, May 11-15, 2014 How to code primal DPG? 17 / 48

Shifted scaled polynomials

˜ P s

n(x′; t′) := ˜

Pn(2x′ − t′; t′) Shifted Scaled Lobatto Polynomials are integrals of shifted scaled Legendre polynomials ˜ Ls

n+1(x′; t′) :=

x′ ˜ P s

n(s′; t′) ds′

and are computed using recursion.

Brown U, May 11-15, 2014 How to code primal DPG? 18 / 48

Edge functions revisited.

Homogenization of polynomials: [ˆ u]n(λ1, λ2) := ˆ un(λ1, λ2; λ1 + λ2) = ˆ un(

λ1 λ1+λ2 , λ2 λ1+λ2 )(λ1 + λ2)n

Edge H1 shape function again: φ = [ˆ u]n(λ1, λ2) Fastforward to edge H(curl) shape function: E = [ˆ u]n(λ1, λ2) (λ1∇λ2 − λ2∇λ1)

• Whitney edge shape function

Critical property: ∇ × E = (n + 2)[ˆ u]n(λ1, λ2)∇λ1 × ∇λ2 (no derivatives of un needed !).

Brown U, May 11-15, 2014 How to code primal DPG? 19 / 48

An analogous property for H(div) shape functions

Face shape function: V = [ˆ u]n(λ1, λ2, λ3) (λ1∇λ2 × ∇λ3 + λ3∇λ1 × ∇λ2 + λ2∇λ3 × ∇λ1) Divergence of V : ∇ ◦ V = (n + 3)[ˆ u]n(λ1, λ2, λ3) [∇λ1, ∇λ1, ∇λ3]

• mixed product

requires no derivatives of un. Orientations are accounted for by swapping affine coordinates (Show the code).

Brown U, May 11-15, 2014 How to code primal DPG? 20 / 48

Tones of details skipped

In the incoming ICES Report: Construction of H1 face and element bubbles, Construction of H(curl) face and element bubbles, Construction of H(div) bubbles, Construction of shape functions for the prism (tensor product of a triangle and interval), etc.

Brown U, May 11-15, 2014 How to code primal DPG? 21 / 48

Tricky pyramid

x =

ξ 1+ζ

y =

η 1+ζ

z =

ζ 1+ζ

ξ =

x 1−z

η =

y 1−z

ζ =

z 1−z

Vertex shape functions: λ1 = (1−x)(1−y)

1+z

= (1−ξ−ζ)(1−η−ζ)

1−ζ

λ2 = x(1−y)

1+z

= ξ(1−η−ζ)

1−ζ

λ3 =

xy 1+z

=

ξη 1−ζ

λ4 = (1−x)y

1+z

= (1−ξ−ζ)η

1−ζ

λ5 =

z 1+z

= ζ

Brown U, May 11-15, 2014 How to code primal DPG? 22 / 48

Tricky pyramid

Some of the shape functions can be constructed using the formulas for the tetrahedron with the affine coordinates replaced with vertex shape functions. For instance, for vertical edge 15 (common to two triangular faces), H1 edge modes: φ = [ˆ u]n(λ1, λ5) H(curl) modes: E = [ˆ u]n(λ1∇λ5 − λ5∇λ1) For edge 12, however, the construction is different as it needs to reconcile the triangular vertical face and horizontal square base face shape functions.

Brown U, May 11-15, 2014 How to code primal DPG? 23 / 48

Tricky pyramid; a stumbling block

For all but one case, we were able to construct the edge and face modes as products of properly extended higher order shape functions defined on the edge or face and the lowest

• rder edge or shape function. This is a unique property of N´

ed´ elec spaces of first type. The stumbling block was the H(div) face 125 modes.

Brown U, May 11-15, 2014 How to code primal DPG? 24 / 48

Pull back maps (Piola transforms) and the original space

  Vξ Vη Vζ   = (1 + z)2   1 + z −x 1 + z −y 1     Vx Vy Vz     Vx Vy Vz   = − 1 (1 + z)3   1 x 1 y −(1 + z)     Vξ Vη Vζ   {V ∈ Qk,k−1,k−1

k+2

× Qk−1,k,k−1

k+2

× Qk−1,k−1,k

k+2

: ∇ · V ∈ Qk−1,k−1,k−1

k+3

}

Brown U, May 11-15, 2014 How to code primal DPG? 25 / 48

Lowest order mode¶

V = 1 (1 + z)3   1 − y − 1

2z

  ∇ · V = − 3 2(1 + z)4 . Higher order mode, first attempt: V =

1 (1+z)3

  (1 − y) uk−1(

x 1+z , z 1+z )

− C

2 z

  , ∇ · V = −

1 (1+z)4

• (uk−1(

x 1+z , z 1+z ) − C)(1 + z) + 3C 2

• where C = uk−1(0, 1).

¶At a first glance, presence of non-zero z-component seems to be unnecessary. With

vanishing Vz, the shape function would still have had the right vanishing properties in the x, y, z

• space. It is only upon computing the divergence that we learn that the third component is

necessary to fit the function into the right space.

Brown U, May 11-15, 2014 How to code primal DPG? 26 / 48

Our solution: add a bubble

The main source of the trouble seems to come from the fact that the H(div) space is a bit “too tied”. The problem can be solved by adding the bubble V = 1 (1 + z)3   z   ∇ · V = 1 − 2z (1 + z)4

• r, in the finite pyramid frame,

V = −ζ  

ξ 1−ζ η 1−ζ

−1   ∇ · V = 1 − 2 ζ 1 − ζ . The spaces are extended accordingly for the higher order element, {V ∈ Qk,k−1,k−1

k+2

× Qk−1,k,k−1

k+2

× Qk−1,k−1,k

k+2

: ∇ · V ∈ Qk−1,k−1,k

k+3

} .

Brown U, May 11-15, 2014 How to code primal DPG? 27 / 48

Higher order modes

V = uk−1(

x 1+z , z 1+z )

(1 + z)3   1 − y  

• r, in the finite pyramid frame,

V = u(ξ, ζ)  

1−η−ζ 1−ζ

  div V = −u(ξ, ζ) 1 − ζ .

Brown U, May 11-15, 2014 How to code primal DPG? 28 / 48

A summary

◮ A logically coherent construction for elements of all shapes and the whole

exact sequence.

◮ Small number of primitives (scaled polynomials) and a hierarchical logic

allowing for quick modifications.

◮ Stand alone package in Fortran 90 (6k lines). ◮ Extended element-to-nodes connectivity needed, standard assembly

procedure.

◮ Partial verification (reproduction of polynomials) done. ◮ Anyone interested in helping with verification ?

Brown U, May 11-15, 2014 How to code primal DPG? 29 / 48

Act Three Primal DPG Method

Brown U, May 11-15, 2014 How to code primal DPG? 30 / 48

Primal DPG method

Standard assumptions: Ω ⊂ I RN Lipschitz domain, Elements:K Faces (Edges):e Skeleton:Γh =

K ∂K

Internal skeleton:Γ0

h = Γh − Γ

Brown U, May 11-15, 2014 How to code primal DPG? 31 / 48

Primal DPG method

Given f ∈ L2(Ω), consider the model problem,

• u

= u0

• n Γ := ∂Ω

−∆u = f in Ω Multiply the PDE with a test function v, integrate over each element K, integrate by parts and sum up over all elements,

• K
• K

∇u · ∇v +

• K
• ∂K

∂u ∂n v =

• K
• f v

The boundary term represents jumps,

• K
• ∂K

∂u ∂n v =

• e⊂Γ0

h

• e

∂u ∂ne [v] +

• e⊂Γ
• e

∂u ∂ne v

Brown U, May 11-15, 2014 How to code primal DPG? 32 / 48

Primal DPG method

This leads to the variational problem: u ∈ H1(Ω), u = u0 on Γ, ˆ t ∈ H−1/2(Γh) (∇u, ∇hv) − ˆ t, vΓh = (f, v) v ∈ H1(Ωh) where H−1/2(Γh) = trace of H(div, Ω) on Γh equipped with the quotient norm. Theorem The variational problem above is well posed with a mesh independent inf-sup constant γ.

• L. Demkowicz and J. Gopalakrishnan. “A primal DPG method without a first order

reformulation”, Comput. Math. Appl., 66(6):1058–1064, 2013

Brown U, May 11-15, 2014 How to code primal DPG? 33 / 48

DPG is a minimum residual method ∗∗

u ∈ U b(u, v) = l(v) v ∈ V ⇔ Bu = l B : U → V ′ Bu, v = b(u, v)

∗∗

◮ J.H. Bramble, R.D. Lazarov, J.E. Pasciak, “A Least-squares Approach Based on a Discrete

Minus One Inner Product for First Order Systems”Math. Comp, 66, 935-955, 1997.

◮ L.D., J. Gopalakrishnan. “A Class of Discontinuous Petrov-Galerkin Methods. Part II:

Optimal Test Functions.”Numer. Meth. Part. D. E., 27, 70-105, 2011.

Brown U, May 11-15, 2014 How to code primal DPG? 34 / 48

DPG is a minimum residual method ∗∗

u ∈ U b(u, v) = l(v) v ∈ V ⇔ Bu = l B : U → V ′ Bu, v = b(u, v)

◮ Minimum residual method: Uh ⊂ U, 1 2Buh − l2 V ′ → min uh∈Uh

∗∗

◮ J.H. Bramble, R.D. Lazarov, J.E. Pasciak, “A Least-squares Approach Based on a Discrete

Minus One Inner Product for First Order Systems”Math. Comp, 66, 935-955, 1997.

◮ L.D., J. Gopalakrishnan. “A Class of Discontinuous Petrov-Galerkin Methods. Part II:

Optimal Test Functions.”Numer. Meth. Part. D. E., 27, 70-105, 2011.

Brown U, May 11-15, 2014 How to code primal DPG? 34 / 48

DPG is a minimum residual method ∗∗

u ∈ U b(u, v) = l(v) v ∈ V ⇔ Bu = l B : U → V ′ Bu, v = b(u, v)

◮ Minimum residual method: Uh ⊂ U, 1 2Buh − l2 V ′ → min uh∈Uh ◮ Riesz operator:

RV : V → V ′, RV v, δv = (v, δv)V is an isometry, RV vV ′ = vV .

∗∗

◮ J.H. Bramble, R.D. Lazarov, J.E. Pasciak, “A Least-squares Approach Based on a Discrete

Minus One Inner Product for First Order Systems”Math. Comp, 66, 935-955, 1997.

◮ L.D., J. Gopalakrishnan. “A Class of Discontinuous Petrov-Galerkin Methods. Part II:

Optimal Test Functions.”Numer. Meth. Part. D. E., 27, 70-105, 2011.

Brown U, May 11-15, 2014 How to code primal DPG? 34 / 48

DPG is a minimum residual method ∗∗

u ∈ U b(u, v) = l(v) v ∈ V ⇔ Bu = l B : U → V ′ Bu, v = b(u, v)

◮ Minimum residual method: Uh ⊂ U, 1 2Buh − l2 V ′ → min uh∈Uh ◮ Riesz operator:

RV : V → V ′, RV v, δv = (v, δv)V is an isometry, RV vV ′ = vV .

◮ Minimum residual method reformulated: 1 2Buh − l2 V ′ = 1 2R−1 V (Buh − l)2 V → min uh∈Uh

∗∗

◮ J.H. Bramble, R.D. Lazarov, J.E. Pasciak, “A Least-squares Approach Based on a Discrete

Minus One Inner Product for First Order Systems”Math. Comp, 66, 935-955, 1997.

◮ L.D., J. Gopalakrishnan. “A Class of Discontinuous Petrov-Galerkin Methods. Part II:

Optimal Test Functions.”Numer. Meth. Part. D. E., 27, 70-105, 2011.

Brown U, May 11-15, 2014 How to code primal DPG? 34 / 48

(Practical) DPG is an easy mixed method††

Introducing, ( R−1

V (Buh − l)

• =:ψ(error representation function)

, R−1

V Bδuh)V = 0

δuh ∈ Uh

††W. Dahmen, Ch. Huang, Ch. Schwab, and G. Welper. “Adaptive Petrov Galerkin methods

for first order transport equations”, SIAM J. Num. Anal. 50(5): 242-2445, 2012

Brown U, May 11-15, 2014 How to code primal DPG? 35 / 48

(Practical) DPG is an easy mixed method††

Introducing, ( R−1

V (Buh − l)

• =:ψ(error representation function)

, R−1

V Bδuh)V = 0

δuh ∈ Uh

• r

ψ = R−1

V (Buh − l)

(ψ, R−1

V Bδuh)V = 0

δuh ∈ Uh

††W. Dahmen, Ch. Huang, Ch. Schwab, and G. Welper. “Adaptive Petrov Galerkin methods

for first order transport equations”, SIAM J. Num. Anal. 50(5): 242-2445, 2012

Brown U, May 11-15, 2014 How to code primal DPG? 35 / 48

(Practical) DPG is an easy mixed method††

Introducing, ( R−1

V (Buh − l)

• =:ψ(error representation function)

, R−1

V Bδuh)V = 0

δuh ∈ Uh

• r

(ψ, δv)V − b(uh, δv) = −l(δv) ∀δv ∈ V b(δuh, ψ) = 0 ∀δuh ∈ Uh

††W. Dahmen, Ch. Huang, Ch. Schwab, and G. Welper. “Adaptive Petrov Galerkin methods

for first order transport equations”, SIAM J. Num. Anal. 50(5): 242-2445, 2012

Brown U, May 11-15, 2014 How to code primal DPG? 35 / 48

(Practical) DPG is an easy mixed method††

Introducing, ( R−1

V (Buh − l)

• =:ψ(error representation function)

, R−1

V Bδuh)V = 0

δuh ∈ Uh

• r

(ψ, δv)V − b(uh, δv) = −l(δv) ∀δv ∈ V b(δuh, ψ) = 0 ∀δuh ∈ Uh In practice, error representation function ψ must be approximated within an “enriched test space” ψh ∈ ˜ Vh ⊂ V (ψh, δvh)V − b(uh, δvh) = −l(δvh) ∀δvh ∈ ˜ Vh b(δuh, ψh) = 0 ∀δuh ∈ Uh

††W. Dahmen, Ch. Huang, Ch. Schwab, and G. Welper. “Adaptive Petrov Galerkin methods

for first order transport equations”, SIAM J. Num. Anal. 50(5): 242-2445, 2012

Brown U, May 11-15, 2014 How to code primal DPG? 35 / 48

Brezzi’s Conditions

◮ inf-sup in the kernel - trivially satisfied ◮ LBB inf-sup condition:

sup

˜ vh∈ ˜ Vh

|b(uh, ˜ vh)| ˜ vhV ≥ γhuhU , is “easy” to satisfy if we take “large” enriched space ˜ Vh ⊂ V . The proof relies on the construction of an appropriate Fortin operator ‡‡ .The Fortin constant affects the ultimate discrete inf-sup constant γh that no longer is equal to γ.

‡‡J. Gopalakrishnan and W. Qiu. “An analysis of the practical DPG method”, Math. Comp.,

2012.

Brown U, May 11-15, 2014 How to code primal DPG? 36 / 48

Brezzi’s Conditions

◮ inf-sup in the kernel - trivially satisfied ◮ LBB inf-sup condition:

sup

˜ vh∈ ˜ Vh

|b(uh, ˜ vh)| ˜ vhV ≥ γhuhU , is “easy” to satisfy if we take “large” enriched space ˜ Vh ⊂ V . The proof relies on the construction of an appropriate Fortin operator ‡‡ .The Fortin constant affects the ultimate discrete inf-sup constant γh that no longer is equal to γ. The approximate mixed problem can be viewed as the discretization

• f two infinite-dimensional mixed problems solved for:

uh, ψ (a-posteriori error estimation), and u, ψ = 0 (a-priori estimates). See Jay’s talk.

‡‡J. Gopalakrishnan and W. Qiu. “An analysis of the practical DPG method”, Math. Comp.,

2012.

Brown U, May 11-15, 2014 How to code primal DPG? 36 / 48

Main point: ψ is Condensed Out Elementwise

uh =

N

• i=1

uiei, ψh ≈

M

• j=1

ψjgj, M >> N Matrix form:

• G

−B BT ψ u

• =
• −b
• where

Gij

• Gram matrix

= (gi, gj), Bij

• expanded stiffness matrix

= b(ei, gj) . For broken test space, G is element block-diagonal, and ψh is eliminated elemenwise, BT G−1B

• DPG stiffness matrix

uh = BT G−1b

• DPG load vector

. DPG can be viewed as a preconditioned least squares method.

Brown U, May 11-15, 2014 How to code primal DPG? 37 / 48

Primal DPG Formulation

Group unknown (watch for the overloaded symbol): uh := ( uh

• field

, ˆ th

• flux

) Mixed system:   G −B1 −B2 BT

1

BT

2

    ψ u ˆ t   =   −b   where B1, B2 correspond to (∇uh, ∇hvh) and −ˆ th, vh, resp. Eliminate ψ to get the DPG system: BT

1 G−1B1

BT

1 G−1B2

BT

2 G−1B1

BT

2 G−1B2

u ˆ t

• =

BT

1 G−1b

BT

2 G−1b

• Brown U, May 11-15, 2014

How to code primal DPG? 38 / 48

Primal DPG method

Neglecting the error steming from the approximation of optimal test functions (computation of residual), we have,

• u − uh2

H1(Ω)

+ˆ t − ˆ th2

H−1/2(Γh)

1/2 ≤ 1

γ

inf

wh,rh

• u − wh2

H1(Ω) + ˆ

t − rh2

H−1/2(Γh)

1/2

• best approximation error

Additionally,

• u − uh2

H1(Ω)

+ˆ t − ˆ th2

H−1/2(Γh)

1/2 ≤ 1

γ

sup

v∈H1(Ωh)

|(∇uh, ∇hv) − ˆ th, vΓh| vH1(Ωh)

• computable residual

Brown U, May 11-15, 2014 How to code primal DPG? 39 / 48

FE discretization

Hexahedral meshes H1 element for field uh: Pp ⊗ Pp ⊗ Pp, Trace of H(div) element: (Pp ⊗ Pp−1 ⊗ Pp−1) × (Pp−1 ⊗ Pp ⊗ Pp−1) × (Pp−1 ⊗ Pp−1 ⊗ Pp) for flux ˆ th, and the enriched element: Pp+∆p ⊗ Pp+∆p ⊗ Pp+∆p, for test function vh. In reported experiments: p = 1, 2, 3, ∆p = 2.

Brown U, May 11-15, 2014 How to code primal DPG? 40 / 48

Smooth solution, uniform refinements

Rectangular domain Ω = (0, 1) × (0, 2) × (0, 1), Smooth solution: u = sin πx sin πy sin πz Boundary condition: u = 0. Residual versus H1 error.

Brown U, May 11-15, 2014 How to code primal DPG? 41 / 48

Manufactured shock solution

BC: u = u0.

Brown U, May 11-15, 2014 How to code primal DPG? 42 / 48

Shock solution, uniform and h-adaptive refinements, p = 1

Convergence history for the residual and H1 error

Brown U, May 11-15, 2014 How to code primal DPG? 43 / 48

Shock solution, uniform and h-adaptive refinements, p = 2

Convergence history for the residual and H1 error

Brown U, May 11-15, 2014 How to code primal DPG? 44 / 48

Shock solution, uniform and h-adaptive refinements, p = 3

Convergence history for the residual and H1 error

Brown U, May 11-15, 2014 How to code primal DPG? 45 / 48

Shock solution, p = 3, Mixed BC

Mixed BC: trace: bottom, top, flux: sides. Convergence history for the residual and H1 error

Brown U, May 11-15, 2014 How to code primal DPG? 46 / 48

Reaction-dominated difusion, p = 2.

• u

= 0

• n Γ

−ǫ2∆u + u = 1 in Ω ǫ = 0.01, left: solution after 7 iterations, right: convergence history

Brown U, May 11-15, 2014 How to code primal DPG? 47 / 48

Convection-dominated difusion, p = 2.

   −ǫ2∆u − u = sin πy sin πz at x = 0 u = 0

• n the rest of Γ

−ǫ2∆u + ∂u

∂x

= 0 in Ω ǫ = 0.01, left: solution after 5 iterations, right: convergence history

Brown U, May 11-15, 2014 How to code primal DPG? 48 / 48

Concluding remarks

◮ Incorporation of orientations into the element shape functions routine

dramatically simplifies assembly procedure and constrained approximation.

◮ Coding DPG within a code supporting the exact sequence elements is very

straightforward.

◮ With ultraweak variational formulation there is no need to construct 3D

conforming discretizations - DPG naturally extends to polyhedral elements.

Thank you

Brown U, May 11-15, 2014 How to code primal DPG? 49 / 48