fractional order Sobolev spaces, Poisson, Signorini, mixed method, bubble finite element stabiliza- tion

SOME NONSTANDARD FINITE ELEMENT ESTIMATES WITH APPLICATIONS TO 3D POISSON AND SIGNORINI PROBLEMS^∗

FAKER BEN BELGACEM^†ANDSUSANNE C. BRENNER^‡

Abstract. In this paper we establish several nonstandard finite element estimates involving fractional order Sobolev spaces, with applications to bubble stabilized mixed methods for the three-dimensional Poisson and Sig- norini problems.

Key words. fractional order Sobolev spaces, Poisson, Signorini, mixed method, bubble finite element stabiliza- tion.

AMS subject classifications. 65N30, 35J85.

1. Introduction. It is well-known that interpolation error estimates play an important role in the analysis of finite element methods. The simplest interpolation error estimates appear in the following form:

m

X

j=0

h^j|u−Πhu|_H^j(Ω)≤Ch^m|u|_H^m(Ω), (1.1)

whereΩis a bounded polyhedral domain inR^d(d= 1,2,3),mis an integer andΠhis an interpolation operator fromH^m(Ω)to the finite element spaceVhassociated with a regular triangulationThofΩof mesh-sizeh, and the seminorm| · |_H^k_(Ω)for a nonnegative integer is defined by

|v|²_Hk(Ω)= X

|α|=k

k∂^αvk²_L2(Ω).

In the case whereΠhis defined locally on each element, since all the seminorms in (1.1) are also local, the estimate (1.1) can be established by a purely local analysis (cf. [12], [8]).

For applications to problems whose solutions are not regular, it is important to have estimates for(u−Πhu) in fractional order Sobolev seminorms. Letk be a nonnegative integer and0< λ <1. The seminorm| · |_Hk+λ(Ω)is defined by (cf. [1], [22])

|v|²_Hk+λ(Ω)= X

|α|=k

Z

Ω

Z

Ω

[∂^αv(x)−∂^αv(y)]²

|x−y|^d+2λ dxdy, (1.2)

and the normk · k_H^k+λ_(Ω)is given by

kvk²_Hk+λ(Ω)=kvk²_Hk(Ω)+|v|²_Hk+λ(Ω)=

k

X

j=0

|v|²_Hj(Ω)+|v|²_Hk+λ(Ω).

∗Received November 10, 2000. Accepted for publication April 11, 2001. Recommended by R. Varga.

†Math´ematiques pour l’Industrie et la Physique, Unit´e Mixte de Recherche CNRS–UPS–INSAT–UT1 (UMR 5640), Universit´e Paul Sabatier, 118 route de Narbonne, 31062 TOULOUSE CEDEX, FRANCE. E-mail:

[email protected]

‡Department of Mathematics, University of South Carolina, Columbia, SC 29208, USA. This research was sup- ported in part by the National Science Foundation under Grant No. DMS-0074246. E-mail:

[email protected]

134

Sometimes such estimates follow easily from the basic estimate (1.1) and the interpola- tion theory of Sobolev spaces (cf. [22]). For example, letΠ = Π¹_hbe the nodal interpolation operator for theP1finite element space and0< λ <1, then (1.1) withm= 2implies that

ku−Π¹_hukL²(Ω)≤Ch²kukH²(Ω) ∀u∈H²(Ω), ku−Π¹_hukH¹(Ω)≤ChkukH²(Ω) ∀u∈H²(Ω), and hence

ku−Π¹_huk_H^λ_(Ω)≤Cλh^2−λkuk_H²_(Ω) ∀u∈H²(Ω).

Replacinguwith(u−p)for an appropriate first order polynomialpso that Z

Ω

∂^α(u−p)dx= 0 for |α| ≤1,

we can then deduce from the preceding estimate and the Friedrichs inequality (cf. [19]) that

|u−Π¹_hu|_H^λ_(Ω)≤Cλh^2−λ|u|_H²(Ω).

Note that throughout this paper we use the symbolC(with or without subscripts) to represent generic positive constants which can take different values at different places.

On the other hand, in the case whereΠ = Π⁰_his the piecewiseL²-orthogonal projection operator fromL²(Ω)to theP0finite element space, the standard estimate (1.1) withm= 1 only gives

ku−Π⁰_huk_L²_(Ω)≤Ch|u|_H¹_(Ω). For0< λ < ¹₂, the estimate

|u−Π⁰_hu|_Hλ(Ω)≤C_λh^1−λ|u|_H¹_(Ω) (1.3)

does not follow from the interpolation theory of Sobolev spaces.

Similarly, for0< λ < ¹₂, the estimate

|u−Π¹_hu|_H^1+λ_(Ω)≤Cλh^1−λ|u|_H²(Ω)

(1.4)

does not follow from (1.1) and interpolation.

Inverse estimates are also important in the analysis of finite element methods. A standard inverse estimate for a Lagrange finite element spaceV_h(Ω)⊂H¹(Ω)takes the form

|v|_H¹_(Ω)≤Ch⁻¹kvk_L²_(Ω) ∀v∈Vh(Ω), (1.5)

and it can be obtained by a purely local analysis (cf. [12], [8]).

Letλ∈(0,1). From (1.5) and the trivial estimate

kvkL²(Ω)≤ kvkL²(Ω) ∀v∈V_h(Ω),

we can obtain by interpolation the following inverse estimate for the seminorm| · |_H^λ_(Ω):

|v|_H^λ(Ω)≤Cλh^−λkvkL²(Ω) ∀v∈Vh(Ω).

However, for0< λ < ¹₂, the inverse estimate

|v|_H^1+λ_(Ω)≤Cλh^−λ|v|_H¹(Ω) ∀v∈Vh(Ω) (1.6)

does not follow from (1.5) and interpolation.

In this paper we will establish certain finite element estimates in fractional order Sobolev seminorms which include in particular (1.3), (1.4) and (1.6), with applications to some three dimensional mixed finite element methods.

The main difficulty in dealing with estimates in fractional order Sobolev seminorms is due to the non-local nature of the definition (1.2). Our key observation is that the analysis of such estimates can be reduced to a purely local one by combining the ideas from [7] and the estimate

Z

Ω\T

1

|x−y|^d+2λ dy≤ Cλ

ρ(x, ∂T)^2λ ∀x∈T, (1.7)

whereT is an element inT_handρ(x, ∂T) = infy∈∂T|x−y|is the distance fromxto the boundary ofT. Furthermore, the local analysis can be handled by the following estimate on a reference domainTˆ:

Z

Tˆ

u²(x)

ρ(x, ∂Tˆ)^2λ dx≤CT ,λˆ kuk²_Hλ( ˆT) ∀u∈H^λ( ˆT) and 0< λ < 1 2. (1.8)

The estimate (1.7) can be obtained easily by using the regularity ofT_hand a direct calcula- tion. The estimate (1.8), which follows from the Hardy inequalities, comes from the theory of Sobolev spaces. A proof of it can be found in either [18] or [15].

The rest of the paper is organized as follows. The nonstandard finite element estimates in fractional order Sobolev norms are proved in Section 2. Applications to three-dimensional mixed finite element methods for the Poisson problem and the Signorini contact problem are then given in Section 3.

2. Nonstandard finite element estimates in fractional order Sobolev norms. LetΩ be a bounded polyhedral domain inR^dford= 1,2,3. LetT_hbe a regular triangulation (cf.

[12], [8]) ofΩ, whereh= maxT∈T_hdiamT. We will first show how certain estimates for the globally defined fractional order Sobolev seminorms can be reduced to local estimates.

LEMMA 2.1. Letλ ∈ (0,1)andw ∈ H^k+λ(Ω). Then the following error estimate holds:

|w|²_Hk+λ(Ω)≤Cλ

X

|α|=k

X

T∈T_h

|∂^αw|²_Hλ(T)+ Z

T

[∂^αw(x)]² ρ(x, ∂T)^2λ dx

.

Proof. By the definition (1.2), we have

|w|²_Hk+λ(Ω)= X

|α|=k

Z

Ω

Z

Ω

[∂^αw(x)−∂^αw(y)]²

|x−y|^d+2λ dxdy

= X

|α|=k

X

T∈T_h

Z

T

Z

T

[∂^αw(x)−∂^αw(y)]²

|x−y|^d+2λ dxdy (2.1)

+ X

T,T⁰∈T_h T6=T⁰

Z

T

Z

T⁰

[∂^αw(x)−∂^αw(y)]²

|x−y|^d+2λ dxdy

.

Note that the first sum inside the bracket on the right-hand side of (2.1) equals P

T∈T_h|∂^αw|²_Hλ(T). The second sum can be estimated following the ideas in [7]:

X

T,T⁰∈T_h T6=T⁰

Z

T

Z

T⁰

[∂^αw(x)−∂^αw(y)]²

|x−y|^d+2λ dxdy≤2

X

T,T⁰∈T_h T6=T⁰

Z

T

Z

T⁰

[∂^αw(x)]²

|x−y|^d+2λ dxdy (2.2)

+ X

T,T⁰∈Th

T6=T⁰

Z

T

Z

T⁰

[∂^αw(y)]²

|x−y|^d+2λ dxdy

.

Let us focus on the first term on the right-hand side of (2.2), since the second one can be worked out in exactly the same way. We have, by (1.7),

X

T,T⁰∈T_h T6=T⁰

Z

T

Z

T⁰

[∂^αw(x)]²

|x−y|^d+2λ dxdy= X

T∈T_h

Z

T

[∂^αw(x)]²Z

Ω\T

1

|x−y|^d+2λ dy dx (2.3)

≤C_λ X

T∈T_h

Z

T

[∂^αw(x)]² ρ(x, ∂T)^2λ dx.

The lemma follows from (2.1)–(2.3).

REMARK2.2. Note that the right-hand side of the estimate in Lemma2.1is in general infinite whenλ∈[¹₂,1). Therefore Lemma2.1is of interest mainly forλ∈(0,¹₂).

Let us now consider estimates for the piecewise constant interpolation operatorΠ⁰_h. Let P0(T)stand for the set of constant functions defined onT ∈ T_h. The P0 finite element subspaceM_h(Ω)is defined by

M_h(Ω) =n

v_h∈L²(Ω) :v_h

_T ∈ P0(T) ∀T ∈T_ho .

The piecewise constant interpolation operatorΠ⁰_his the orthogonal projection fromL²(Ω) intoMh(Ω)and is defined as follows:

Z

Ω

(v−Π⁰_hv)ψhdx= 0 ∀v∈L²(Ω), ψh∈Mh(Ω).

(2.4)

THEOREM2.3. For anyλ∈(0,¹₂)andµ∈[λ,1], the following error estimate holds:

|v−Π⁰_hv|_Hλ(Ω)≤C_λh^µ−λ|v|_H^µ_(Ω) ∀v∈H^µ(Ω).

(2.5)

Proof. Setv_h= Π⁰_hv. From Lemma2.1we have

|v−vh|²_Hλ(Ω)≤Cλ

X

T∈T_h

|v−vh|²_Hλ(T)+ Z

T

[(v−v_h)(x)]² ρ(x, ∂T)^2λ dx

. (2.6)

Hence the proof of (2.5) is reduced to a local estimate which can be handled by (1.8) and the usual scaling argument as follows.

LetKˆ be the reference element and consider the transformationxˆ 7→ x = BTxˆ+b, whereBT is an invertible matrix. By (1.2), (1.8), Friedrichs’ inequality and the regularity of T_h, we obtain

|v−vh|²_Hλ(T)+ Z

T

[(v−vh)(x)]² ρ(x, ∂T)^2λ dx

≤ kB_T⁻¹k^d+2λ|detB_T|²|ˆv−ˆv_h|²_{Hλ( ˆT)} +kB_T⁻¹k^2λ|detBT|

Z

Tˆ

[(ˆv−ˆvh)(ˆx)]² ρ(ˆx, ∂Tˆ)^2λ dˆx

≤CλkB⁻¹_T k^2λ|detBT| kˆv−vˆhk²_Hλ( ˆT)

≤CλkB⁻¹_T k^2λ|detBT| |ˆv|²_{Hλ( ˆT)}

≤C_λkB⁻¹_T k^2λ|detB_T| |ˆv|²_{Hµ( ˆT)}

≤C_λkB⁻¹_T k^2λ|detB_T|kB_Tk^d+2µ|detB⁻¹_T |²|v|²_Hµ(T)

≤Cλh^2(µ−λ)|v|²_Hµ(T). In view of (2.6), we deduce that

|v−vh|²_Hλ(Ω)≤Cλh^2(µ−λ) X

T∈T_h

|v|²_Hµ(T)≤Cλh^2(µ−λ)|v|²_Hµ(Ω),

and the lemma follows.

By a straight-forward argument (or cf. [13]), we also have

kv−Π⁰_hvkL²(Ω)≤C_λh^µ|v|H^µ(Ω) ∀v∈H^µ(Ω), λ≤µ≤1, which together with Theorem2.3imply

kv−Π⁰_hvk_H^λ(Ω)≤Cλh^µ−λ|v|H^µ(Ω) ∀v∈H^µ(Ω), λ≤µ≤1.

(2.7)

For the applications in Section 3 we will need interpolation error estimates in the dual spaces of the Sobolev spaces. The operatorΠ⁰_hcan be extended to the dual spaceH^λ(Ω)⁰for λ∈(0,¹₂)as follows:

hv−Π⁰_hv, ψhiλ,Ω = 0, ∀v∈H^λ(Ω)⁰, ψh∈Mh(Γ), (2.8)

whereh·,·iλ,Ω is the duality pairing betweenH^λ(Ω)⁰ andH^λ(Ω)that generalizes the inner product(·,·)_L²_(Ω).

THEOREM2.4. For anyλ∈(0,¹₂)andµ∈[λ,1], the following error estimate holds: kv−Π⁰_hvkH^µ(Ω)⁰ ≤Cλh^µ−λkvk_H^λ(Ω)⁰ ∀v∈H^λ(Ω)⁰.

Proof. Letv ∈ H^λ(Ω)⁰be arbitrary. By duality and the definitions (2.4) and (2.8), we can write

kv−Π⁰_hvk_H^µ(Ω)⁰ = sup

w∈H^µ(Ω)

hv−Π⁰_hv, wiµ,Ω

kwk_H^µ(Ω)

= sup

w∈H^µ(Ω)

hv, w−Π⁰_hwiλ,Ω

kwk_H^µ(Ω)

. (2.9)

From (2.7), we have

hv, w−Π⁰_hwiλ,Ω≤ kvk_H^λ_(Ω)⁰kw−Π⁰_hwk_H^λ_(Ω)≤Cλh^µ−λkvk_H^λ_(Ω)⁰kwk_H^µ(Ω), and the theorem follows from (2.9).

REMARK2.5. Theorem2.3and Theorem2.4remain valid on ad-dimensional(d= 1,2) polyhedral surfaceΓ, with essentially identical proofs. Note that the spaceH^λ(Γ)⁰can also be written asH^−λ(Γ).

Next, we consider the nodal interpolation operatorΠ¹_hfromC( ¯Ω)to aC⁰Lagrange (or tensorial) finite element spaceVhassociated withT_h(cf. [12], [8]).

THEOREM2.6. Letλ∈(0,¹₂). Then we have, ford= 1,2,

|v−Π¹_hv|_H^1+λ_(Ω)≤Cλh^µ−λ|v|_H^1+µ_(Ω) ∀v∈H^1+µ(Ω), λ≤µ≤1, (2.10)

and ford= 3,

|v−Π¹_hv|_H1+λ(Ω)≤Cλ,µh^µ−λ|v|_H^1+µ_(Ω) ∀v∈H^1+µ(Ω), 1

2 < µ≤1.

(2.11)

Proof. We will follow the notation introduced in Theorem2.3. First we consider the cases ofd= 1andd= 2. LetΠ¹_hv=vh. From Lemma2.1we have

|v−vh|²_H1+λ(Ω)≤C X

T∈T_h

|v−vh|²_H1+λ(T)+ Z

T

|∇v(x)− ∇vh(x)|² ρ(x, ∂T)^2λ dx

. (2.12)

Using scaling, (1.2), (1.8) and the regularity ofT_h, we have

|v−vh|²_H1+λ(T)+ Z

T

|∇v(x)− ∇vh(x)|² ρ(x, ∂T)^2λ dx

≤ kB_T⁻¹k^d+2λkBTk²|detBT|²|ˆv−vˆh|²_H1+λ( ˆT)

+kB_T⁻¹k^2λkBTk²|detBT| Z

Tˆ

|∇ˆv(x)− ∇ˆvh(x)|² ρ(ˆx, ∂T)ˆ^2λ dˆx (2.13)

≤CλkB⁻¹_T k^2λkBTk²|detBT|

|ˆv−vˆh|²_H1( ˆT)+|ˆv−ˆvh|²_H1+λ( ˆT)

. From the equivalence of norms on finite dimensional vector spaces and Sobolev’s in- equality (cf. [1], [19], [22]), it is easy to see that

|ˆv−vˆh|²_H1( ˆT)+|ˆv−vˆh|²_H1+λ( ˆT)≤Ch

|ˆv|²_H1( ˆT)+|ˆv|²_H1+λ( ˆT)+kˆvhk²_L∞

( ˆT)

i

≤Ch

|ˆv|²_H1( ˆT)+|ˆv|²_H1+λ( ˆT)+kˆvk²_L∞

( ˆT)

i (2.14)

≤C_λkˆvk²

H^1+λ( ˆT).

Combining (1.2), (2.13), (2.14), the Bramble-Hilbert lemma (cf. [6], [13]) and scaling, we have

|v−v_h|²_H1+λ(T)+ Z

T

|∇v(x)− ∇vh(x)|² ρ(x, ∂T)^2λ dx

≤C_λkB_T⁻¹k^2λkBTk²|detB_T|

× inf

p∈P1( ˆT)

h|(ˆv−p)−(ˆv−p)h|²_H1( ˆT)+|(ˆv−p)−(ˆv−p)h|²_{H1+λ( ˆT)}i

≤CλkB_T⁻¹k^2λkBTk²|detBT| inf

p∈P1( ˆT)

kˆv−pk²_{H1+λ( ˆT)} (2.15)

≤CλkB_T⁻¹k^2λkBTk²|detBT||ˆv|²_{H1+λ( ˆT)}

≤CλkB_T⁻¹k^2λkBTk²|detBT||ˆv|²_{H1+µ( ˆT)}

≤C_λkB_T⁻¹k^2λkB_Tk²|detB_T|kB_Tk^2d+2µkB⁻¹_T k²|detB⁻¹_T |²|v|²_Hµ(T)

≤Cλh^2(µ−λ)|v|²_Hµ(T).

The estimate (2.10) follows from (2.12) and (2.15).

The proof of (2.11) is similar, except that (2.14) must be replaced by

|ˆv−ˆvh|²_H1( ˆT)+|ˆv−vˆh|²_H1+λ( ˆT)≤Cµkˆvk²_H1+µ( ˆT).

The theorem now follows.

REMARK2.7. The results of [7] can be recovered from (2.10) by takingµ=λ.

Finally, we turn to inverse estimates involving fractional order Sobolev norms. For this purpose we will assume that the triangulationT_his quasi-uniform (cf. [12], [8]).

THEOREM2.8. Letλ∈(0,¹₂)andθ∈[0, λ]. Then the following estimate holds:

|v|_H^λ(Ω)≤C_λh^θ−λ|v|_H^θ(Ω) ∀v∈M_h(Ω).

Proof. We will follow the notation in Theorem2.3. Letv ∈Mh(Ω)be arbitrary. From Lemma2.1we have

|v|²_Hλ(Ω)≤C X

T∈T_h

|v|²_Hλ(T)+ Z

T

[v(x)]² ρ(x, ∂T)^2λ dx

. (2.16)

Using the equivalence of norms on a finite dimensional vector space, we obtain, by scaling, (1.2), (1.8) and the quasi-uniformity ofT_h,

|v|²_Hλ(T)+ Z

T

[v(x)]²

ρ(x, ∂T)^2λ dx≤ kB_T⁻¹k^d+2λ|detB_T|²|ˆv|²_{Hλ( ˆT)} +kB⁻¹_T k^2λ|detB_T|

Z

Tˆ

[ˆv(ˆx)]² ρ(ˆx, ∂T)ˆ ^2λ dˆx

≤CλkB⁻¹_T k^2λ|detBT|kˆvk²_L2( ˆT)

≤C_λkB⁻¹_T k^2λ|detB_T||detB⁻¹_T |kvk²_L2(T)

≤C_λh^−2λkvk²_L2(T). Combining the preceding estimate and (2.16) we have

|v|_H^λ_(Ω)≤Cλh^−λkvk_L²(Ω) ∀v∈Mh(Ω).

(2.17)

In other words the theorem holds forθ= 0.

The proof forθ ∈ (0, λ]is more complicated. For T ∈ T_h, we denote byσ_T the collection of elements inT_hwhich share at least one vertex withT, i.e.,

σT ={T⁰∈T_h: T ∩T⁰6=∅}.

The domainST is defined by

ST = [

T⁰∈σT

T⁰. (2.18)

We have

|v|²_Hλ(Ω)= X

T,T⁰∈T_h

Z

T

Z

T⁰

[v(x)−v(y)]²

|x−y|^d+2λ dxdy

= X

T,T⁰∈T_h T⁰∈σT

Z

T

Z

T⁰

[v(x)−v(y)]²

|x−y|^d+2λ dxdy (2.19)

+ X

T,T⁰∈T_h T⁰6∈σT

Z

T

Z

T⁰

[v(x)−v(y)]²

|x−y|^d+2λ dxdy.

There is an easy estimate for the second sum on the right-hand side of (2.19):

X

T,T⁰∈T_h T⁰6∈σT

Z

T

Z

T⁰

[v(x)−v(y)]²

|x−y|^d+2λ dxdy

≤Ch^2(θ−λ) X

T,T⁰∈T_h T⁰6∈σT

Z

T

Z

T⁰

[v(x)−v(y)]²

|x−y|^d+2θ dxdy (2.20)

≤Ch^2(θ−λ)|v|²_Hθ(Ω).

In view of (1.2) and (2.18), we can bound the first sum on the right-hand side of (2.19) by

X

T,T⁰∈T_h T⁰∈σT

Z

T

Z

T⁰

[v(x)−v(y)]²

|x−y|^d+2λ dxdy≤ X

T∈T_h

|v|²_Hλ(S_T). (2.21)

LetSˆ_T be a reference domain with unit diameter which is similar toS_T,F be the affine transformation that mapsSˆT toST, andˆv=v◦F be the pull-back ofvtoSˆT. We obtain, by applying (2.17) toSˆT and using the Bramble-Hilbert lemma (cf. [13]),

|ˆv|_Hλ( ˆST)= inf

p∈P0( ˆST)

|ˆv−p|_Hλ( ˆST)

(2.22)

≤C_λ inf

p∈P0( ˆST)

kˆv−pk_L2( ˆST)≤C_λ|ˆv|_Hθ( ˆST). Combining (2.22) with a scaling argument, we have

|v|_H^λ(ST)≤C_λh^θ−λ|v|_H^θ(ST), which together with (2.21) imply

X

T,T⁰∈T_h T⁰∈σT

Z

T

Z

T⁰

[v(x)−v(y)]²

|x−y|^d+2λ dxdy≤C_λh^2(θ−λ) X

T∈T_h

|v|²_Hθ(ST)

(2.23)

≤C_λh^2(θ−λ)|v|²_Hθ(Ω). The case forθ∈(0, λ]now follows from (2.19), (2.20) and (2.23).

In exactly the same way one can prove the following theorem for aC⁰ Lagrange (or tensorial) finite element spaceV_h(Ω).

THEOREM2.9. Letλ∈(0,¹₂)andθ∈[0, λ]. Then the following estimate holds:

|v|_H^1+λ_(Ω)≤Cλh^θ−λ|v|_H^1+θ_(Ω) ∀v∈Vh(Ω).

3. Bubble stabilization of three-dimensional mixed finite element methods. In this section we apply Theorem2.4to establish optimal error estimates for the numerical solution of some second order elliptic partial differential equations by the bubble stabilized finite el- ement method, when the exact solution is not regular. First, we present the discretization of the 3D Poisson problem, where the Dirichlet condition is dualized `a la Babuˇska. The second application deals with the approximation of the 3D unilateral contact problem known as the Signorini problem.

3.1. The Poisson problem with Lagrange multipliers. LetΩbe a bounded polyhedral domain inR³with boundaryΓ =∂Ω. Givenf ∈L²(Ω)andg ∈H¹²(Γ), the problem of interest is the Poisson equation:

−∆u=f inΩ, (3.1)

u=g onΓ.

(3.2)

By using Lagrange multipliers to enforce the Dirichlet condition (cf. [2]), we can formu- late (3.1)–(3.2) as the following saddle point problem:

Find(u, ϕ)∈H¹(Ω)×H¹²(Γ)⁰such that Z

Ω

∇u· ∇v dx+hϕ, vi¹

2,Γ= Z

Ω

f v dx ∀v∈H¹(Ω), (3.3)

hψ, ui¹

2,Γ=hψ, gi¹

2,Γ ∀ψ∈H¹²(Γ)⁰. (3.4)

LetT^Ω

h be a regular triangulation made of elements that are tetrahedra with a maximum sizeh(the extension to the hexahedra does not create any technical difficulty). The trace of T^Ω

h on the boundaryΓresults in a regular (2D) triangulationT^Γ

h made of triangular elements which are the faces of tetrahedral elements inT^Ω

h . LetYh(Ω)be the space defined by Yh(Ω) =n

vh∈C⁰(Ω) :vh

κ∈ P1(κ) ∀κ∈T^Ω)

h

o,

whereP1(κ)is the set of all affine functions overκ. In the stabilized finite element approach, the discrete spaceY_h(Ω)is enriched by cubic bubble functions defined on the boundaryΓ.

Let{x1, x2, x3}be the vertices of the triangleT ∈ T^Γ

h which is a face of the tetrahedral elementκ_T inT^Ω

h . The vertices ofκ_T are(xi)1≤i≤4andλ_i (1≤i≤4)is the barycentric coordinate associated withxi. The bubble function we need to use is defined to be

ϕ_T(x) = 60

|T|λ1(x)λ2(x)λ3(x) ∀x∈κ_T,

and extended by zero elsewhere. Then the locally stabilized finite element space is given by Xh(Ω) =Yh(Ω)⊕ M

T∈T^Γ

h

RϕT

.

The approximate Lagrange multipliers are piecewise constant functions with respect to the meshT^Γ

h , i.e.,

M_h(Γ) =n

ψ_h∈L²(Γ) :ψ_h

T ∈ P0(T) ∀T ∈T^Γ

h

o, and the discrete problem for (3.3)–(3.4) is:

Find(uh, ϕh)∈Xh(Ω)×Mh(Γ)such that Z

Ω

∇u_h· ∇v_hdx+hϕ_h, v_hi¹

2,Γ= Z

Ω

f v_hdx ∀v_h∈X_h(Ω), (3.5)

hψh, uhi¹

2,Γ=hψh, gi¹

2,Γ ∀ψh∈Mh(Γ).

(3.6)

The properties that allow for existence and uniqueness results are the coercivity of the form(uh, vh)7→ (∇(·),∇(·))_L²(Ω)on a subspace ofXh(Ω)and the so-calledinf-supcon- dition of the form (vh, ψh) 7→ hψh, vhi¹

2,Γ on the spacesXh(Ω) ×Mh(Γ). There is no

particular difficulty in checking that the seminorm| · |_H¹(Ω)is equivalent to theH¹-norm in the space

nv_h∈X_h(Ω) :hψ_h, v_hi¹

2,Γ= 0 ∀ψ_h∈M_h(Γ)o .

This proves the coercivity of the first form on this space. Moreover, following the lines of [3], the spacesX_h(Ω)andM_h(Γ)satisfy the Babuˇska-Brezzi condition (also known as the inf-supcondition):

ψh∈Minfh(Γ) sup

vh∈Xh(Ω)

hψh, vhi¹

2,Γ

kvhk_H¹_(Ω)kψhk

H¹²(Γ)⁰

> γ,

where the constantγdoes not depend onh. Therefore, using the saddle point theory (cf. [9]), one can prove that problem (3.5)–(3.6) has a unique solution(uh, ϕh)∈ Xh(Ω)×Mh(Γ) which satisfies the following abstract error estimate:

ku−u_hk_H¹_(Ω)+kϕ−ϕ_hk

H¹²(Γ)⁰ ≤ (3.7)

C

vh∈Xinfh(Ω)ku−vhk_H¹(Ω)+ inf

ψh∈Mh(Γ)kϕ−ψhk

H¹²(Γ)⁰

.

In the following theorem we apply Theorem2.4to derive from (3.7) the convergence rate of our mixed finite element solution of (3.5)–(3.6).

THEOREM 3.1. Assume that the exact solution u of Poisson’s problem belongs to H^1+λ(Ω) (0< λ < ¹₂). Then the following error estimate holds:

ku−uhkH¹(Ω)+kϕ−ϕhk

H¹²(Γ)⁰ ≤Ch^λ(kuk_H^1+λ(Ω)+kfkL²(Ω)).

Proof. It is sufficient to observe that, sinceu∈H^1+λ(Ω)and∆u∈L²(Ω),ϕ= ^∂u_∂n ∈ H^12−λ(Γ)⁰with

kϕkH^12−λ(Γ)⁰ ≤C(kuk_H^1+λ_(Ω)+kfk_L²(Ω)).

We can then use Theorem2.4to obtain

ψh∈Minfh(Γ)kϕ−ψhk

H¹²(Γ)⁰ ≤Ch^λkϕk

H^12−λ(Γ)⁰ ≤Ch^λ(kuk_H^1+λ_(Ω)+kfk_L²_(Ω)).

The bound

vh∈Xinf_h(Ω)ku−v_hk_H¹_(Ω)≤Ch^λkuk_H1+λ(Ω)

can be obtained by using a finite element interpolation for non-smooth functions (cf. [21], [5]).

REMARK 3.2. In two dimensions, there is no need for stabilization since the natural spacesY_h(Ω)andM_h(Γ) (built as described above with obvious modification)are compat- ible regarding theinf-supcondition(cf. [2], [20]). In three dimensions this condition is lost and is restored by bubble stabilization (cf. [10], [11]). Note, however, that even for non- regular two-dimensional solutions, the one-dimensional result in Theorem2.4is needed for proving optimal convergence results.

3.2. The unilateral contact variational inequality. Assume again thatΩis a bounded polyhedral domain inR³. The boundaryΓ =∂Ωis a union of three non-overlapping sections Γu,ΓgandΓC. The partΓuof nonzero measure is subject to the Dirichlet conditions, while onΓgthe Neumann condition is prescribed, andΓC is the candidate to be in contact with a rigid obstacle. To avoid technicalities arising from the special Sobolev spaceH

1 2

00(ΓC), we assume thatΓuandΓCdo not touch.

For given dataf ∈L²(Ω)andg∈H¹²(Γg)⁰, the Signorini problem consists of finding usuch that

−∆u=f inΩ, (3.8)

u= 0 onΓu, (3.9)

∂u

∂n =g onΓg, (3.10)

u≥0, ∂u

∂n ≥0, u∂u

∂n = 0 onΓC, (3.11)

wherenis the outward unit normal of∂Ω. This model is currently encountered in the condi- tioning field (whereuis a temperature) and in the hydrostatic domain (whereuis a pressure).

Sometimes, for practical reasons, one may want to use the mixed formulation where the conditionϕ = ^∂u_∂n ≥ 0is taken into account explicitly. In this approach the space for the displacementuis the subspaceH₀¹(Ω,Γu)ofH¹(Ω)consisting of functions that vanish on Γu, and the flux (normal derivative)ϕ= ^∂u_∂n

ΓC belongs to the closed convex set M(ΓC) =n

ψ∈H¹²(ΓC)⁰ :ψ≥0o .

Here the nonnegativity of a distribution ψ ∈ H¹²(ΓC)⁰ is to be understood in the sense thathψ, χi¹

2,ΓC ≥ 0for any nonnegativeχ ∈H¹²(ΓC). Then(u, ϕ)is the solution of the following mixed variational inequality:

Find(u, ϕ)∈H₀¹(Ω,Γu)×M(ΓC)such that Z

Ω

∇u· ∇v dx+hϕ, vi¹

2,ΓC = Z

Ω

f v dx ∀v∈H₀¹(Ω,Γu), (3.12)

hψ−ϕ, ui¹

2,ΓC ≥0 ∀ψ∈M(ΓC).

(3.13)

A complete analysis of this mixed problem is provided in [17] (cf. also [16]) where an existence and uniqueness result is proven (cf. Theorem 3.14 therein). Moreover we have the following stability estimate:

kuk_H¹(Ω)+kϕk

H¹²(ΓC)⁰ ≤C(kfk_L²(Ω)+kgk

H¹²(Γg)⁰).

It is useful to note that ifK(Ω)is the convex cone K(Ω) =n

v∈H₀¹(Ω,Γu) :v≥0o ,

thenu∈K(Ω)is also the unique solution of the following primal problem:

Findu∈K(Ω)such that Z

Ω

∇u· ∇(v−u)dx≥ Z

Ω

f(v−u)dx+hg, v−ui¹

2,Γg ∀v∈K(Ω).

(3.14)

REMARK3.3. The Neumann and Signorini conditions (3.10) and (3.11) are taken into account in a weak sense in the primal variational inequality (3.14). Indeed, we have

h∂u

∂n, vi¹

2,Γ− hg, vi¹

2,Γg ≥ 0 ∀v∈H

1 2

00(Γ,Γu) and v _Γ

C ≥0, (3.15)

h∂u

∂n, ui¹

2,Γ− hg, ui¹

2,Γg = 0, (3.16)

where H

1 2

00(Γ,Γu) is the subspace of H¹²(Γ) consisting of functions that vanish on Γu. Roughly speaking, (3.15) says that ^∂u_∂n =gonΓgand_∂n^∂u ≥0onΓCwhile (3.16) expresses the saturation conditionu^∂u_∂n = 0onΓC.

The approximation of the variational inequality (3.12)–(3.13) is obtained by generalizing the two-dimensional bubble stabilized mixed finite elements in [3] to three dimensions. The finite element tools are the same as those introduced in the previous section. We assume moreover that the meshT^Ω

h is compatible with the partition of the boundary, meaning that the trace of it toΓu,Γgand toΓCresults in two dimensional triangulations. The triangulation ofΓCis denoted byT^C

h .

Taking into account the Dirichlet boundary condition, we define Yh(Ω) =n

vh∈C⁰(Ω) :vh

κ∈ P1(κ) ∀κ∈T^Ω

h and vh

ΓC = 0o .

The finite element space whereuhis computed is then given by Xh(Ω) =Yh(Ω)⊕ M

T∈T^C

h

RϕT

,

and the convex cone for the discrete Lagrange multipliers onΓCis taken to be Mh(ΓC) =n

ψh∈L²(ΓC) :ψh

T ∈ P0(T) ∀T ∈T^C

h and ψh≥0o . We are now ready to set the discrete mixed variational inequality:

Find(uh, ϕh)∈Xh(Ω)×Mh(ΓC)such that Z

Ω

∇u_h· ∇v_hdx+hϕ_h, v_hi¹

2,ΓC = Z

Ω

f v_hdx ∀v_h∈X_h(Ω), (3.17)

hψh−ϕh, uhi¹

2,ΓC ≥0 ∀ψh∈Mh(ΓC).

(3.18)

Again the availability of the Babuˇska-Brezzi condition

ψh∈Minfh(ΓC) sup

vh∈Xh(Ω)

hψh, vhi¹

2,ΓC

kvhk_H¹(Ω)kψhk

H¹²(ΓC)⁰

> γ

allows us to prove existence, uniqueness and stability results.

The analysis of the discretization error is based on the saddle point theory for variational inequalities (cf. [17], [16]). The methodology is to first obtain the convergence rate on the primal variableuby analyzing (3.14) and its approximation. We have therefore to write down a variational problem foruhafter suppressing the Lagrange multiplier. Let us then introduce the closed convex cone

K_h(Ω) =n

v_h∈X_h(Ω) :hψ_h, v_hi¹

2,ΓC = Z

ΓC

ψ_hv_hdΓ≥0 ∀ψ_h∈M_h(ΓC)o .

It is easy to see thatuhis also the unique solution of the following variational inequality:

Finduh∈Kh(Ω)such that Z

Ω

∇u_h· ∇(v_h−u_h)dx≥ Z

Ω

f(vh−u_h)dx+hg, v_h−u_hi¹

2,Γg

(3.19)

∀v_h∈K_h(Ω).

Note thatK_h(Ω) 6⊂K(Ω)and hence problem (3.19) is a nonconforming discretization of (3.14). Applying Falk’s lemma (cf. [14]) to (3.19) yields the following error estimate, the proof of which can be found for instance in [4].

LEMMA3.4. The following error estimate holds:

ku−u_hk²_H1(Ω)≤Ch

vh∈Kinfh(Ω) ku−v_hk²_H1(Ω)+< ∂u

∂n, v_h>¹

2,Γ−hg, v_hi¹

2,Γg

(3.20)

+ inf

v∈K(Ω) h∂u

∂n, v−uhi¹

2,Γ− hg, v−uhi¹

2,Γg

i .

The first infimum of the bound given in Lemma 3.2 is the approximation error and the boundary term involved there is specifically generated by the discretization of variational inequalities. The second infimum is the consistency error, the price for the “variational crime”

due to the nonconformity of the approximation. These two errors will be studied separately.

LEMMA 3.5. Assume that for someλ(0 < λ < ¹₂)we haveg∈ H^12−λ(Γg)⁰and that the exact solutionuof Signorini’s problem belongs toH^1+λ(Ω). Then the following estimate holds:

vh∈Kinfh(Ω) ku−vhk²_H¹_(Ω)+< ∂u

∂n, vh>¹

2,Γ−hg, vhi¹

2,Γg

≤ Ch^2λh

kuk_H1+λ(Ω)+kfk_L²_(Ω)+kgk

H^12−λ(Γg)⁰

ikuk_H1+λ(Ω).

Proof. It suffices to choosevh∈Yh(Ω)to be the Lagrange interpolant ofu. It is checked thatvh∈Kh(Ω)sincevh

ΓC ≥0. The estimate is derived following [4].

LEMMA 3.6. Assume that for someλ(0 < λ < ¹₂)we haveg∈ H^12−λ(Γg)⁰and that the exact solutionuof Signorini’s problem belongs toH^1+λ(Ω). Then the following estimate holds:

v∈K(Ω)inf h∂u

∂n, v−uhi¹

2,Γ− hg, v−uhi¹

2,Γg

≤ Ch^λh

kuk_H^1+λ_(Ω)+kfk_L²(Ω)+kgk

H^12−λ(Γg)⁰

i ku−uhk_H¹(Ω)+h^λkuk_H^1+λ_(Ω) .

Proof. First of all, observe that sinceg ∈ H^12−λ(Γg)⁰ and ^∂u_∂n ∈H^12−λ(Γ)⁰, a density argument and (3.15) imply

h∂u

∂n, vi¹

2−λ,Γ− hg, vi¹

2−λ,Γg ≥0 ∀v∈H^12−λ(Γ) and v _Γ

C ≥0.

It follows that ifψ_h = Π⁰_h(^∂u_∂n)∈ M_h(Γ)is the piecewise constant interpolant of ^∂u_∂n, then we have

ψ_h T = 1

|T|h∂u

∂n,1Ti¹

2−λ,Γ = 1

|T|hg,1Ti¹

2−λ,Γg ∀T ∈T^C

h ,