HDG METHODS FOR SECOND-ORDER ELLIPTIC PROBLEMS (Numerical Analysis : New Developments for Elucidating Interdisciplinary Problems II)

HDG METHODS FOR SECOND‐ORDER ELLIPTIC PROBLEMS

ISSEI OIKAWA

ABSTRACT. In this_article,wereviewthehybrid(hybridìzedorhybridizable)discon‐

tinuous Galerkin _(HDG) method based on aclassical _hybridfinite element method

for second‐order _{elliptic problems.} Our HDG method was firstlyobtained by sta‐

bilizingthe_{simplified hybrid displacement} method. _Optimal error estimatesin the energy andL^{2} norms were _proved for the Poisson_equation. The method was ex‐

tendedtoconvection‐diffusion_{problems by introducing}akind ofanupwindterm. It

wasverified mathematicallyand numerically that the methodisrobustevenin the

convection‐dominated_case,where the standard finite element method fails duetoits

numerical_instability.

1. INTRODUCTION

In recent years,

hybrid

(hybridized

or

hybridizable)

discontinuous Galerkin

(HDG)

methods have been

_investigated

and

_applied

to various

_problems.

The usual discon‐ tinuous Galerkin

_(DG)

method utilizes two _types of numerical fluxes todealwith the

discontinuity

ofan

approximate

solution _{u_{h}} on inter‐element boundaries. Inthe HDG

method,

a numerical trace

ûh

is introduced to

_approximate

the trace of a solution

besides u_{h}, which isa new unknown and may be called the

hybrid

unknown.

Thenumber of

_degrees

of freedom

_(DOF)

oftheDG methodis much

_larger

than that of the standard finite element method.

_By

the static

_{condensation,}

that _is,

_eliminating

the

_hybrid

unknown

_ûh

_by

u_{h}, weobtainadiscretized

equation

interms of

only

on

ûh.

Asa

result,

the number of DOF oftheHDGmethodcanbe

considerably reduced,

which

isthe main

_advantage

of theHDG methodoverthe DG method. Wenotethat the HDG method has remarkablefeatures besidesthe above

_advantage,

suchas_{superconvergence}

properties andvarious connections withother numerical methods

_{(nonconforming}

and

mixed finite element

_methods,

etc

The HDG method was

firstly

introduced

by

Cockburn et al.

_[10],

in which the

_hy‐

bridization ofthe localdiscontinuous Galerkin

_(LDG)

method

_{(cf. [3])}

is successful to

unify

the formulations of various

_hybrid

methods. An overview of the HDG methods

was

already

provided

in

[10],

andwerefer thereaders toit as a_surveypaper.

Inthis

_article,

werevisitandreviewadifferent

hybridization

of the DG method based on aclassical

hybrid

finite element method.

Hybridization

of the finite element method

was

early

_{proposed by}

Pian in 1964

(cf. [36])

and

_by

de

_Veubeque

in 1965

_{[11, 43].}

_Later,

the

_hybrid

_displacement

methodwas

proposed by Tong

in 1970

[41],

inwhich the

_hybrid

displacement

and

_Lagrange

_multiplier

areintroducedas newunknownsoninter‐element boundaries. The

_simplified

_{hybrid displacement method,}

where the

_{Lagrange multiplier}

is taken to be the normal

_gradient

ofu_{h}, was also

investigated by

Kikuchi and Ando

[19,

21, 22, 20,

23, 18, 24,

25].

Those methods were

partially successful, however, they

suffered from numerical

_instability.

Decades

_later,

in

_[26,

_{30, 35,}

_31],

stabilized methods

were

developed

for linear

elasticity problems

and the Poisson

_equation.

The

instability

was overcome

by introducing

the stabilization

technique

oftheinterior

_penalty

method

[2],

which is describedinSection2. Numerical resultsarenot shown in this

_article,

_see, e.g.

[26,

30,

35,

31].

For

_stationary

convection‐diffusion

_problems,

the HDG method has been

_developed

and there have been _many

_published

_{papers, for}

_example,

see

[27,

13, 7, 8,

32, 14,

6,

38].

Here we focus on the _present author’s work

[32]

because the

resulting

HDG methods listed above are not so different from each other. In the HDG

method,

\mathrm{a} convection‐diffusion _equation is

_decomposed

into diffusive and convective _parts, and

they

are discretized

separately.

The diffusive _partcan be discretized inthe same_way asthePoisson

equation.

The convective_partis discretized

by newly

introducing

akind

ofan

upwind

term. The

_key

idea of

_devising

the

_upwind

term is to switch u_{h} and

ûh

according

tothe outflow and inflow inter‐element boundaries. InSection

_3,

we are

going

to statethe

_upwind

scheme

_proposed

_by

the_presentauthor

_[32]

for convection‐diffusion

problems.

Numerical results willbe

_presented

to validate the

_{stability ,of}

the scheme.

2. HDG METHOD FOR DIFFUSION PROBLEMS

Let

_{$\Omega$\subset \mathbb{R}^{n}(n=2,3)}

be abounded

polygonal

or

polyhedral

domain and

f\in L^{2}( $\Omega$)

be a

_given

function. We consider the Poisson

_equation

with

_{homogeneous boundary}

condition:

(2.1a)

- $\Delta$ u=f

in

_$\Omega$,

(2.1b)

u=0 on \partial $\Omega$.

The HDG methodcan alsobe

applied

tothe

_problems

with

_{non‐homogeneous}

Dirich‐ let and Neumann

_boundary

_conditions,

but we here consider

only

the

homogeneous

Dirichlet

_boundary

condition for

_simplicity.

2.1. Notation. Let

_{\{T_{h}\}_{h>0}}

be a

family

ofmeshes ofthe domain $\Omega$. The

subscript

h

stands for the mesh size h

_{:=\displaystyle \max_{K\in T_{h}}\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{m}(K)}

. We assume that

\{T_{h}\}_{h}

satisfies the

meshesconsistof

_only

_triangles

ortetrahedrons. We alsoassumethat

\{T_{h}\}_{h}

satisfies the local

_{quasi‐uniformity}

_[17],i.e.,

there exists aconstant C such that diam

_{(K)/\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{m}(e)\leq}

C for any K \in

_{T_{h}}

and

_edge

e \subset \partial K.

Throughout

the

article,

the

symbol

Cdenotes

a

generic

constant

_independent

of h. The set of all

edges

of

T_{h}

is denoted

by

\mathcal{E}_{h}

=

\{e\subseteq\partial K : K\in T_{h}\}

. The skeletonof

T_{h}

, defined

by

\displaystyle \bigcup_{K\in T_{h}}\partial K

, is denoted

by

the same

symbol \mathcal{E}_{h}

. Wedefine

L_{D}^{2}(\mathcal{E}_{h})=

{

\hat{v}\in L^{2}(\mathcal{E}_{h})

:\hat{v}=0 on \partial $\Omega$

}.

Wewilluse the standard

notation of the Sobolev spaces

[1],

such as

H^{m}(D)

,

W^{m,p}(D)

,

\Vert \Vert_{m,D}

=

\Vert

.

\Vert_{H^{m}(D)},

|\cdot|_{m,D}=|\cdot|_{H^{m}(D)}

for an

integer

m and adomain D. The

piecewise

orbroken Sobolev

spaces are introduced:

H^{m}(T_{h})

:=

\{v \in L^{2}( $\Omega$) : v|_{K} \in H^{m}(K) \forall K \in T_{h}\}

. The inner

products

aredenoted

by

(u, v)_{T_{h}}

:=\displaystyle \sum_{K\in T_{h}}\int_{K}

uvdx,

_{\{u, v\rangle_{\partial T_{h}}}

:=\displaystyle \sum_{K\in T_{h}}\int_{\partial K}

uvds.

The finite element_spacesfor

_{approximating}

uanditstrace ûare denoted

by W_{h}

and

M_{h},

respectively.

We

impose

the

homogeneous

Dirichlet

boundary

conditionon

M_{h}

,i.e.,

assume

M_{h}\subset L_{D}^{2}(\mathcal{E}_{h})

. In usual cases, the finite element spacesare set tobe piecewise

polynomial

\cdot

spaces ofsame

degree;

W_{h}=P_{k}(T_{h})

and

M_{h}=P_{k}(\mathcal{E}_{h})

.

Recently,

it turned

outthat

_optimal

_convergencescanbe achieved

by setting

W_{h}=P_{k+1}(T_{h})

,

M_{h}=P_{k}(\mathcal{E}_{h})

and

_taking

an L^{2}

‐projection

inthe stabilization_term,see

[28,

_33,

_34,

_39,

_{40, 29, 9,}

37].

2.2. The scheme. We

_give

the formulation ofthe HDG method

_proposed

in

_{[30, 35].}

The methodis

_equivalent

tothe IP‐H method defined in

_[10],

and wewillhere call it so. The IP‐H method is asfollows: find u_{h}\in W_{h} and

ûh

\in

_Mh

such that

(2.2)

a_{h}^{d}(u_{h}

,ûh; v_{h},

\hat{v}_{h})=(f, v_{h})_{ $\Omega$}

\forall v_{h}\in W_{h},

\hat{v}_{h}\in M_{h},

where

a_{h}^{d}(u_{h},\hat{u}_{h|v_{h},\hat{v}_{h})} := (\nabla u_{h}, \nabla v_{h})_{T_{h}}

+

\langle n

.

\nabla u_{h},

\hat{v}_{h}

-v_{h}\rangle_{\partial T_{h}}

+

\langle n

.

\nabla v_{h},

\hat{u}_{h}

-u_{h}\rangle_{\partial T_{h}}

+

_{S_{h}(u_{h},\hat{u}_{h1}v_{h}, \hat{v}_{h})}

,

S_{h}(u_{h},\hat{u}_{h1}v_{h}, \hat{v}_{h}) := \langle $\tau$(\hat{u}_{h} -u_{h}) , \hat{v}_{h} -v_{h}\rangle_{\partial T_{h}}.

Here $\tau$ is astabilization_parameter, which is

_usually

set tobe

_{$\tau$=$\tau$_{0}/h_{e}}

, where $\tau$_{0} is a positive constant, h_{e}

:=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{m}(\mathrm{e})

for an

edge

eand nisthe unit outward normalvector to\partial K. It can be

proved

that the scheme is coercive if$\tau$_{0} is set to be

sufficiently large.

In

_general,

too

_large

$\tau$ is

likely

to

_spoil

the

_{discontinuity}

of the

_approximate

solutions.

Therefore,

we should select a moderate value for $\tau$_{0}. The HDG method coercive for

any

positive

$\tau$_{0} was

already

obtained;

the LDG‐H method

[10]

and the HDG method

using

a

lifting

_operator

[31].

As will be shown

later,

both the methods are

essentially

In thenext _section, we are

_going

todescribe howthe method is derived.

2.3. Derivation. Let K \in

_{T_{h}}

and u \in

H^{2}( $\Omega$)

.

Multiplying

(2.1a)

by

a test function

v\in H^{2}(K)

and

_integrating

it

_by

_parts overK, weget

(\nabla u;\nabla v)_{K}-\{n\cdot\nabla u, v\}_{\partial K}=(f, v)_{K}.

Summing

the above over

K\in T_{h} yields

(2.3)

(\nabla u, \nabla v)_{T_{h}}-\langle n\cdot\nabla u, v\rangle_{\partial T_{h}}=(f, v)_{ $\Omega$}.

Wenow introducea

hybrid

function

\hat{v}\in L_{D}^{2}(\mathcal{E}_{h})

as atest function. Since n\cdot\nabla uand \hat{v}

are both

single‐valued

on

\mathcal{E}_{h}

and\hat{v}vanishes on \partial $\Omega$, thetransmission conditionfollows:

(2.4)

\{n\cdot\nabla u, $\gamma$ v_{\partial T_{h}}=0.

Adding

this into

_(2.3)

and

_symmetrizing

_it,we obtain

(2.5)

(\nabla u, \nabla v)_{T_{h}}+\{n\cdot\nabla u, \hat{v}-v\rangle_{\partial T_{h}}+\langle n\cdot\nabla v, \hat{u}-u\rangle_{\partial T_{h}}=(f, v)_{ $\Omega$},

where û is the trace of u. Since the scheme is still not stable in

general,

we add the

stabilization or

penalty

term $\tau$

{û—u, \hat{v}-v\rangle_{\partial T_{h}}

. Thus weobtain

(2.2).

Remark. In

_(2.2),

_taking

\hat{v}_{h}\equiv 0

on

\mathcal{E}_{h}

and v_{h}\equiv 0on

$\Omega$\backslash K

forsome

K\in T_{h}

, wehave

(\nabla u_{h}, \nabla v_{h})_{K}-\langle n\cdot\nabla u_{h}, v_{h}\rangle_{\partial K}-\langle n

.

\nabla v_{h},

u_{h}\rangle_{\partial K}+\langle $\tau$ u_{h}, v_{h}\rangle_{\partial K}

= (_{f}, v_{h})_{K} + (\hat{u}_{h}, $\tau$ v_{h} - n \nabla v_{h}\rangle_{\partial K},

which

_implies

that

_{u_{h}|_{K}}

canbe determińed

by only

ûh

|\partial

K. There is no direct connec‐ tion between

_{u_{h}|_{K}}

and

_{u_{h}|_{K'}}

for distinct elements

_{K, K'\in T_{h}}

, and

they,

arelinked

only

through

thenumerical trace

_ûh

_{|_{\partial K\cap\partial K'}}

. It enablesustodo theso‐called static condensa‐

tion,

i.e.,

the constructionofalinear_systemintermsof

_{only ûh by}

_{element‐by‐element}

elimination ofu_{h}.

2.4. Error estimates. We _present the outline oferror

analysis

for the IP‐Hmethod. The _energynormis defined

by

\Vert|(v_{h},\hat{v}_{h})\Vert|_{d}^{2}:=\Vert\nabla v_{h}\Vert_{T_{h}}^{2}+\Vert h_{e}^{-1/2}(\hat{v}_{h}-v_{h})\Vert_{\partial T_{h}}^{2},

where

\Vert\nabla v_{h}\Vert_{T_{h}}^{2}:=(\nabla v_{h}, \nabla v_{h})_{T_{h}},

\Vert h_{e}^{-1/2}(\hat{v}_{h}-v_{h})\Vert_{\partial T_{h}}^{2}:=\langle h_{e}^{-1}(\hat{v}_{h}-v_{h})

,

We assume the _{approximation property:} for any v \in

H^{k+1}( $\Omega$)

and its trace \hat{v}, there

exists aconstant C

_independent

of h suchthat

\displaystyle \inf_{v_{h}\in W_{h},\hat{v}_{h}\in M_{h}}\Vert|(v-v_{h}, \hat{v}-\hat{v}_{h})\Vert|_{d}\leq Ch^{k}|v|_{k+1, $\Omega$}.

The

_following

fundamental

_properties

onthe bilinear form

a_{h}^{d}

)

hold. Lemma 1. The

_following

hold.

(1) (Consistency)

Letu bethe exact solution of

_(2.1)

and û denote the traceofu.

Then

a_{h}^{d} (u, \^{u}; v_{h}, \hat{v}_{h})=(f, v_{h})_{ $\Omega$} \forall v_{h}\in W_{h}, \hat{v}_{h}\in M_{h}.

(2) (Boundedness)

There exists aconstant Csuch that

|a_{h}^{d}(w_{h},\hat{w}_{h};v_{h},\hat{v}_{h})|\leq C\Vert|(w_{h},\hat{w}_{h})\Vert|_{d}\Vert|(v_{h}, \hat{v}_{h})\Vert|_{d}

\forall w_{h}, v_{h}\in W_{h},

\hat{w}_{h},

\hat{v}_{h}\in M_{h}.

(3) (Coercivity)

There exists aconstant C such that

a_{h}^{d}(v_{h}, \hat{v}_{h};v_{h},\hat{v}_{h}) \geq C\Vert|(v_{h},\hat{v}_{h})\Vert|_{d}^{2} \forall v_{h}\in W_{h}, \hat{v}_{h}\in M_{h}.

Proof.

Thefull

_proof

was

firstly

given

inthe_presentauthor’s Master’s thesis

[30],

which,

however,

is written in

_Japanese.

_So,

wereferto

[32].

\square

Remark. To obtainanerror

_estimate,

weneednot

only

the

consistency

and boundedness

forW_{h} and.

_{M_{h}}

but also those for

_{H^{2}( $\Omega$)}

.

However,

weomitteditbecause thosearenot

directly

usedinthis article andwehavetointroduce the

_auxiliary

normandnotations. From theabove

_lemma,

the

_{following optimal}

errorestimates canbe deduced. Theorem 2. Assume that

_{\{T_{h}\}_{h}}

satisfies the chunkiness condition and the

_quasi‐

uniformity

and that the

_{approximation}

_property holds. Let u be the solution of

(2.1)

and û denotethetrace ofu. If

u\in H^{k+1}( $\Omega$)

_, then

(2.6)

\Vert|(u-u_{h}, \^{u} - \hat{u}_{h})\Vert|_{d}\leq Ch^{k}|u|_{k+1, $\Omega$}.

By

Aubin‐Nitsche’s

_trick,

we have

2.5.

_Nonsymmetric

schemes. The third term in

_(2.5)

is called a consistent term because

_{\langle n\cdot\nabla v, \hat{u}-u\rangle_{\partial T_{h}}=0}

holds for the exact solutionu andits traceû. Wenote

that the

_nonsymmetric

version ofthe scheme

_{[30, 35]}

can also be

deduced,

for a real number s,

a_{h}^{d}

(u_{h}

,ûh; v_{h},

\hat{v}_{h})=(\nabla u_{h}, \nabla v_{h})_{T_{h}}+\{n\cdot\nabla u_{h}, \hat{v}_{h}-v_{h}\rangle_{\partial T_{h}}

+S\langle n \nabla v_{h}, \hat{u}_{h} -u_{h}\rangle_{\partial T_{h}} +S_{h}(u_{h},\hat{u}_{h};v_{h}, \hat{v}_{h})

.

The IP‐H

_(symmetric)

scheme

_(2.2)

is included as s = 1. We call the scheme with

s\neq 1

the

_nonsymmetric

scheme.

_Although

an

optimal

H^{1}‐errorestimate for the non‐

symmetric

scheme was

proved

as well as the symmetric

scheme,

an

optimal

L^{2}‐error

estimate could not be

_proved

because ofthe lack of the

_{adjoint consistency.}

Notethat the order of convergence in the L^{2} norm is _greaterthanor

equal

tothat of the _energy norm, which follows

easily

from the fact that the energy norm is _stronger than the

L^{2} norm. In

[30, 35],

it was shown

by

numerical

experiments

that the L^{2}‐orders of

convergence are

actually suboptimal

when the

degrees

of

polynomials

are even. For

odd

_polynomials,

the

_optimal

_{convergence in}the L^{2} norm was observed insome cases.

TheL^{2}

_{suboptimality}

ofthe

_nonsymmetric

DGmethodwere

investigated

in

[15,

12,

4],

whereas there are few studies for the HDG method. For the DG

_method,

in

_[15],

the

suboptimal

convergence was in fact demonstrated in the one and two dimensions in

special

caseswhere themeshandexact solutionare

carefully designed.

It

might

be the

case forthe HDG method.

2.6. A

_lifting

operator. In

_[31]

, a

lifting

operator was introduced and the HDG

method

_using

it wasalso

proposed.

Thelocal

lifting

_operator

R_{h}^{\partial K}

:

L^{2}(\partial K)\rightarrow W_{h}(K)^{n}

is defined

_by,

for

_{\hat{ $\mu$}\in L^{2}(\partial K)}

,

(R_{h}^{\partial K}(\hat{ $\mu$}), w)_{K}=\langle\hat{ $\mu$},

w.

n\rangle_{\partial K} \forall w\in W_{h}(K)^{n}.

For

_{$\mu$\in H^{1}(K)}

_, wedefine

R_{h}^{\partial K}( $\mu$)=R_{h}^{\partial K}( $\mu$|_{\partial K})

.

The

_(global)

_lifting

_operator

_{R_{h}}

:

\displaystyle \prod_{K\in T_{h}}L^{2}(\partial K)

\rightarrow

_{W_{h}^{n}}

is defined

_by

_{R_{h}( $\mu$)|_{K}}

=

R_{h}( $\mu$|_{\partial K})

for all

_{K\in T_{h}}

. Note that the

lifting

operator satisfies

(R_{h}( $\mu$), w)_{T_{h}}=\displaystyle \{ $\mu$, w\cdot n)_{\partial T_{h}}=\sum_{K\in T_{h}}\langle $\mu$|_{\partial K}, w\cdot n\rangle_{\partial K} \forall w\in W_{h}^{n}

for

_{$\mu$\displaystyle \in\prod_{K\in T_{h}}L^{2}(\partial K)}

.

TheIP‐H method

_using

the

_lifting

_operator, which is

_going

to be called the IPL‐H method inthis

_article,

read as: find

u_{h}\in W_{h}

and

ûh

\inMh such that

(2.7)

Thanks tothe additional stabilization term

_(Rh(ûh—

u_{h}

), R_{h}(\hat{v}_{h}-v_{h}))_{T_{h}}

, the scheme

is coercive for all $\tau$>0. Wenotethat the scheme

(2.7)

can be rewrittenas

(\nabla u_{h}+R_{h} (\hat{u}_{h} -u_{h}), \nabla v_{h}+R_{h}(\hat{v}_{h}-v_{h}))_{T_{h}}+S_{h}(u_{h},\hat{u}_{h};v_{h},\hat{v}_{h})=(f, v_{h})_{ $\Omega$}.

2.7.

_Equivalence

between the IPL‐H and LDG‐H methods. We are

going

to

show that the IPL‐H method is

_essentially

_equivalent

to the LDG‐H method. In the LDG‐H method

_{(cf. [10]),}

the mixed formulationof

_(2.1)

is considered:

q+\nabla u=0, \nabla\cdot q=f.

Let _{V_{h}} be the finite_{element space} for

_{approximating}

q. TheLDG‐H method reads as:

find

_{q_{h}\in V_{h},}

_{u_{h}\in W_{h}} and

_{\^{u} h\in M_{h}}

such that

(2.8a)

(q_{h}, v)$\tau$_{h}

—

(u_{h}, \nabla v)T_{h}

+

_{\langle\hat{u}_{h},}

v

n)_{\partial T_{h}}

=

0,

\forall v_{\in V_{h},}

(2.8b)

-(q_{h}, \nabla w)_{T_{h}} + \{\hat{q}_{h} n, w\rangle_{\partial T_{h}} = (f, w)_{ $\Omega$}, \forall w \in Wh,

(2.8c)

\langle\hat{q}_{h} n, \hat{w}\rangle_{ $\theta$ T_{h}} = 0, \forall\hat{w} \in M_{h},

where the numerical flux

_{\hat{q}_{h}}

is defined

_by

_q\hat{}

h =

qh + $\tau$

(uh—ûh)n.

Proposition

3. The IPL‐H method is

_equivalent

tothe LDG‐H method with

_{V_{h}=W_{h}^{n}.}

Proof.

To

_{begin with, integrating}

_(2.8a)

_by

_parts,we have

(2.9)

(q_{h}, v)_{T_{h}}+(\nabla u_{h}, v)_{\partial T_{h}}+ \langle

û h—

u_{h},

v\cdot n\rangle_{\partial T_{h}}=0.

Substituting

v=\nabla w into

_(2.9)

_yields

(q_{h}, \nabla w)$\tau$_{h}

+

_{(\nabla u_{h}, \nabla w)_{\partial T_{h}}}

+

\langle\hat{u}_{h}

—

u_{h},n

\nabla w\rangle_{\partial T_{h}}

= 0.

From this and

_(2.8b),

it follows that

(\nabla u_{h}, \nabla w)$\tau$_{h} + \langle\hat{u}_{h} - u_{h}, n \nabla w\rangle_{\partial T_{h}} + \{\hat{q}_{h} n, w\rangle_{\partial T_{h}} = (f, w)_{ $\Omega$}.

By

(2.8c)

and the definition of

_{\hat{q}_{h}}

, we have

(\nabla u_{h}, \nabla w)$\tau$_{h} + \langle\hat{u}_{h} - u_{h}, n. \nabla w\rangle_{\partial T_{h}} + \langle(q_{h} +T (u_{h} -\hat{u}_{h})) n, W -\hat{W}\}_{\partial T_{h}}

(2.10)

= (f, w)_{ $\Omega$} \forall w \in W_{h}, \hat{w} \in M_{h}.

Substituting

v=R_{h}(w-\hat{w})

into

_(2.9),

wededuce

\langle q_{h}\cdot n, w-\hat{w}\rangle_{\partial T_{h}}=(q_{h}, R_{h}(w-\hat{w}))_{T_{h}}

=-(\nabla u_{h}, R_{h}(w-\hat{w}))_{T_{h}}- \langle

ûh—uh,

R_{h}(w-\hat{w})\cdot n\rangle_{\partial T_{h}}

=-(n\cdot\nabla u_{h}, w-\hat{w}\rangle_{\partial T_{h}}+(R_{h} (ûh—uh), R_{h}(\hat{w}-w))_{T_{h}}.

3. HDG METHOD FOR CONVECTION‐DIFFUSION PROBLEMS We consider the

_stationary

convection‐diffusion _equation:

- $\epsilon \Delta$ u+b\cdot\nabla u+cu=f

in

_$\Omega$,

u=0 on

\partial $\Omega$,

where $\epsilon$>0 isthe diffusion

_coefficient,

_{b\in W^{1,\infty}( $\Omega$)^{n}}

is

_{divergence‐free,}

c isa

positive

constant. We are concerned with the convection‐dominated case, i.e. $\epsilon$ \ll

_{\Vert b\Vert_{\infty},}

because the standard finite element method _gets unstable in such a case as is well known.

3.1. The

_upwind

scheme. The convection‐diffusion

_problem

is

_decomposed

into dif‐

fusiveand convective_parts. In theHDG

_method,

_they

are discretized

separately.

The HDG method

_proposed

in

_[32]

isasfollows: find

u_{h}\in W_{h}

and

ûh

\in

_Mh

such that

$\epsilon$ a_{h}^{d}(u_{h},\hat{u}_{h1}v_{h}, \hat{v}_{h}) +a_{h}^{c}(u_{h}, \hat{u}_{h1}v_{h}, \hat{v}_{h}) = (f, v_{h})_{ $\Omega$} \forall v_{h} \in W_{h}, \hat{v}_{h} \in M_{h},

where

a_{h}^{\mathrm{c}}

(u_{h}

,ûh;v_{h},

\hat{v}_{h}):=(b\cdot\nabla u_{h}, v_{h})_{T_{h}}+\{\^{u} h

-u_{h},

(b\cdot n)\hat{v}_{h}\rangle_{\partial T_{h}^{+}}

(3.1)

+

_\{

ûh —

u_{h},

(b\cdot n)v_{h}\}_{\partial T_{h}^{-}}+(cu_{h}, v_{h})_{ $\Omega$}.

The inner

_products

are defined

by

\{w,

_{v\rangle_{\partial T_{h}}\pm}

:=

\displaystyle \sum_{K\in T_{h}} \displaystyle \int_{\partial K^{\pm}}

wvds,

where

\partial K^{-}:=\{x\in\partial K : b(x)\cdot n(x)<0\},

\partial K^{+}:=\{x\in\partial K : b(x)\cdot n(x)\geq 0\},

see also

_Figure

1. We note

_that,

beforethe

_upwind

scheme

_(3.1)

was

proposed

in

[32],

theHDG discretization for the convective_part was

already proposed

in

[10, 13].

The scheme of

_{[10, 13]}

was as follows:

\overline{a}_{h}^{c}(u_{h},\hat{u}_{h1}v_{h}, \hat{v}_{h})

:=

-(u_{h}, b \nabla v_{h})$\tau$_{h}

+

_{\langle(b n)u_{h}}

,vh

-\hat{v}_{h}\rangle_{\partial T_{h}}+

+\{

_{(b. n)}

ûh,

v_{h}-\hat{v}_{h}\rangle_{\partial T_{h}^{-}}+(cu_{h}, v_{h})_{ $\Omega$}.

Wecanshow that both theschemes are

equivalent

to each other.

Proposition

4. It holds that

_{a_{h}^{\mathrm{c}}}

_{(u_{h}}

,

ûh;

v_{h},

\hat{v}_{h}) =\overline{a}_{h}^{c}(u_{h}

,

ûh;

v_{h},\hat{v}_{h})

for allu_{h}, v_{h} \in W_{h}

FIGURE 1. Inflow and outflow element boundaries.

Proof.

We assume that c = 0

only

for

simplicity. Integrating by

parts and

_recalling

that b is

_{divergence‐free,}

we have

(b \nabla u_{h}, v_{h})=-(u_{h}, b . \nabla v_{h})+\{(b . n)u_{h}, v_{h}\rangle_{\partial T_{h}^{+}}+\{(b . n)u_{h}, v_{h}\rangle_{\partial T_{h}^{-}}.

The bilinear form

_{a_{h}^{c}}

canbe rewritten as

a_{h}^{c}

(u_{h}

,

ûh;

v_{h},\hat{v}_{h})=-(u_{h}, b\cdot\nabla v_{h})+

+

=:I_{3} =:I_{4}

andwe have

I_{1}+I_{3}=\langle u_{h}, (b\cdot n)(v_{h}-\hat{v}_{h})\}_{\partial T_{h}^{+}}+\{\hat{u}_{h}, (b\cdot n)\hat{v}_{h}\rangle_{\partial T_{h}^{+}}=:I_{5}+I_{6},

I2+I4 =

\langleûh,

(b\cdot n)v_{h}\}_{\partial T_{h}^{-}}

=_:I7.

Using

the transmission condition

\{

_{(b. n)}

ûh,

\hat{v}_{h}\rangle_{\partial T_{h}}=\langle(b. n)

ûh,

\hat{v}_{h}\rangle_{\partial T_{h}^{+}}+\langle(b. n)

ûh,

\hat{v}_{h}\rangle_{\partial T_{h}^{-}}

=0,

we_get

I_{6}+I_{7}=\{

_{(b. n)}

ûh,

v_{h}-\hat{v}_{h}\rangle_{\partial T_{h}^{-}}.

Thus,

itfollows that

I_{1}+I_{2}+I_{3}+I_{4}=\{(b\cdot n)u_{h},

_{v_{h}-\hat{v}_{h}\rangle_{\partial T_{h}^{+}}+\{(b}

._n

)û

h,

v_{h}-\hat{v}_{h}\}_{\partial T_{h}^{-}},

which

_completes

the

_proof.

3.2. Error

_analysis.

We _quote theoretical results

_proved

in

_[42];

the

_coercivity,

inf‐ sup condition of

\overline{a}_{h}^{c}

(\cdot, a_{h}^{c} )

as

_proved

inthe

_previous

_section)

and statetheerror estimates. Thenorm

corresponding

to theconvectiveterm isdefined

_by

\Vert|(v_{h},\hat{v}_{h})\Vert|_{c}^{2} :=c\Vert u_{h}\Vert_{L^{2}( $\Omega$)}^{2}+\Vert h_{K}^{1/2}b\cdot\nabla v_{h}\Vert_{T_{h}}^{2}+\Vert|b\cdot n|^{1/2}(\hat{v}_{h}-v_{h})\Vert_{\partial T_{h}}^{2}.

Lemma 5. The

_following

hold.

(1) (Coercivity) [42,

Lemma

_4.2]

For all

_{v_{h}\in W_{h}}

and

_{\hat{v}_{h}\in M_{h}}

,wehave the

equality

\displaystyle \overline{a}_{h}^{c}(v_{h}, \hat{v}_{h};v_{h}, \hat{v}_{h})=c\Vert v_{h}\Vert_{L^{2}( $\Omega$)}^{2}+\frac{1}{2}\Vert|b\cdot n|^{1/2}(\hat{v}_{h}-v_{h})\Vert_{\partial T_{h}}^{2}.

(2) (Inf‐sup stability) [42,

Lemma

_4.5]

There exists aconstant Csuch that

C\displaystyle \Vert|(v_{h}, \hat{v}_{h})\Vert|_{c}\leq\sup_{w_{h}\in W_{h},\hat{w}_{h}\in M_{h}} \Vert|(w_{h},\hat{w}_{h})\Vert|_{c}

\overline{a}_{h}^{c}(v_{h}.., \hat{v}_{h};w_{h},\hat{w}_{h}) \forall v_{h}\in W_{h}, \hat{v}_{h}\in M_{h}.

From this and Lemma

_1,

the

_following

error estimate can be

obtained,

which is an

improved

result shown in

_[32,

Theorem

_3].

Theorem 6. If

_{u\in H^{k+1}( $\Omega$)}

, then wehave

$\epsilon$^{1/2}(

u—uh,

u-\hat{u}_{h})\Vert|_{d}+\Vert|(u-u_{h}, u-\hat{u}_{h})\Vert|_{c}\leq C($\epsilon$^{1/2}+h^{1/2})h^{k}\Vert u\Vert_{k+1}.

In

_particular,

if $\epsilon$ is smaller than h and

\{T_{h}\}_{h}

is

quasi‐uniform,

then the errors in the

streamline and L^{2} norms arebounded as

\Vert b\cdot\nabla(u-u_{h})\Vert_{T_{h}} \leq Ch^{k}\Vert u\Vert_{k+1, $\Omega$},

\Vert u-u_{h}\Vert_{L^{2}( $\Omega$)}\leq Ch^{k+1/2}\Vert u\Vert_{k+1, $\Omega$}.

Remark. The above error estimates are no

longer

useful when $\epsilon$ is very small since

\Vert u\Vert_{k+1, $\Omega$}

may

depend

on

negative

_powers of $\epsilon$. In

[32],

to show the robustness of the

HDG

_method,

itwasshownthat the HDGsolution isclosetothesolution of the reduced

problem.

3.3. Numericalresults. To validate the

_stability

ofthe HDGmethod in the convection‐ dominated_case, we

provide

numericalresults and _comparethe HDG method with the standard finite element method and the streamline

_upwind

Petrov‐Galerkin

_(SUPG)

method. Thetest

_problem

is as follows.

- $\epsilon \Delta$ u+(1,0)^{T}\cdot\nabla u=1

in $\Omega$

_:=(0,1)^{2},

u=0 on \partial $\Omega$.

When $\epsilon$is very

small,

inthis

problem,

boundary layers

_appear

along

thecharacteristic

boundary

(near

y=0

and

_y=1

₎

and the outflow

_boundary

_(near

x=1

_),

which causes

We fixed the mesh size h \approx

1/10

_{and the stabilization parameter $\tau$_{0}} = 8

and,

an unstructured

_triangular

mesh was used. All _computations were carried out with

$\Gamma$ \mathrm{r}\mathrm{e}\mathrm{e} $\Gamma$ \mathrm{e}\mathrm{m}++

_[16]

. In

Figure

2, the numerical solutions are

displayed

for $\epsilon$ =

10^{-k}(k

=

1,3, 5,

7).

The standard finite elementsolutionsbreak downduetonumerical

instability

when

_{$\epsilon$\leq 10^{-3}}

. The SUPG method seemstobe robust even for

very small $\epsilon$,

however,

overshoot appearsnearthe outflow

boundary.

Weobservethat the HDGmethod isro‐ bustandfree fromthe overshoot

_phenomena

unlike theSUPG method. Wecan alsosee that the HDG solutions_getcloser to the solution ofthe reduced

_problem,

_{u(x, y)=x,}

as $\epsilon$ tends tozero.

REFERENCES

[1]

R. A. ADAMS ANDJ. J. F. _FOURNIER, Sobolevspaces (2nded AcademicPress, Amsterdam,

2003.

[2]

D. N. _ARNOLD, An interior_{penalty finite}element method with discontinuouselements,SIAM J.

Numer,Anal., 19_(1982),pp. 742‐760.

[3]

D. N._ARNOLD, $\Gamma$.BREZZI, B.COCKUURN,ANDL. D.MARINI, Unified analysis ofdiscontinuous

Galerkin methods_{for elliptic problems,} SIAM J. Numer._Anal., 39 _(2001),pp.1749‐1779.

[4] B. AYUSODE_Dios,F. _BREZZI, O._HAVLE, AND L. D. _MARINI, L^{2}‐estimates_forthe DG IIPG‐0

scheme,Numer. Meth. Partial Differ._Equ.,28 _(2012),pp. 1440‐1465.

[5]

S. C. BRENNERANDL. R. _SCOTT, The mathematical_{theory of finite}element methods

_(3rd

ed

Springer,New_York, 2008.

[6]

H._CHEN,J.LI,ANDW.QIU,RobustaposterzorierrorestimatesforHDG methodforconvection‐ diffusion equations, IMA J. Numer._Anal., 36 _(2016),pp.437‐462.

[7]

Y. CHENAND B. _{COCKBURN, Analysis of variable‐degree}HDG methods_{for convection‐diffusion}

equations. Part I: General_{nonconforming meshes,} IMA J. Numer. _Anal., 32 _{(2012), pp.} 1267‐ 1293.

[8]

—,Analysis of variable‐degreeHDG methodsfor convection‐diffusionequations. Part II: Semi‐

matching nonconforming meshes, Math. _Comp.,83 _(2014),pp.87−111.

[9]

B. _COCKBURN,Static_{condensation, hybridization,} andthe_devisingofthe HDG_methods,Lecture

Notesin_{Computational}Science andEngineering, 114 _(2016),pp. 129‐177.

[10]

B. _COCKBURN, J. _{GOPALAKRISHNAN,} AND R. _{LAZAROV, Unified hybridization of}discontinuous

Galerkin, mixed_{J} and continuous Galerkin methods_forsecond order _{elliptic problems,} SIAM J.

Numer._Anal., 47 _(2009),pp. 1319‐1365.

[11]

B. F. DE _{VEUBEKE, Displacement}and _equilibrium models in the_finite element method, O. C.

Zienkiewicz and G. S. Holister _(eds) Streess_Analysis, John_{Wiley, (1965),}pp. 145‐197.

[12] V.DOLEJŠí AND O. HAVLE, TheL^{2}‐optimality ofthe IIPG methodforodd_{degrees of polynomial}

approximationin ID,J.Sci. Comput.,42 _(2010),pp. 122‐143.

[13]

H. EGGER AND J. _SCHöBERL, A _hybridmixed discontinuous_{galerkin finite‐element}method_for

convection‐diffusion problems,IMA J. Numer._Anal., 30 _(2010),pp. 1206‐1234.

[14]

G. _FU, W. _Qiu, AND W. _ZHANG,An _{analysis of}HDG methods_forconvection‐dominated_diffu‐ sion_problems,ESAIM:M2AN,49 _(2015),pp.225‐256.

[15]

J. GUZMÁNANDB. _{RIVIÈRE,Sub‐optimalconvergenceof non‐symmetric}discontinuous Galerkin

methods_forodd_{polynomial approximations,}J. Sci. _Comput.,40_(2009), pp.273‐280.

[16]

$\Gamma$. HECHT, Newdevelopmentin FreeFem++

,J. Numer. Math.,20(2012),pp. 251‐265.

[17]

F. _{KIKUCHI, Rellich‐type} discrete _{compactnessfor}some discontinuous Galerkin _{FEM, Japan} J.

Indust. _{Appl. Math.,}29 _(2012),pp. 269‐288.

[18]

$\Gamma$. KIKUCHIANDY. ANDO, Convergence of simplified hybrid displacementmethodfor plate bend‐

ing,J. Fac. _Eng. Univ._{Tokyo (B),}31 _(1972),pp.693‐713.

[19]

—, Finite elementanalysis ofvibrationofthin elasticplates by simplified hybrid displacement

method, J.Nucl. Sci. _Technol.,9_(1972),pp. 189‐190.

[20] —,Lumped finite element basesforbeams andplates,J. Nucl. Sci.Technol.,9(1972),pp.749‐ 751.

[21] —,A newvariationalfunctional forthefinite‐elementmethod andits _applicationto_plate and

shell_problems,Nucl. _{Eng. Des.,} 21 _(1972),pp.95‐113.

[22]

—, Simplified hybrid displacement method applied to plate buckling problems, J. Nucl. Sci.

Technol., 9 _(1972),pp. 497‐499.

[23]

—, Somefinite element solutionsfor plate bending problems by simplified hybrid displacement

method, Nucl._{Eng. Des.,} 23_(1972), pp. 155‐178.

[24]

—, Application of simplified hybrid displacementmethodtolarge deflection analysis ofelastic‐

plastic platesand_shells, J. Fac. _Eng.Univ._{Tokyo (B),}32 _(1973),pp. 117‐135.

[25] —, Simplifiedhybrid displacementmethodforlinearfinite‐element analysis of general shells,

Int. J. Pres. Ves._Pip.,2 _(1974),pp. 155‐164.

[26]

F. _KIKUCHI,K. _ISHII, ANDI._OIKAWA, DiscontinuousGalerkinFEMof hybrid displacementtype

‐

development of polygonal elements−, Theor._Appl.Mech. _Japan,57_(2009),pp. 395‐404.

[27]

R. J. LABEURAND G. N. WELLS, A Galerkin_interfacestabilisation method_forthe advection‐

diffusionand_{incompressible}Navier‐Stokes equations,Comput.Methods_Appl.Mech._Engrg., 196

(2007), pp.4985‐5000.

[28] C. _{LEHRENFELD, Hybrid} discontinuous Galerkin methods _{for incompressible flow problems,}

diploma thesis, MathCCES_/IGPM,RWTH_{Aachen, (2010).}

[29]

C. LEHRENFELD AND J. _{SCHöBERL, High} order _{exactly divergence‐free hybrid} discontinuous

Galerkin methods_{for unsteady incompressible flows, Comput.} Methods_Appl.Mech._Engrg., 307

(2016),pp.339‐361.

[30]

I. _OIKAWA, A discontinuous Galerkin_finite element method _{of hybridtype (in Japanese,} trans‐

lated_title),Master’s thesis: Graduate School of MathematicaìSciences,The_Universityof_Tokyo,

(2008).

[31]

—, Hybr $\iota$ dized discontinuous Galerkin method with lifting operator, JSIAM Lett., 2 (2010),

pp. 99‐102.

[32] —, Hybridizeddiscontinuous Galerkin methodfor convection‐diffusion problems, JapanJ. In‐

dust. Appl. Math.,31 _(2014),pp. 335‐354.

[33]

—,A hybridizeddzscontinuous Galerkin method with reducedstabilization,J.Sci. Comput., 65

(2015),pp. 327‐340.

[34]

—, Analysis ofa reduced‐order HDG methodforthe Stokes _equations, J. Sci. _Comput., 67

[35] I.OIKAWAANDF. _KIKUCHI,Discontinuous Galerkin FEMof hybndtype,JSIAMLett.,2 _(2010),

pp.49‐52.

[36]

T. H. H. PIANANDC.‐C. _{WU, Hybrid}and_{incompatible finite}element_methods,vol.4,Chapman

& Hall_{/\mathrm{C}\mathrm{R}\mathrm{C}}, 2006.

[37]

W. QIU, J. _SHEN, AND K. _SHI, An HDG methodfor linear _elasticity with _{strong symmetric}

stresses,to appear.

[38] W. _QIU AND K. _SHI, An HDG method_for convection _{diffusion equation,} J. Sci. Comput., 66

(2016),pp.346‐357.

[39]

—,An HDG methodforconvectiondiffusionequation,J.Sci.Comput.,66(2016),pp.346‐357.

[40]

—,A superconvergentHDG methodfortheincompressibleNavier‐Stokes equationsongeneral polyhedral meshes,IMA J. Numer._Anal.,36 _(2016),pp. 1943‐1967.

[41] P. TONG,New_{displacement hybrid finite}element models_forsolid_continua,Int. J. Numer. Meth.

Eng.,2 _(1970),pp.73‐83.

[42]

G. N. _{WELLS,Analysis of}aninterfacestabilizedfinite element method: the _{advection‐diffusion‐}

reaction_equation, SIAM J. Numer.Anal.,49 _(2011),pp.87‐109.

[43] O. C. ZIENKIEWIC\mathrm{Z}^{\cdot}, Displacement and equilibrium models in thefinite element method byB.

Fraeijs de Veubeke, Chapter 9, Pages 145‐197ofStressAnalysis, Edited _by O. C. Zienkiewicz

andG. S. Holister,Published_byJohn _WileyớSons, 1965,Int. J. Numer. Meth.Eng.,52_(2001),

pp. 287‐342.

GRADUATE SCHOOLOFFUNDAMENTAL SCIENCE AND_ENGINEERING, FACULTYOFSCIENCE AND

ENGINEERING, WASEDA _UNIVERSITY, 3‐4‐1 _{OKUBO, SHINJUKU‐KU, TOKYO, 169‐8555,} JAPAN

FIGURE 2. Solutions of the HDG

_(left),

SUPG

_(center)

iand

standard finite element

_(right)

_methods,

where

_{$\epsilon$=10^{-k}(k=1,3,5,7)}

from_top to bottom.