Discontinuous Galerkin Methods for an
Elliptic Variational Inequality of 4th-Order
Fei Wang,1
Weimin Han,2
Jianguo Huang3
and
Tianyi Zhang4
Abstract. Discontinuous Galerkin (DG) methods are studied for solving an
elliptic variational inequality of 4th-order. Numerous discontinuous Galerkin
schemes for the Kirchhoff plate bending problem are extended to the variational
inequality. Numerical results are presented to illustrate convergence orders of the
different methods.
Keywords. Variational inequality of 4th-order, discontinuous Galerkin method,
Kirchhoff plates
AMS Classification. 65N30, 49J40
1
Introduction
In this paper, we introduce and study several discontinuous Galerkin (DG) methods for
solving an elliptic variational inequality of 4th-order.
1.1
Discontinuous Galerkin methods
Discontinuous Galerkin methods form an important family of nonconforming finite element
methods for boundary value or initial-boundary value problems of hyperbolic, parabolic
and elliptic partial differential equations. We refer to [10] for a historical account about
DG methods. Discontinuous Galerkin methods use piecewise smooth functions without
much global smoothness to approximate the problem solution, and connect the information
between two neighboring elements by the so-called numerical traces. The practical interest
in DG methods is due to their flexibility in mesh design and adaptivity, in that they allow
elements of arbitrary shapes, irregular meshes with hanging nodes, and the discretionary local
1
Department of Mathematics, Pennsylvania State University, State College, PA 16802, USA.
Department of Mathematics & Program in Applied Mathematical and Computational Sciences, University of Iowa, Iowa City, Iowa 52242, USA
3
Department of Mathematics, and MOE–LSC, Shanghai Jiao Tong University, Shanghai 200240, China;
Division of Computational Science, E-Institute of Shanghai Universities, Shanghai Normal University, China
4
Department of Mathematics, University of Iowa, Iowa City, Iowa 52242, USA
2
1
shape function spaces. In addition, the increase of the locality in discretization enhances the
degree of parallelizability.
There are two basic approaches to construct DG methods for linear elliptic boundary
value problems. The first approach is to choose an appropriate bilinear form that contains
penalty terms to penalize jumps across neighboring elements to make the scheme stable. The
second approach is to choose appropriate numerical fluxes to make the method consistent,
conservative and stable. In [1] and [2], Arnold, Brezzi, Cockburn, and Marini provided a
unified error analysis of DG methods for linear elliptic boundary value problems of 2nd-order
and succeeded in building a bridge between these two families of DG methods, establishing
a framework to understand their properties, differences and the connections between them.
In [21], numerous DG methods were extended for solving elliptic variational inequalities of
2nd-order, and a priori error estimates were established, which are of optimal order for linear
elements. DG methods for the Signorini problem and a quasistatic contact problem were
also studied in [22] and [23], respectively. In this paper, we study DG methods to solve an
elliptic variational inequality of 4th-order for the Kirchhoff plates. The novel difficulty in
constructing stable DG methods for such problems is caused by their high order of four.
The major known DG methods for the biharmonic equation in the literature are primal
DG methods, namely variations of interior penalty (IP) methods ([4, 5, 7, 12, 16, 17, 18,
20]). Fully discontinuous IP methods, which cover meshes with hanging nodes and locally
varying polynomial degrees, ideally suited for hp-adaptivity, were investigated systematically
in [16, 17, 18, 20] for biharmonic problems. In [12], a C 0 IP formulation was introduced
for Kirchhoff plates and quasi-optimal error estimates were obtained for smooth solutions.
Unlike fully discontinuous Galerkin methods, C 0 type DG methods do not “double” the
degrees of freedom at element boundaries. To consider the C 0 IP method under a weak
regularity assumption on the solution, a rigorous error analysis was presented in [7]. A
weakness of this method is that the penalty parameter can not be precisely quantified a
priori, and the penalty parameter must be chosen suitably large to guarantee stability.
However, a large penalty parameter has a negative impact on accuracy. Based on this
observation, a C 0 DG (CDG) method was introduced in [24], where the stability condition
can be precisely quantified. In [15], a consistent and stable CDG method, called the LCDG
method, was derived for the Kirchhoff plate bending problem, which can be viewed as an
extension of the LDG method studied in [8, 9]. We will extend these three methods and
propose two other new CDG methods to solve the elliptic variational inequality of 4th-order.
1.2
Kirchhoff plate bending problem
We now describe a Kirchhoff plate bending problem. Let Ω ⊂ R2 be a bounded polygonal
domain with boundary Γ. The boundary value problem of a clamped Kirchhoff plate under
2
a given scaled vertical load f ∈ L2 (Ω) is (cf. [19])
( P
2
i,j=1 Mij,ij (u) + f = 0 in Ω,
u = ∂ν u = 0 on Γ,
where
Mij (u) := −(1 − κ)∂ij u − κ
2
X
∂kk uδij ,
(1.1)
1 ≤ i, j ≤ 2,
k=1
δij is the usual Kronecker delta, κ ∈ (0, 0.5) denotes the Poisson ratio of an elastic thin plate
occupying the region Ω and ν stands for the unit outward normal vector on Γ. As in [15],
we introduce an auxiliary matrix-valued function by
σ := −(1 − κ)∇2 u − κ tr(∇2 u)I,
(1.2)
where I is the identity matrix of order 2 and tr(·) is the trace operation on matrices. Here,
we denote the gradient of v by ∇v and the Hessian of v by ∇2 v, i.e.,
∂11 v ∂12 v
2
t
∇ v := ∇(∇v) = ∇((∂1 v, ∂2 v) ) =
.
∂21 v ∂22 v
Then, the problem (1.1) can be rewritten as

κ
1

2


 1 − κ σ − 1 − κ2 (trσ)I = −∇ u in Ω,
− ∇ · (∇ · σ) = f in Ω,



 u = ∂ u = 0 on Γ.
(1.3)
ν
For a vector-valued function v = (v1 , v2 )t and a matrix-valued function σ = (σij )2×2 , we
define their divergences by
∇ · v := v1,1 + v2,2 ,
∇ · σ := (σ11,1 + σ21,2 , σ12,1 + σ22,2 )t .
We denote the normal and tangential components of a vector v on the boundary by vν = v ·ν
and v τ = v − vν ν. Similarly, for a tensor σ, we define its normal component σν = σν · ν
and tangential component σ τ = σν − σν ν. We have the decomposition formula
(σν) · v = (σν ν + σ τ ) · (vν ν + v τ ) = σν vν + σ τ · v τ .
For two matrices τ and σ, their double dot inner product and corresponding norm are
P
σ : τ = 2i,j=1 σij τij and |τ | = (τ : τ )1/2 .
The following result is very useful for the analysis of DG methods, which can be verified
directly by integration by part.
3
Lemma 1.1 Let D be a bounded domain with a Lipschitz boundary ∂D. For a symmetric
matrix-valued function τ and a scalar function v, the following two identities hold
Z
Z
Z
Z
2
v ∇ · (∇ · τ ) dx =
∇ v : τ dx −
∇v · (τ n) ds +
v n · (∇ · τ ) ds,
D Z
DZ
∂D
∂D
Z
2
∇ v : τ dx = −
∇v · (∇ · τ ) dx +
∇v · (τ n) ds,
D
D
∂D
whenever the terms appearing on both sides of the above identities make sense. Here n is
the unit outward normal to ∂D.
Multiplying the second equation in (1.3) by a test function v ∈ H02 (Ω) and noticing
v = ∂ν v = 0, we get the following equation by applying Lemma 1.1,
Z
Z
2
− σ : ∇ v dx =
f v dx.
(1.4)
Ω
Ω
With the definition of σ, the weak formulation of the problem (1.3) can be derived from
(1.4) as follows:
Find u ∈ H02 (Ω) : a(u, v) = (f, v) ∀ v ∈ H02 (Ω),
(1.5)
where the bilinear form is
Z
a(u, v) =
∆u ∆v + (1 − κ) (2 ∂12 u ∂12 v − ∂11 u ∂22 v − ∂22 u ∂11 v) dx,
(1.6)
Ω
and the linear form is
Z
(f, v) =
f v dx.
Ω
In this paper, we consider an elliptic variational inequality (EVI) of the 4th-order for
Kirchhoff plates ([11]):
Find u ∈ K :
a(u, v − u) ≥ (f, v − u) ∀ v ∈ K.
(1.7)
Here,
K = v ∈ H 2 (Ω) ∩ H01 (Ω) : ∂ν v ≥ 0 on Γ .
(1.8)
Applying the standard theory on elliptic variational inequalities (e.g., [3, 13]), we know the
problem (1.7) has a unique solution u ∈ K. This variational inequality describes a simply
supported plate. The displacement u is held fixed on the boundary, and the rotation of the
plate is unilateral on the boundary.
In error analysis of numerical solutions for the problem (1.7), we need to take advantage
of pointwise relations satisfied by the solution u.
4
Proposition 1.2 Assume the solution of the problem (1.7) has the regularity u ∈ H 3 (Ω).
Then,
−∇ · (∇ · σ) = f
σ τ = 0, σν ≤ 0, ∂ν u ≥ 0, σν ∂ν u = 0
a.e. in Ω,
a.e. on Γ. (1.9)
Proof. Note that σ is defined by (1.2). Then σ ∈ H 1 (Ω)2×2 . We rewrite (1.7) as
Z
−σ : ∇2 (v − u) − f (v − u) dx ≥ 0 ∀ v ∈ K.
Ω
Take v = u ± ϕ for any ϕ ∈ C0∞ (Ω) to obtain
Z
Z
2
− σ : ∇ ϕ dx =
f ϕ dx ∀ ϕ ∈ C0∞ (Ω).
Ω
Ω
Thus,
−∇ · (∇ · σ) = f
in the sense of distribution.
Since f ∈ L2 (Ω), we deduce that ∇ · (∇ · σ) ∈ L2 (Ω) and
−∇ · (∇ · σ) = f
a.e. in Ω.
This is the first relation in (1.9).
Since ∇ · σ ∈ L2 (Ω)2 and ∇ · (∇ · σ) ∈ L2 (Ω), we have
Z
Z
− ∇ · (∇ · σ) v dx = (∇ · σ) · ∇v dx ∀ v ∈ C0∞ (Ω).
Ω
Ω
From this relation, we obtain
Z
Z
− ∇ · (∇ · σ) v dx = (∇ · σ) · ∇v dx ∀ v ∈ H01 (Ω).
Ω
Ω
Therefore, for any v ∈ H01 (Ω) ∩ H 2 (Ω),
Z
Z
− ∇ · (∇ · σ) v dx = (∇ · σ) · ∇v dx
Ω
ZΩ
Z
= (σν) · ∇v ds − σ : ∇2 v dx,
Γ
Ω
i.e.,
Z
Z
f v dx −
a(u, v) =
Ω
(σν) · ∇v ds ∀ v ∈ H01 (Ω) ∩ H 2 (Ω).
Γ
Recalling the inequality (1.7), we then have
Z
− ∇(v − u) · (σν) ds ≥ 0 ∀ v ∈ K.
Γ
5
(1.10)
In (1.10), we choose v = 0 and 2u in turn to obtain
Z
∇u · (σν) ds = 0.
(1.11)
Γ
Hence,
Z
∇v · (σν) ds ≤ 0 ∀ v ∈ K.
(1.12)
Γ
By the arbitrariness of v ∈ K, it follows that σ τ = 0. Then we get
Z
σν ∂ν v ds ≤ 0 ∀ v ∈ K.
(1.13)
Γ
By the arbitrariness of ∂ν v ≥ 0 on Γ for v in K, we have σν ≤ 0 a.e. on Γ. Back to (1.11),
we further deduce σν ∂ν u = 0 a.e. on Γ. So the relations on the boundary Γ in (1.9) hold.
Throughout the paper, we assume the solution of the problem (1.7) has the regularity
u ∈ H 3 (Ω). In [14, pp. 323–327], one can find regularity results u ∈ H 3 (Ω) for solutions of
some variational inequalities of 4th-order.
The rest of the paper is as follows. In Section 2, we introduce some notations and
derive some C 0 discontinuous Galerkin methods for solving the Kirchhoff plate bending
problem, and extend them to solve the elliptic variational inequality of 4th-order. In Section
3, consistency of the CDG methods, boundedness and stability of the bilinear forms are
presented. In the final section, we report results from a numerical example.
2
2.1
DG methods for Kirchhoff plate problem
Notations
We introduce some notations frequently used later on. For a given function space B, let
(B)2×2
:= {τ ∈ (B)2×2 : τ t = τ }. Given a bounded set D ⊂ R2 and a positive integer m,
s
H m (D) is the usual Sobolev space with the corresponding norm k · km,D and semi-norm
| · |m,D , which are abbreviated by k · km and | · |m , respectively, when D is chosen as Ω. k · kD
is the norm of the Lebesgue space L2 (D). We assume Ω is a polygonal domain and denote
by {Th }h a family of triangulations of Ω, with the minimal angle condition satisfied. Let
hT = diam(T ) and h = max{hT : T ∈ Th }. For a triangulation Th , let Eh be the union of
all edges. We have the non-overlapping decomposition Eh = Ehi ∪ Ehb , where Ehi ⊂ Eh is the
union of all interior edges, i.e., the union of all edges in Eh that do not lie on Γ, and similarly,
Ehb ⊂ Eh is the union of the edges on Γ. For any e ∈ Eh , denote by he its length. Related to
6
the triangulation Th , let
n
o
2×2
Σ := τ ∈ L2 (Ω) s : τij |T ∈ H 1 (T ) ∀ T ∈ Th , i, j = 1, 2 ,
V := v ∈ H01 (Ω) : v|T ∈ H 2 (T ) ∀ T ∈ Th .
The corresponding finite element spaces are
n
o
2×2
Σh := τ h ∈ L2 (Ω) s : τhij |T ∈ Pl (T ) ∀ T ∈ Th , i, j = 1, 2 ,
Vh := vh ∈ H01 (Ω) : vh |T ∈ P2 (T ) ∀ T ∈ Th .
Here, for a triangle T ∈ Th , Pl (T ) and P2 (T ) are the polynomial spaces of degrees l and 2,
respectively, with l = 0, 1. Note that we have the following property
∇2h Vh ⊂ Σh ,
κ
1
Σh −
(trΣh ) I ⊂ Σh ,
1−κ
1 − κ2
(2.1)
where ∇2h Vh |T := ∇2 (Vh |T ) for any T ∈ Th . ∆h v is defined by the relation ∆h v = ∆v on
any element T ∈ Th . Considering the CDG methods for the variational inequality (1.7), we
introduce the corresponding finite element set
Kh := {vh ∈ Vh : ∂ν vh ≥ 0 at all vertex nodes on Γ} .
On each element T , vh |T is a quadratic polynomial function, so ∂ν vh is a linear polynomial
on each edge. With the constraint that ∂ν vh ≥ 0 at all vertex nodes on Γ, we know
∂ν vh ≥ 0 on Γ.
(2.2)
For a function v ∈ L2 (Ω) with v|T ∈ H m (T ) for all T ∈ Th , define broken norm and
seminorm by
X
1/2
X
1/2
2
2
kvkm,h =
kvkm,T
, |v|m,h =
|v|m,T
.
T ∈Th
T ∈Th
The above symbols are used in a similar manner when v is a vector or matrix-valued function.
Throughout this paper, C denotes a generic positive constant independent of h and other
parameters, which may take different values at different occurrences. To avoid writing these
constants repeatedly, we use “x . y” to mean that “x ≤ Cy”. For two vectors u and v,
u ⊗ v is a matrix with ui vj as its (i, j)-th component.
Consider two elements T + and T − with a common edge e ∈ Ehi and let n+ and n− be
their outward unit normals on e. For a scalar-valued function v, set its restriction on T ±
7
by v ± = v|T ± . Similarly, for a matrix-valued function τ , write τ ± = τ |T ± . Then define
averages and jumps on e ∈ Ehi as follows:
1
{v} = (v + + v − ), [v] = v + n+ + v − n− ,
2
1
{∇v} = (∇v + + ∇v − ), [∇v] = ∇v + · n+ + ∇v − · n− ,
2
1
{τ } = (τ + + τ − ), [τ ] = τ + n+ + τ − n− .
2
For e ∈ Ehb , the above definitions are modified:
{v} = v, [v] = vν,
{∇v} = ∇v, [∇v] = ∇v · ν,
{τ } = τ , [τ ] = τ ν.
The jump J·K of the vector ∇v is
1
J∇vK = (∇v + ⊗ n+ + n+ ⊗ ∇v + + ∇v − ⊗ n− + n− ⊗ ∇v − ) on e ∈ Ehi ,
2
1
J∇vK = (∇v ⊗ ν + ν ⊗ ∇v) on e ∈ Ehb .
2
2×2
Define a global lifting operator r i : (L2 (Ehi ))s → Σh by
Z
Z
2×2
φ : {τ } ds ∀ τ ∈ Σh , φ ∈ L2 (Ehi ) s .
r i (φ) : τ dx = −
Ω
(2.3)
Ehi
2×2
Moreover, for each e ∈ Eh , introduce a local lifting operator r e : (L2 (e))s → Σh by
Z
Z
2×2
r e (φ) : τ dx = − φ : {τ } ds ∀ τ ∈ Σh , φ ∈ L2 (e) s .
(2.4)
e
Ω
It is easy to check that the following identity holds
X
2×2
r i (φ) =
r e (φ|e ) ∀ φ ∈ L2 (Ehi ) s ,
e∈Ehi
so we have
kr i (φ)k2 = k
X
r e (φ|e )k2 ≤ 3
e∈Ehi
X
e∈Ehi
8
kr e (φ|e )k2 .
(2.5)
2.2
Discontinuous Galerkin formulations
We first present the derivation of a general primal formulation of CDG methods for the
problem (1.5). By the first equation in (1.3), the first relation in (1.9) and Lemma 1.1, we
have
Z Z
Z
κ
1
σ:τ−
trσ trτ dx =
∇u · (∇ · τ ) dx −
∇u · (τ nT ) ds
1−κ
1 − κ2
T
T
∂T
for any smooth second-order tensor-valued function τ , and
Z
Z
Z
Z
2
− f v dx =
∇ v : σ dx −
∇v · (σnT ) ds +
T
T
∂T
nT · (∇ · σ)v ds
∂T
for any smooth scalar-valued function v. Thus, consider a CDG approximate solution
(σ h , uh ) ∈ Σh × Vh governed by
Z
Z
Z κ
1
c h · (τ h nT ) ds,
σh : τ h −
trσ h trτ h dx =
∇uh · (∇ · τ h ) dx −
∇u
2
1
−
κ
1
−
κ
T
∂T
T
(2.6)
Z
Z
Z
∇2h vh : σ h dx −
∇vh · (b
σ h nT ) ds, (2.7)
− f vh dx =
T
T
∂T
[
for all (τ h , vh ) ∈ Σh × Vh and all T ∈ Th . Here we take ∇
· σ h = 0 in the last equation for
sake of simplicity. To derive numerous CDG methods, we first introduce an identity. For a
scalar function v and a symmetric matrix-valued function τ , smooth on each element of the
partition Th , after a direct manipulation, we have
XZ
XZ
XZ
∇v · (τ nT ) ds =
[τ ] · {∇v} ds +
{τ } : J∇vK ds.
(2.8)
T ∈Th
∂T
e∈Ehi
e
e∈Eh
e
We now sum the equations (2.6) and (2.7) over all T ∈ Th and apply Lemma 1.1 and
(2.8) to obtain
Z Z
Z
1
κ
2
c h } · [τ h ] ds
{∇uh − ∇u
σh : τ h −
trσ h trτ h dx = − ∇h uh : τ h dx +
2
i
1
−
κ
1
−
κ
Ω
Ω
Eh
Z
c h K : {τ h } ds,
+
J∇uh − ∇u
(2.9)
Eh
Z
Z
Z
2
− f vh dx =
∇h vh : σ h dx −
{∇vh } · [b
σ h ] ds
Ω
Ω
Ehi
Z
−
J∇vh K : {b
σ h } ds.
(2.10)
Eh
9
Taking τ h = (1 − κ)∇2h vh + κ tr (∇2h vh ) I in (2.9), we have
Z
Z
Z
2
2
2
∇h vh : σ h dx = − (1 − κ)∇h uh : ∇h vh dx − κ tr ∇2h uh tr ∇2h vh dx
Ω
Ω
ZΩ
c h } · (1 − κ)[∇2 vh ] + κ[tr ∇2 vh ] ds
+
{∇uh − ∇u
h
h
Ei
Zh
c h K : (1 − κ){∇2 vh } + κ tr {∇2 vh } I ds.
J∇uh − ∇u
+
h
h
Eh
Combining the last equation and (2.10), we obtain an equation which does not rely on σ h
explicitly:
Z
Bh (uh , vh ) =
f vh dx ∀ vh ∈ Vh ,
(2.11)
Ω
where
Z
Bh (uh , vh ) :=
κ)∇2h uh
Z
∇2h vh
(1 −
:
dx + κ tr ∇2h uh tr ∇2h vh dx
ΩZ
Ω
c h − ∇uh } · (1 − κ)[∇2 vh ] + κ[tr ∇2 vh ] ds
+
{∇u
h
h
Ei
Zh
c h − ∇uh K : (1 − κ){∇2 vh } + κ tr {∇2 vh } I ds
+
J∇u
h
h
ZEh
Z
+
{∇vh } · [b
σ h ] ds +
J∇vh K : {b
σ h } ds.
Ehi
(2.12)
Eh
CDG methods for the problem (1.5) can be obtained from (2.11)–(2.12) by proper choices
c h.
b h and ∇u
of numerical traces σ
The relations (2.11)–(2.12) are also the starting point for designing CDG methods for
solving the variational inequality (1.7) through choice of suitable numerical traces to guarantee consistency and stability. For example, taking

c h = {∇uh } on e ∈ Eh ,

∇u



η
b h = −(1 − κ){∇2h uh } − κ tr {∇2h uh } I + J∇uh K on e ∈ Ehi ,
σ
he



 σ
b
b τ = 0, σ
bhν ≤ 0, σ
bhν ∂ν uh = 0 on e ∈ Eh ,
we obtain from (2.11)–(2.12) that
(1)
B1,h (uh , vh )
Z
Z
f vh dx −
=
Ω
Γ
10
σ
bhν ∂ν vh ds,
(2.13)
where
(1)
B1,h (uh , vh )
Z
=
Z
dx + κ tr ∇2h uh tr ∇2h vh dx
:
(1 −
Ω
ΩZ
−
J∇uh K : (1 − κ){∇2h vh } + κ tr {∇2h vh } I ds
Ei
Zh
J∇vh K : (1 − κ){∇2h uh } + κ tr {∇2h uh } I ds
−
Ei
Zh
ηh−1
+
e J∇uh K : J∇vh K ds.
κ)∇2h uh
∇2h vh
(2.14)
Ehi
Here η is a function, which is equal to a constant ηe on each e ∈ Ehi , with {ηe }e∈Ehi having a
uniform positive bound from above and below. In (2.13), let vh = wh − uh with wh ∈ Kh .
Since σ
bhν ≤ 0, σ
bhν ∂ν uh = 0 on e ∈ Ehb , we obtain
Z
(1)
B1,h (uh , wh − uh ) ≥
f (wh − uh ) dx ∀ wh ∈ Kh .
(2.15)
Ω
For a compact formulation, we can use lifting operator r i to get
Z
(1)
B2,h (uh , vh ) = (1 − κ)∇2h uh : ∇2h vh + r i (J∇vh K) dx
ΩZ
+ κ tr ∇2h uh tr ∇2h vh + r i (J∇vh K) dx
ZΩ
+ r i (J∇uh K) : (1 − κ)∇2h vh + κ tr ∇2h vh I dx
ZΩ
ηh−1
+
e J∇uh K : J∇vh K ds.
(2.16)
Ehi
This is the C 0 interior penalty (IP) formulation. A similar C 0 IP method was studied in [7].
The two formulas (2.14) and (2.16) are equivalent on the finite element spaces Vh , so
either form can be used to compute the finite element solution uh . In this paper, we give a
(1)
priori error estimates strictly based on the first formula B1,h . Because of the equivalence of
(1)
these two formulations on Vh , we will prove the stability for the second formula B2,h on Vh ,
(1)
which ensures the stability for the first formulation B1,h on Vh . This comment is valid for
the rest of the CDG methods.
We mow introduce four more CDG methods for the variational inequality (1.7). The
methods are all of the form (2.15), and so we will only list the corresponding bilinear form.
Comparing with the DG methods for the second order elliptic problem, we can give the
11
C 0 non-symmetric interior penalty (NIPG) formulations,
Z
Z
(2)
2
2
B1,h (uh , vh ) = (1 − κ)∇h uh : ∇h vh dx + κ tr ∇2h uh tr ∇2h vh dx
ΩZ
Ω
J∇uh K : (1 − κ){∇2h vh } + κ tr {∇2h vh } I ds
+
Ei
Zh
J∇vh K : (1 − κ){∇2h uh } + κ tr {∇2h uh } I ds
−
Ei
Zh
ηh−1
+
e J∇uh K : J∇vh K ds,
Ehi
or equivalently,
(2)
B2,h (uh , vh )
Z
=
(1 − κ)∇2h uh : ∇2h vh + r i (J∇vh K) dx
ΩZ
+ κ tr ∇2h uh tr ∇2h vh + r i (J∇vh K) dx
ZΩ
− r i (J∇uh K) : (1 − κ)∇2h vh + κ tr ∇2h vh I dx
ZΩ
+
ηh−1
e J∇uh K : J∇vh K ds.
Ehi
That is, we solve the variational inequality
Z
(2)
f (wh − uh ) dx ∀ wh ∈ Kh .
B1,h (uh , wh − uh ) ≥
(2.17)
Ω
Using the local lifting operator r e , we can give the third example. Taking

c h = {∇uh } on e ∈ Eh ,

∇u



 σ
b h = − (1 − κ){∇2h uh } − κ tr({∇2h uh })I − (1 − κ){r i (J∇uh K)} − κ{tr(r i (J∇uh K))I}

− (1 − κ){η r e (J∇uh K)} − κ{η tr(r e (J∇uh K))}I on e ∈ Ehi ,




b τ = 0, σ
σ
bhν ≤ 0, σ
bhν ∂ν uh = 0 on e ∈ Ehb ,
12
we get from (2.12) that
Z
Z
(3)
2
2
B1,h (uh , vh ) = (1 − κ)∇h uh : ∇h vh dx + κ tr ∇2h uh tr ∇2h vh dx
ΩZ
Ω
2
−
J∇uh K : (1 − κ){∇h vh } + κ tr {∇2h vh } I ds
Ei
Zh
−
J∇vh K : (1 − κ){∇2h uh } + κ tr {∇2h uh } I ds
Ei
Zh
+ r i (J∇vh K) : (1 − κ)r i (J∇uh K) + κ tr(r i (J∇uh K))I dx
Ω
XZ
+
η ((1 − κ)r e (J∇uh K) : r e (J∇vh K) + κ tr(r e (J∇uh K))tr(r e (J∇vh K))) dx,
e∈Ehi
Ω
or equivalently,
Z
(3)
B2,h (uh , vh ) = (1 − κ) ∇2h uh + r i (J∇uh K) : ∇2h vh + r i (J∇vh K) dx
ΩZ
+ κ tr ∇2h uh + r i (J∇uh K) tr ∇2h vh + r i (J∇vh K) dx
Ω
XZ
+
η ((1 − κ)r e (J∇uh K) : r e (J∇vh K) + κ tr(r e (J∇uh K))tr(r e (J∇vh K))) dx,
e∈Ehi
Ω
which is the CDG formulation proposed in [24]. That is, we solve the variational inequality
Z
(3)
B1,h (uh , wh − uh ) ≥
f (wh − uh ) dx ∀ wh ∈ Kh .
(2.18)
Ω
With the choice of

c h = {∇uh } on e ∈ Eh ,

∇u



 σ
b h = − (1 − κ){∇2h uh } − κ tr({∇2h uh })I − (1 − κ){η r e (J∇uh K)}

− κ{η tr(r e (J∇uh K))}I on e ∈ Ehi ,




b τ = 0, σ
σ
bhν ≤ 0, σ
bhν ∂ν uh = 0 on e ∈ Ehb ,
we obtain
(4)
B1,h (uh , vh )
Z
=
κ)∇2h uh
Z
(1 −
:
dx + κ tr ∇2h uh tr ∇2h vh dx
ΩZ
Ω
−
J∇uh K : (1 − κ){∇2h vh } + κ tr {∇2h vh } I ds
Ei
Zh
−
J∇vh K : (1 − κ){∇2h uh } + κ tr {∇2h uh } I ds
Ehi
XZ
+
η ((1 − κ)r e (J∇uh K) : r e (J∇vh K) + κ tr(r e (J∇uh K))tr(r e (J∇vh K))) dx,
e∈Ehi
∇2h vh
Ω
13
or equivalently,
Z
(4)
B2,h (uh , vh ) = (1 − κ)∇2h uh : ∇2h vh + r i (J∇vh K) dx
ΩZ
+ κ tr ∇2h uh tr ∇2h vh + r i (J∇vh K) dx
ZΩ
+ r i (J∇uh K) : (1 − κ)∇2h vh + κ tr ∇2h vh I dx
Ω
XZ
+
η ((1 − κ)r e (J∇uh K) : r e (J∇vh K) + κ tr(r e (J∇uh K))tr(r e (J∇vh K))) dx,
e∈Ehi
Ω
which is the CDG formulation extended from the DG method of [6] for elliptic problem of
second order. That is, we solve the variational inequality
Z
(4)
f (wh − uh ) dx ∀ wh ∈ Kh .
(2.19)
B1,h (uh , wh − uh ) ≥
Ω
Choosing

c h = {∇uh } on e ∈ Eh ,

∇u



 σ
b h = − (1 − κ){∇2h uh } − κ tr({∇2h uh })I − (1 − κ){r i (J∇uh K)} − κ{tr(r i (J∇uh K))I}
i

+ ηh−1

e J∇uh K on e ∈ Eh ,



b τ = 0, σ
σ
bhν ≤ 0, σ
bhν ∂ν uh = 0 on e ∈ Ehb ,
we get the LCDG method ([15]),
Z
Z
(5)
2
2
B1,h (uh , vh ) := (1 − κ)∇h uh : ∇h vh dx + κ tr ∇2h uh tr ∇2h vh dx
ΩZ
Ω
−
J∇uh K : (1 − κ){∇2h vh } + κ tr {∇2h vh } I ds
Ei
Zh
−
J∇vh K : (1 − κ){∇2h uh } + κ tr {∇2h uh } I ds
Ei
Zh
+ r i (J∇vh K) : (1 − κ)r i (J∇uh K) + κ tr(r i (J∇uh K))I dx
ZΩ
+
ηh−1
e J∇uh K : J∇vh K ds,
Ehi
or equivalently,
(5)
B2,h (uh , vh )
Z
:=
(1 − κ) ∇2h uh + r i (J∇uh K) : ∇2h vh + r i (J∇vh K) dx
ΩZ
+ κ tr ∇2h uh + r i (J∇uh K) tr ∇2h vh + r i (J∇vh K) dx
ZΩ
+
ηh−1
e J∇uh K : J∇vh K ds.
Ehi
14
That is, we solve the variational inequality
Z
(5)
B1,h (uh , wh − uh ) ≥
f (wh − uh ) dx ∀ wh ∈ Kh .
(2.20)
Ω
In the following sections, we will study CDG methods for the EVI (1.7) defined as follows:
Find uh ∈ Kh such that
Bh (uh , vh − uh ) ≥ (f, vh − uh ) ∀ vh ∈ Kh ,
(2.21)
(j)
where the bilinear form Bh (w, v) = B1,h (w, v) with j = 1, · · · , 5.
3
Consistency, boundedness and stability
We present some properties of the five DG methods introduced in Section 2. First, we
address the consistency of the methods (2.21).
Lemma 3.1 Assume u ∈ H 3 (Ω) is the solution of the problem (1.7). For all the five CDG
(j)
methods Bh (w, v) = B1,h (w, v) with j = 1, · · · , 5, we have
Bh (u, vh − u) ≥ (f, vh − u)
∀ vh ∈ Kh .
Proof. Noting J∇uK = 0 on each edge e ∈ Ehi , we use (1.2) to get
Z
Z
2
2
Bh (u, vh − u) = (1 − κ)∇ u : ∇h (vh − u) dx + κ tr ∇2 u tr ∇2h (vh − u) dx
ΩZ
Ω
−
J∇(vh − u)K : (1 − κ)∇2 u + κ tr ∇2 u I ds
Ehi
Z
XZ
2
=−
σ : ∇h (vh − u) dx +
J∇(vh − u)K : σ ds.
T ∈Th
Ehi
T
Using Lemma 1.1 and noticing [σ] = 0 on each edge e ∈ Ehi , we have
XZ
XZ
2
σ : ∇h (vh − u) dx = −
∇(vh − u) · (∇ · σ) dx
T ∈Th
T
T
T ∈Th
+
XZ
T ∈Th
∇(vh − u) · (σnT ) ds
∂T
Z
=−
Z
∇(vh − u) · (∇ · σ) dx +
Eh
Ω
15
J∇(vh − u)K : σ ds.
Combining the above two equations,
Z
Z
Bh (u, vh − u) =
∇(vh − u) · (∇ · σ) dx − J∇(vh − u)K : σ ds.
Ω
Γ
Since vh − u = 0 on Γ and ∂ν vh ≥ 0 on Γ for all vh ∈ Kh , we use Lemma 1.1 and (1.9) to
obtain
Z
Z
Bh (u, vh − u) = − (vh − u)∇ · (∇ · σ) dx − ∇(vh − u) · (σν) ds
Γ
Z
Z Ω
= f (vh − u) dx − σν ∂ν vh ds
Γ
ZΩ
≥ f (vh − u) dx.
Ω
So the stated result holds.
Let V (h) := Vh + H01 (Ω) ∩ H 3 (Ω) and define two mesh-dependent energy norms by
X
X
2
2
2
h2T |v|23,T , v ∈ V (h).
|v|2∗ = |v|22,h +
h−1
kJ∇vKk
,
9v9
=
|v|
+
e
0,e
∗
T ∈Th
e∈Ehi
To show these formulas define norms, we only need prove that |v|∗ = 0 and v ∈ V (h) imply
v = 0. From |v|2,h = 0, we have v|T ∈ P1 (T ) and so ∇v is piecewise constant. Let e be a
common edge of two elements T + and T − . From kJ∇vKk0,e = 0, we obtain (∇v)+ = (∇v)− .
Thus, ∇v is constant in Ω and so v ∈ P1 (Ω). Since v = 0 on Γ, we conclude v = 0 in Ω.
Before presenting boundedness and stability results of the bilinear forms, we give a useful
estimate for the lifting operator r e .
Lemma 3.2 For any v ∈ V (h) and e ∈ Ehi ,
2
2
−1
2
C1 h−1
e kJ∇vKk0,e ≤ kr e (J∇vK)k0,h ≤ C2 he kJ∇vKk0,e .
(3.1)
Proof. The second inequality was proved in [15]. For v ∈ H 3 (Ω), J∇vK = 0 on e ∈ Ehi . So
we only need to consider the case v ∈ Vh . By the formula between (4.4) and (4.5) in [2], we
know5
∗
2
−1
2
2
h−1
∀ ϕ ∈ [P1 (e)]2 ,
(3.2)
e kϕk0,e . kre (ϕ)k0,Ω . he kϕk0,e
where the lifting operator re∗ : (L2 (e))2 → Wh is defined by
Z
Z
∗
re (v) · wh dx = − v · {wh } ds, ∀wh ∈ Wh .
Ω
5
The definition of
re∗
e
is given, please check.
16
n
o
2
2
Here, Wh := wh ∈ (L (Ω)) : whi |K ∈ Pl (K), ∀ K ∈ Th , i = 1, 2 .
For two matrix-valued functions φ = (φij )2×2 and τ = (τij )2×2 , let φ1 = (φ11 , φ21 )t ,
φ2 = (φ12 , φ22 )t , τ 1 = (τ11 , τ21 )t , τ 2 = (τ12 , τ22 )t , so that φ = (φ1 , φ2 ), τ = (τ 1 , τ 2 ). Then
Z
Z
Z
Z
r e (φ) : τ dx = − φ : {τ } ds = − φ1 · {τ 1 } ds − φ2 · {τ 2 } ds
Ω
eZ
Z e
Z e
=
re∗ (φ1 ) · τ 1 dx + re∗ (φ2 ) · τ 2 dx = (re∗ (φ1 ), re∗ (φ2 )) : τ dx,
Ω
Ω
Ω
for all τ ∈ σ h . So r e (φ) = (re∗ (φ1 ), re∗ (φ2 )), kr e (φ)k20,Ω = kre∗ (φ1 )k20,Ω + kre∗ (φ2 )k20,Ω , and
−1
2
2
∗
2
∗
2
2
2
h−1
e kφk0,e = he (kφ1 k0,e + kφ2 k0,e ) . kre (φ1 )k0,Ω + kre (φ2 )k0,Ω = kr e (φ)k0,Ω .
Let φ = J∇vK, the first inequality is proved.
From (3.1) and (2.5), we have
kr i (J∇vK)k20,h = k
X
r e (J∇vK)k20,h ≤ 3C2
e∈Ehi
X
2
h−1
e kJ∇vKk0,e .
e∈Ehi
(j)
For the boundedness of the primal forms B1,h with j = 1, · · · , 5, first notice that
ktr(τ )k0,h . kτ k0,h . By the Cauchy-Schwarz inequality and Lemma 3.2, we get the following inequalities:
Z
∇2h w : ∇2h v dx ≤ |w|2,h |v|2,h ,
(3.3)
Ω
1/2 

Z
r i (J∇wK) : r i (J∇vK) dx . 
Ω
X
2 
h−1
e kJ∇wKk0,e
Ehi
e∈Ehi
,
(3.4)
1/2
1/2 
X
2 
h−1
e kJ∇wKk0,e
e∈Ehi
η r e (J∇wK) : r e (J∇vK) dx . sup ηe 
e∈Ehi
X

2 
h−1
e kJ∇vKk0,e
, (3.5)
e∈Ehi
1/2 

Ω
2 
h−1
e kJ∇vKk0,e
e∈Ehi

ηh−1
e J∇wK : J∇vK ds ≤ sup ηe
XZ
e∈Ehi

e∈Ehi

Z
1/2
X
X
e∈Ehi
2 
h−1
e kJ∇wKk0,e
1/2
X

2 
h−1
e kJ∇vKk0,e
e∈Ehi
(3.6)
17
.
2
2
Using the trace inequality k∇2 vk0,e . h−1
e |v|2,K + he |v|3,K with e an edge of K, we have
Z
XZ
2
J∇wK : {∇h v} ds =
J∇wK : {∇2h v} ds
Ehi
e∈Ehi
e
1/2 

≤
X
2 
h−1
e kJ∇wKk0,e
e∈Ehi
X

he k{∇2h v}k20,e 
e∈Ehi
1/2

.
1/2
X
he−1 kJ∇wKk20,e 
!1/2
X
(|v|22,T + h2T |v|23,T )
.
(3.7)
T ∈Th
e∈Ehi
The inequalities (3.3) and (3.7) are needed by all bilinear forms. For the CDG methods
(j)
with the bilinear form B1,h , j = 1, 2, 5, the inequality (3.5) is needed. The inequality (3.4)
(j)
(j)
is needed by the formulas B1,h with j = 3, 5. The methods with the bilinear forms B1,h ,
j = 3, 4, need the inequality (3.6). Then we have the following result.
(j)
Lemma 3.3 (Boundedness) Let Bh = B1,h with j = 1, · · · , 5. Then
Bh (w, v) . 9w 9 9v 9
∀ (w, v) ∈ V (h) × V (h).
(3.8)
(j)
(j)
For stability over Vh , note that 9v9 = |v|∗ for any v ∈ Vh . Formulations B1,h and B2,h
(j)
are equivalent on Vh , so we just need to prove the stability for B2,h based on | · |∗ . We use
the Cauchy-Schwarz inequality and Lemma 3.2 to get
Z
Z
Z
2
(1)
2
2
2
B2,h (v, v) = (1 − κ) ∇h v : ∇h v dx + κ
tr(∇h v) dx + 2(1 − κ) ∇2h v : r i (J∇vK) dx
Ω
Ω
Z
Z Ω
2
ηh−1
+ 2κ tr ∇2h v tr (r i (J∇vK)) dx +
e |J∇vK| ds
Ei
Ω
h
1
2
2
2
2
≥ (1 − κ)|v|2,h + κk∆h vk0,h − (1 − κ) |v|2,h + kr i (J∇vK)k0,h
X
2
2
2
− κ k∆h vk0,h + ktr (r i (J∇vK)) k0,h + η0
h−1
e kJ∇vKk0,e
e∈E i
h
X
3(1 − κ)C2
2
2
≥ (1 − )(1 − κ)|v|2,h + η0 −
− 6C2 κ
h−1
e kJ∇vKk0,e ,
i
e∈Eh
18
where 0 < < 1 is a constant and C2 is the generic positive constant in (3.1). Therefore,
stability is valid for C 0 IP method when mine∈Ehi ηe = η0 > 3(1 − κ)C2 + 6C2 κ = 3(1 + κ)C2 .
Z
Z
Z
2
(2)
2
2
2
2
B2,h (v, v) = (1 − κ)∇h v : ∇h v dx + κ tr(∇h v) dx +
ηh−1
e (J∇vK) ds
Ω
Ehi
Ω
≥ (1 − κ)|v|22,h + η0
X
2
h−1
e kJ∇vKk0,e .
e∈Ehi
So stability is valid for the C 0 NIPG method for any η0 > 0. This property is the reason
(2)
(2)
why the method with the bilinear form B2,h is useful even though B2,h is not symmetric.
Z
(4)
2
2
B2,h (v, v) ≥ (1 − κ)|v|2,h + κk∆h vk0,h + 2(1 − κ) ∇2h v : r i (J∇vK) dx
Ω
Z
X
+ 2κ ∆h v tr (r i (J∇vK)) dx + η0
(1 − κ)kr e (J∇vK)k20,h + κktr(r e (J∇vK))k20,h
Ω
e∈Ehi
1
2
≥ (1 −
+
− (1 − κ)
+ kr i (J∇vK)k0,h − κk∆h vk20,h
X
X
2
− κktr (r i (J∇vK)) k20,h + η0 C1 (1 − κ)
h−1
kJ∇vKk
+
η
κ
ktr(r e (J∇vK))k20,h
0
e
0,e
κ)|v|22,h
κk∆h vk20,h
|v|22,h
e∈E i
e∈Ehi
h
3C2 X −1
2
he kJ∇vKk20,e
≥ (1 − )(1 − κ)|v|2,h + (1 − κ) η0 C1 −
e∈Ehi
X
+ (η0 κ − 3κ)
ktr(r e (J∇vK))k20,h .
e∈Ehi
Since C2 > C1 , η0 > 3 is guaranteed from η0 > 3C2 /C1 . Thus, stability is valid for this CDG
formulation when η0 > 3C2 /C1 . For the method of Wells-Dung corresponding to the form
(3)
(5)
B2,h and the LCDG method corresponding to the form B2,h , stability can be analyzed by a
similar argument (cf. [24] and [15], respectively), with η0 > 0.
(j)
Lemma 3.4 (Stability) Let Bh = B2,h with j = 1, · · · , 5. Assume
mini ηe > 3 (1 + κ) C2 for j = 1
e∈Eh
and
mini ηe > 3 C2 /C1 for j = 4,
e∈Eh
with C1 and C2 the constants in the inequality (3.1). Then,
9v92 . Bh (v, v)
∀ v ∈ Vh .
(3.9)
It follows from Lemma 3.4 that under the stated conditions, the stability property is also
(j)
valid for B1,h with j = 1, · · · , 5.
19
4
Numerical results
In this section, we present a numerical example on the five methods studied in solving the
elliptic variational inequality (1.7). To solve the discretized variational inequality, we use
the optimization toolbox in MATLAB for the associated quadratic optimization problem.
Let Ω = (−1, 1) × (−1, 1), κ = 0.3. Choose the right hand side function to be
f (x, y) = 24(1 − x2 )2 + 24(1 − y 2 )2 + 32(3x2 − 1)(3y 2 − 1).
We use uniform triangulations Th of the region Ω and piecewise quadratic polynomials, i.e.,
Vh = {vh ∈ H01 (Ω) : vh |T ∈ P2 (T ) ∀ T ∈ Th }.
Since the domain is rectangular, the outward normal is not defined at the four corner points,
(−1, −1), (−1, 1), (1, −1), and (1, 1). So in terms of the restriction in the optimization
problem, at each of the four corners, we specify the constraint ∂ν vh ≥ 0 twice (corresponding
to the two sides intersecting at the corner).
Tables 1–5 provide numerical solution errors in the energy norm 9 · 9 for the five DG
methods discussed in this paper. Since the true solution of the variational inequality (1.7) is
not known, we use the numerical solution corresponding to the meshsize h = 0.0625 as the
true solution in computing the errors.
Table 1: Error 9u − uh 9 for C 0 IP method (2.15)
h
η=1
η = 10
η = 100
η = 1000
1
6.1327
4.4264
2.7621
1.5772
0.5
3.6713
2.0188
1.1607
0.6544
0.25
1.7921
1.0202
0.5767
0.3245
0.125
0.8921
0.5076
0.2860
0.1628
Table 2: Error 9u − uh 9 for NIPG method (2.17)
h
η=1
η = 10
η = 100
η = 1000
1
5.5212
4.3318
2.7521
1.5766
0.5
3.2523
2.0225
1.1610
0.6545
20
0.25
1.7682
1.0215
0.5767
0.3245
0.125
0.8987
0.5082
0.2860
0.1628
Table 3: Error 9u − uh 9 for Wells-Dung DG formulation (2.18) as in [24]
h
η=1
η = 10
η = 100
η = 1000
1
5.1659
4.2043
2.9413
1.7224
0.5
3.1617
2.3618
1.4600
0.9311
0.25
2.3018
1.6532
0.9954
0.5570
0.125
1.6206
1.1411
0.7116
0.3836
Table 4: Error 9u − uh 9 for Baker-DG formulation (2.19) as in [5]
h
η=1
η = 10
η = 100
η = 1000
1
5.1509
4.1924
2.8889
1.7181
0.5
3.1566
2.3418
1.4448
0.9300
0.25
2.3016
1.6425
0.9818
0.5569
0.125
1.6204
1.1408
0.6774
0.3837
Table 5: Error 9u − uh 9 for LCDG method (2.20)
h
η=1
η = 10
η = 100
η = 1000
1
5.4865
4.2776
2.7463
1.5762
0.5
3.0717
1.9847
1.1584
0.6543
0.25
1.7165
1.0158
0.5764
0.3244
0.125
0.8831
0.5071
0.2859
0.1608
We observe that all the five DG methods perform well, and the methods 1, 2, and 5
converge faster than the methods 3 and 4. For the methods 1, 2, and 5, the numerical
convergence order in the norm 9u − uh 9 is around 1, whereas for the methods 3 and 4, the
numerical convergence order in the norm 9u − uh 9 is around 1/2.
References
[1] D.N. Arnold, F. Brezzi, B. Cockburn, and L.D. Marini, Discontinuous Galerkin methods
for elliptic problems, in Discontinuous Galerkin Methods: Theory, Computation and
21
Applications, B. Cockburn, G.E. Karniadakis, and C.-W. Shu, eds., Lecture Notes in
Comput. Sci. Engrg. 11, Springer-Verlag, New York, 2000, 89–101.
[2] D.N. Arnold, F. Brezzi, B. Cockburn, and L.D. Marini, Unified analysis of discontinuous
Galerkin methods for elliptic problems, SIAM J. Numer. Anal. 39 (2002), 1749–1779.
[3] K. Atkinson and W. Han, Theoretical Numerical Analysis: A Functional Analysis
Framework, third edition, Springer-Verlag, New York, Texts in Applied Mathematics,
Volume 39, 2009.
[4] I. Babu˘ska and M. Zl´amal, Nonconforming elements in the finite element method with
penalty, SIAM J. Numer. Anal. 10 (1973), 863–875.
[5] G.A. Baker, Finite element methods for elliptic equations using nonconforming elements,
Math. Comp. 31 (1977), 44–59.
[6] F. Bassi, S. Rebay, G. Mariotti, S. Pedinotti, and M. Savini, A high-order accurate
discontinuous finite element method for inviscid and viscous turbomachinery flows, in
Proceedings of 2nd European Conference on Turbomachinery, Fluid Dynamics and Thermodynamics, R. Decuypere and G. Dibelius, eds., Technologisch Instituut, Antwerpen,
Belgium, 1997, 99–108.
[7] S.C. Brenner and L. Sung, C 0 interior penalty methods for fourth order elliptic boundary
value problems on polygonal domains, J. Sci. Comput. 22/23 (2005), 83–118.
[8] P. Castillo, B. Cockburn, I. Perugia, and D. Sch¨otzau, An a priori error analysis of the
local discontinuous Galerkin method for elliptic problems, SIAM J. Numer. Anal. 18
(2000), 1676–1706.
[9] B. Cockburn, Discontinuous Galerkin methods, ZAMM Z. Angew. Math. Mech. 83
(2003), 731–754.
[10] B. Cockburn, G.E. Karniadakis, and C.-W. Shu, eds., Discontinuous Galerkin Methods. Theory, Computation and Applications, Lecture Notes in Comput. Sci. Engrg. 11,
Springer-Verlag, New York, 2000.
[11] G. Duvaut and J.-L. Lions, Inequalities in Mechanics and Physics, Springer-Verlag,
Berlin, 1976.
[12] G. Engel, K. Garikipati, T. Hughes, M. Larson, L. Mazzei, and R. Taylor, Continuous/discontinuous finite element approximations of fourth-order elliptic problems in
structural and continuum mechanics with applications to thin beams and plates, and
strain gradient elasticity, Comput. Methods Appl. Mech. Engrg. 191 (2002), 3669–3750.
22
[13] R. Glowinski, Numerical Methods for Nonlinear Variational Problems, Springer-Verlag,
New York, 1984.
[14] R. Glowinski, J.-L. Lions, and R. Tr´emoli`eres, Numerical Analysis of Variational Inequalities, North-Holland, New York, 1981.
[15] J. Huang, X. Huang, and W. Han, A new C 0 discontinuous Galerkin method for Kirchhoff plates, Comput. Methods Appl. Mech. Engrg. 199 (2010), 1446–1454.
[16] I. Mozolevski and P.R. B¨osing, Sharp expressions for the stabilization parameters in
symmetric interior-penalty discontinuous Galerkin finite element approximations of
fourth-order elliptic problems, Comput. Methods Appl. Math. 7 (2007), 365–375.
[17] I. Mozolevski and E. S¨
uli, A priori error analysis for the hp-version of the discontinuous
Galerkin finite element method for the biharmonic equation, Comput. Methods Appl.
Math. 3 (2003), 596–607.
[18] I. Mozolevski, E. S¨
uli, and P.R. B¨osing, hp-version a priori error analysis of interior
penalty discontinuous Galerkin finite element approximations to the biharmonic equation, J. Sci. Comput. 30 (2007), 465–491.
[19] J.N. Reddy, Theory and Analysis of Elastic Plates and Shells, second edition, CRC
Press, New York, 2007.
[20] E. S¨
uli and I. Mozolevski, hp-version interior penalty DGFEMs for the biharmonic
equation, Comput. Methods Appl. Mech. Engrg. 196 (2007), 1851–1863.
[21] F. Wang, W. Han, and X. Cheng, Discontinuous Galerkin methods for solving elliptic
variational inequalities, SIAM Journal on Numerical Analysis 48 (2010), 708–733.
[22] F. Wang, W. Han, and X. Cheng, Discontinuous Galerkin methods for solving Signorini
problem, IMA Journal of Numerical Analysis 31 (2011), 1754–1772.
[23] F. Wang, W. Han, and X. Cheng, Discontinuous Galerkin methods for solving the
quasistatic contact problem, Numerische Mathematik (2013), DOI 10.1007/s00211-0130574-0.
[24] G.N. Wells and N.T. Dung, A C 0 discontinuous Galerkin formulation for Kirchhoff
plates, Comput. Methods Appl. Mech. Engrg. 196 (2007), 3370–3380.
23

第1章 序論