The geometric basis of mimetic spectral approximations

(1)

The geometric basis of mimetic spectral

approximations

Marc Gerritsma

Ren´

e Hiemstra, Jasper Kreeft, Artur Palha and Pedro Rebelo

Faculty of Aerospace Engineering

TU Delft, The Netherlands

Discretization Methods for Polygonal and Polyhedral Meshes

(2)

Introduction Physics & Geometry Integral values Some results

Geometry Orientation

Relation between physical variables and geometric objects

All physical variables are related to geometric objects

Mass

M is associated to a

volume

, V .

Average density

ρ =

¯

M_V

. ρ = lim

V →0M_V

mathematically justified, but physically

questionable.

(3)

Relation between physical variables and geometric objects

All physical variables are related to geometric objects

Flux

F is associated to a

area

, A.

Average flux density

F /A. lim

A→0F_A

mathematically justified, but physically

question-able.

(4)

Geometry Orientation

Relation between physical variables and geometric objects

All physical variables are related to geometric objects

Velocity

F is associated to a

curve

, C.

The velocity, v, can be measured by recording the position, r, of a particle at two

consecutive time instants, t

1

and t

2

. These positions are related to the velocity by

r(t

2

) − r(t

1

) =

Z

t2 t1

dr

dt

dt =

Z

t2 t1

v dt

This relation is

exact

and then we approximate ‘the’ velocity by

v ≈

r(t

2

) − r(t

1

)

t

2

− t

1

(5)

Relation between geometric objects

Boundary operator

The most important operator in mimetic methods is the

boundary operator ∂

∂ : k-dim −→ (k − 1)-dim

(6)

Geometry

Orientation

Relation between geometric objects

Boundary operator

The most important operator in mimetic methods is the

boundary operator ∂

∂ : k-dim −→ (k − 1)-dim

(7)

Relation between geometric objects

Boundary operator

The most important operator in mimetic methods is the

boundary operator ∂

∂ : k-dim −→ (k − 1)-dim

(8)

Geometry

Orientation

Orientation and type of orientation

Orientation and sense of orientation

Every geometric object can be

oriented in two ways

. For instance, in a surface we define

a

sense of rotation

, either clockwise or counter clockwise

(9)

Orientation and type of orientation

Orientation and sense of orientation

Every geometric object can be

oriented in two ways

. For instance, in a surface we define

a

sense of rotation

, either clockwise or counter clockwise

(10)

Geometry

Orientation

Orientation and type of orientation

∂ and ?∂?

Let

?

denote the

operator which switches between inner- and outer-orientation

∂

∂ *

*

Then we have the operations:

(11)

Orientation and type of orientation

∂ and ?∂?

Let

?

denote the

operator which switches between inner- and outer-orientation

∂

∂ *

*

Then we have the operations:

(12)

Geometry

Orientation

Oriented dual cell complexes

Double boundary complex

In 3D we have

points

,

curves

,

surfaces

and

volumes

Outer Orientation

∂ ∂ ∂

*

∂ ∂ ∂

(13)

Matrix representation of boundary operator

Set of points:







P

1

P

2

P

3

P

4







P

₁

P

₂

P

₃

P

₄

(14)

Geometry

Orientation

Matrix representation of boundary operator

Set of lines:







∂L

1

∂L

2

∂L

3

∂L

4







=







−1

1

0

0 −1

1 −1

0

1

0

0 −1

0

1 











P

1

P

2

P

3

P

4







P

1 P

2 P

₃

P

₄

L

₃

L

4 L

2 L

₁

(15)

Matrix representation of boundary operator

Surface:

∂∂S

1

=

1 −1

−1

1 





∂L

1

∂L

2

∂L

3

∂L

4







P

1 P

2 P

₃

P

₄

L

₃

L

4 L

2 L

₁

S

1

(16)

Geometry

Orientation

Matrix representation of boundary operator

Surface:

∂∂S

1

=

1 −1

−1

1 





∂L

1

∂L

2

∂L

3

∂L

4







=

1 −1

−1

1 





−1

1

0

0 −1

1 −1

0

1

0

0 −1

0

1 





|

{z

}

=h

0

i







P

1

P

2

P

3

P

4







P

1 P

2 P

₃

P

₄

L

₃

L

4 L

2 L

₁

S

1

(17)

Geometry

Orientation

Final remarks geometric objects

Topological vs metric-dependent operations

(18)

Geometry

Orientation

Final remarks geometric objects

Topological vs metric-dependent operations

The boundary operator

∂

is

topological operator

,

?

operator is

metric-dependent

.

Nilpotency of ∂ and ?∂?

Application of the boundary operator twice always yields the empty set:

∂ ◦ ∂ ≡ 0

1 2 3 4 5 6 7 9 ₁₁ 10 ₁₂ 8 2 3 4 6 5

(19)

Assigning a value to geometric objects

k-chains and k-cochains

A basic k-dimensional object will be called a

k-cell

, τ

k

. A collection of oriented k-cells

is called a

k-chain

, c

k

. The space of all k-chains will be denoted by C

k

The operation which assigns a value to a physical quantity associated with a geometric

object is called a

k-cochain

, c

k

_:

c

k

: C

k

−→ R

⇐⇒

D

c

k

, c

k

E

∈ R

(20)

Assigning values to geometric objects

Mother of all equations The ugly stepmother

Assigning a value to geometric objects

k-cochains and integration

c

k

: C

k

−→ R

⇐⇒

D

c

k

, c

k

E

∈ R

In the continuous setting in 3D this should be compared to

k = 0 , point : f (P ) ,

k = 1 , curve :

Z

C

a(x, y, z) dx + b(x, y, z) dy + c(x, y, z) dz ,

k = 2 , surface :

Z

S

P (x, y, z) dydz + Q(x, y, z) dzdx + R(x, y, z) dxdy ,

k = 3 , volume :

Z

V

ρ(x, y, z) dxdydz .

The expression underneath the integral sign is called a

differential k-form, a

(k)

_.

D

a

(k)

, Ω

k

E

:=

Z

Ωk

a

(k)

∈ R

(21)

Assigning a value to geometric objects

k-cochains and integration

c

k

: C

k

−→ R

⇐⇒

D

c

k

, c

k

E

∈ R

In the continuous setting in 3D this should be compared to

k = 0 , point :

f (P )

,

k = 1 , curve :

Z

C

a(x, y, z) dx + b(x, y, z) dy + c(x, y, z) dz

,

k = 2 , surface :

Z

S

P (x, y, z) dydz + Q(x, y, z) dzdx + R(x, y, z) dxdy

,

k = 3 , volume :

Z

V

ρ(x, y, z) dxdydz

.

The expression underneath the integral sign is called a

differential k-form, a

(k)

_.

D

a

(k)

, Ω

k

E

:=

Z

Ωk

a

(k)

∈ R

(22)

Mother of all equations The ugly stepmother

Assigning a value to geometric objects

k-cochains and integration

Both integration of differential forms and duality pairing between cochains and chains

is a

metric-free

operation

D

c

k

_{, c}

k

E

∈ R

⇐⇒

D

a

(k)

_{, Ω}

k

E

∈ R

c

k

_{, c}

(23)

Cell complex ⇔ computational grid

Topological mesh

If we glue volumes, surfaces, lines and points together we obtain a so-called

cell-complex

.

Manifold

Cell complex

(24)

Mother of all equations

The ugly stepmother

The Mother of all equations

The coboundary operator

Duality pairing between chains and cochains allows us to define

the adjoint of the

boundary operator δ

D

δc

k

, c

k+1

E

:=

D

c

k

, ∂c

k+1

E

The coboundary operator

maps k-cochains into (k + 1)-cochains

:

δ : C

k

−→ C

k+1

(25)

The Mother of all equations

The coboundary operator

D

δc

k

, c

k+1

E

:=

D

c

k

, ∂c

k+1

E

(26)

The ugly stepmother

The Mother of all equations

The coboundary operator

D

δc

k

, c

k+1

E

:=

D

c

k

, ∂c

k+1

E

Let C be an arbitrary curve going from the point A to the point B

k = 0 :

Z

C

grad φ d~

s =

Z

∂C

φ = φ(B) − φ(A)

(27)

The Mother of all equations

The coboundary operator

D

δc

k

, c

k+1

E

:=

D

c

k

, ∂c

k+1

E

Let S be a surface bounded by ∂S then

k = 1 :

Z

S

curl ~

A d ~

S =

Z

∂S

~

A · d~

s

(28)

The ugly stepmother

The Mother of all equations

The coboundary operator

D

δc

k

, c

k+1

E

:=

D

c

k

, ∂c

k+1

E

Let V be a volume, bounded by ∂V then

Z

V

div ~

F dV =

Z

∂V

~

F · d ~

S

(29)

The ugly stepmother

The Mother of all equations

The coboundary operator

D

δc

k

, c

k+1

E

:=

D

c

k

, ∂c

k+1

E

Duality pairing and the boundary operator DEFINE the coboundary operator!

I.e. grad, curl and div are defined through the topological relations and are therefore

coordinate-free

and

metric-free

.

(30)

The ugly stepmother

The Mother of all equations

The coboundary operator

D

δc

k

, c

k+1

E

:=

D

c

k

, ∂c

k+1

E

Duality pairing and the boundary operator DEFINE the coboundary operator!

I.e. grad, curl and div are defined through the topological relations and are therefore

coordinate-free

and

metric-free

.

If we choose basis functions for our numerical method, the

basis functions should cancel

from the equations

. There

cannot be an explicit dependence on the basis functions

.

The same topological relations

hold for low order methods and high order methods

.

(31)

The Mother of all equations

The coboundary operator

D

δc

k

, c

k+1

E

:=

D

c

k

, ∂c

k+1

E

At the continuous level, in terms of differential forms, this relation is given by the

generalized Stokes Theorem

Z

Ω_k+1

dω

(k)

:=

Z

∂Ω_k+1

ω

(k)

(32)

Assigning values to geometric objects Mother of all equations

The ugly stepmother

The ’Hodge-

?’ operator

The ’Hodge-?’ operator

Remember that ? was the operator which switches between

inner- and outer orientation

.

We can also write down a formal adjoint of this operation

D

?c

k

, c

n−k

E

:=

D

c

n−k

, ?c

k

E

The ? operator applied to k-dimensional geometric objects turns them into

(n −

k)-dimensional geometric objects with the other type of orientation

.

The ? operator applied to k-cochains turns them into

(n − k)-cochains acting on

ge-ometric objects of the other orientation

.

The ? operator is

metric-dependent

and can

therefore not be described in purely topological terms

(33)

The ’Hodge-

?’ operator

The ’Hodge-?’ operator

Remember that ? was the operator which switches between

inner- and outer orientation

.

We can also write down a formal adjoint of this operation

D

?c

k

, c

n−k

E

:=

D

c

n−k

, ?c

k

E

The ? operator applied to k-dimensional geometric objects turns them into

(n −

k)-dimensional geometric objects with the other type of orientation

.

The ? operator applied to k-cochains turns them into

(n − k)-cochains acting on

ge-ometric objects of the other orientation

. The ? operator is

metric-dependent

and can

therefore not be described in purely topological terms

(34)

The ugly stepmother

δ

?

_{= ?δ?}

Recall that

?∂? : C

k

−→ C

k+1 Inner Orientation Outer Orientation ∂ ∂ ∂ * ∂ ∂ ∂ * * * * *

(35)

The ugly stepmother

δ

?

_{and grad, curl and div}

δ

?

_{also represents the grad, curl and div}

δ

?

: C

k

−→ C

k−1 Inner Orientation Outer Orientation grad ∂ curl ∂ div ∂ curl * grad* div *

Note that in contrast to δ, δ

?

_{is a}

_{metric-dependent}

_{version of grad, curl and div and}

can therefore

NOT

be the same as the topological grad, curl and div. We will make

this difference explicit by

grad

∗

, curl

∗

and

div

∗

.

(36)

The ugly stepmother

Laplace-Hodge operator

Inner Orientation Outer Orientation grad ∂ curl ∂ div ∂ curl * grad* div *

The scalar Laplace operator acting on

outward oriented points

is given by

(37)

Laplace-Hodge operator

The vector Laplace operator acting on

outward oriented lines

is given by

[−grad div

∗

+ curl

∗

curl] ~

A

(38)

The ugly stepmother

Laplace-Hodge operator

The vector Laplace operator acting on

outward oriented surfaces

is given by

(39)

Laplace-Hodge operator

The vector Laplace operator acting on

outward oriented volumes

is given by

−div grad

∗

ρ

(40)

The ugly stepmother

Laplace-Hodge operator

(41)

DeR-Introduction Physics & Geometry Integral values Some results

The ugly stepmother

Metric

How do we discretize the

metric-dependent part

?

a

k

, b

k

Ω

:=

Z

Ω

a

k

∧ ?b

k

→

finite element methods

In both cases we need a

continuous reconstruction

,

I

, from our discrete values associated

with geometric objects to continuous functions, i.e.

(42)

The ugly stepmother

Metric

Reduction

Let the

reduction operator

be defined by

R : Λ

k

_{(Ω) −→ C}

k

_(Ω)

D

Ra

(k)

_{, τ}

k

E

:=

Z

τk

a

(k)

Rd = δR

Λ

k

_{−−−−−→ Λ}

d k+1



y

R



y

R

C

k

_{−−−−−→ C}

δ k+1

(43)

Metric

Reduction

Let the

reduction operator

be defined by

R : Λ

k

_{(Ω) −→ C}

k

_(Ω)

D

Ra

(k)

_{, τ}

k

E

:=

Z

τk

a

(k)

Rd = δR

Λ

k

_{−−−−−→ Λ}

d k+1



y

R



y

R

C

k

_{−−−−−→ C}

δ k+1

(44)

The ugly stepmother

Metric

Reconstruction

The

reconstruction operator

needs to satisfy

I : C

k

_{(Ω) −→ Λ}

k h

(Ω) ⊂ Λ

k

_(Ω)

dI = Iδ

and

R ◦ I ≡ I

C

k

_{−−−−−→ C}

δ k+1



y

I



y

I

Λ

k

_{−−−−−→ Λ}

d k+1

Spectral element basis functions which satisfy these relations are called

mimetic spectral

elements

(45)

Metric

Reconstruction

The

reconstruction operator

needs to satisfy

I : C

k

_{(Ω) −→ Λ}

k h

(Ω) ⊂ Λ

k

_(Ω)

dI = Iδ

and

R ◦ I ≡ I

C

k

_{−−−−−→ C}

δ k+1



y

I



y

I

Λ

k

_{−−−−−→ Λ}

d k+1

Spectral element basis functions which satisfy these relations are called

mimetic spectral

elements

(46)

The ugly stepmother

Metric

Discretization ⇔ projection

We define the

projection operator

as

π := I ◦ R

The commutation relations ensure that

dπ

= dIR = IδR = IRd =

πd

Λ

k d

−−−−−→ Λ

k+1



y

π



y

π

Λ

k h d

−−−−−→ Λ

k+1_h

NOTE

: This only holds for the topological grad, curl and div!

NOT

for grad

∗

_{, curl}

∗

_or

div

∗

_.

(47)

The ugly stepmother

Metric

Discretization ⇔ projection

We define the

projection operator

as

π := I ◦ R

The commutation relations ensure that

dπ

= dIR = IδR = IRd =

πd

Λ

k

−−−−−→ Λ

d k+1



y

π



y

π

Λ

k h d

−−−−−→ Λ

k+1_h

(48)

The ugly stepmother

Metric

Discretization ⇔ projection

We define the

projection operator

as

π := I ◦ R

The commutation relations ensure that

dπ

= dIR = IδR = IRd =

πd

Λ

k

−−−−−→ Λ

d k+1



y

π



y

π

Λ

k h d

−−−−−→ Λ

k+1_h

(49)

Mimetic spectral elements

Basis functions 1D

-1 -0.5 0 0.5 1 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 GLL nodal interpolation li (x ) x -1 -0.5 0 0.5 1 -2 -1 0 1 2 3 4 5 x ei (x ) GLL edge interpolation

In 1D we only have

points

and

line segments

, so we use

nodal Lagrange interpolation :

h

i

(x

j

) = δ

ij

Edge interpolation :

Z

xj xj−1

e

i

(x) = δ

i,j

,

e

i

(x) = −

i−1

X

k=0

dh

k

(x)

(50)

The ugly stepmother

Comparison with higher order RT-elements

Stokes problem

10−6 10−4 10−2 100 10−6 10−4 10−2 100 10−3 10−2 10−1 100 10−6 10−4 10−2 100 10−3 10−2 10−1 100 10−6 10−4 10−2 100 ||w-w ||L L 0 2 h ||u-u || H L 1 ||u-u || L L 1 2 h h || (w-w ) ||L L 1 2 h d h ||p-p || L L 2 2 h h 2 3 4 1_ 2 3_ 2 3_ 2 2 2 3 3 3 Raviart−Thomas, N=2 Mimetic Spectral, N=2 Mimetic Spectral, N=3

(51)

Comparison with higher order RT-elements

Stokes problem

10−2 10−1 10−10 10−8 10−6 10−4 10−2 100 10−2 10−1 10−10 10−8 10−6 10−4 10−2 100 10−2 10−1 10−10 10−8 10−6 10−4 10−2 100 N=2 N=4 N=6 N=8 N=2 N=4 N=6 N=8 N=2 N=4 N=6 N=8 ||u-u || H L 1 h ||w-w ||H L 0 h h ||p-p || L L 2 2 h h h nor mal v elocit y - tangen tial v elocit y

tangential velocity - pressure

tangen tial v or ticit y - nor mal v elocit y

tangential vorticity - pressure

2 4 6 8 2 4 6 8 2 4 6 8

(52)

The ugly stepmother

How to avoid grad

∗

, curl

∗

and div

∗

Integration by parts

Finite element methods remove the metric-dependent vector operations through

inte-gration by parts

(da

k

, b

k+1

) = da

k

∧ ?b

k+1

_{= (−1)}

k+1

_a

k

_{∧ d ? b}

k+1

₌

a

k

∧ ?d

∗

b

k+1

= (a

k

, d

∗

b

k+1

)

Vector operations

In conventional vector operations this reads (without boundary)

(53)

Laplace-Hodge operator

The scalar Laplace operator acting on

outward oriented points

is given by

−div

∗

gradφ = f

(54)

The ugly stepmother

Laplace-Hodge operator

The scalar Laplace operator acting on

outward oriented points

is given by

(55)

Laplace-Hodge operator

The scalar Laplace operator acting on

outward oriented points

is given by

(gradφ, gradψ) + b.i. = (f, ψ)

(56)

The ugly stepmother

Laplace-Hodge operator

The vector Laplace operator acting on

outward oriented volumes

is given by

(57)

Laplace-Hodge operator

The vector Laplace operator acting on

outward oriented volumes

is given by

~

q = grad

∗

ρ

div~

q = f

(58)

The ugly stepmother

Laplace-Hodge operator

The vector Laplace operator acting on

outward oriented volumes

is given by

(~

q, ~

p) − (grad

∗

ρ, ~

p) = 0

(59)

Laplace-Hodge operator

The vector Laplace operator acting on

outward oriented volumes

is given by

(~

q, ~

p) + (ρ, div~

p) + b.i. = 0

(60)

The ugly stepmother

Laplace-Hodge operator

(61)

Resonant Cavity problem (benchmark case)

Eigenvalue problem (borrowed from our neighbors)

Maxwell equations with unit coefficients and zero force functions.

∇ ×

∇ × ~

E

= λ ~

E

on

Ω = [0, π]

2

(62)

Introduction Physics & Geometry Integral values Some results Eigenvalue problems Mimetic hp-adaptation Harmonic forms Stokes flow

Resonant Cavity problem (benchmark case)

Results

∆~

u = λ~

u,

div~

u = 0,

Ω = [0, π]

2

Not solvable with standard FEM / SEM,

see [Boffi, Acta Numerica 2010].

8 10 13 17 18 20 25 26 29 32 value of eigenvalue Exact eigenvalues 8 10 13 17 18 20 25 26 29 32 value of eigenvalue Result for N=20, c=0.2

(63)

Resonant Cavity problem (benchmark case)

Results

∆~

u = λ~

u,

div~

u = 0,

Ω = [0, π]

2

Not solvable with standard FEM / SEM,

see [Boffi, Acta Numerica 2010].

10−1 ₁₀0 10−10 10−5 100 h error eigenvalues

h−Convergence of first nine non−zero eigenvalues for N=4 and c=0.2 1 1 2 4 4 5 5 8 9 8

(64)

(65)

Resonant cavity in L-shaped domain

Dirichlet boundary conditions

(66)

(67)

(68)

Eigenfunctions on torus

(69)

(70)

(71)

(72)

(73)

(74)

(75)

Stokes problem

(curl curl

∗

− grad

∗

div) u + grad

∗

p = f

div u = 0

(76)

Stokes problem

ω − curl

∗

u = 0

curl ω + grad

∗

p = f

div u = 0

(77)

Stokes problem

|u| 1.0 0.8 0.6 0.4 0.2 0.0 w 5 4 3 2 1 0 -1 -2 -3 -4 -5 y 0.00 -0.02 -0.04 -0.06 -0.08 -0.10 div u 9.0x10 5.4 1.8 -1.8 -5.4 -9.0 -15 p 5 4 3 2 1 0 -1 -2 -3 -4 -5

(78)

Stokes problem

|u| 1.0 0.8 0.6 0.4 0.2 0.0 w 5 4 3 2 1 0 -1 -2 -3 -4 -5

(79)

Stokes problem

X Y -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 3 -0.5 0 0.5 |u|4.20 3.15 2.10 1.05 0.00 X Y -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 3 -0.5 0 0.5 70 35 0 -35 -70 w X Y -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 3 -0.5 0 0.5 150 75 0 -75 -150 p X Y -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 3 -0.5 0 0.5

(80)

Stokes problem

X Y Z X Y Z DivU 5 x 10 3 1 -14

(81)