This paper on the a posteriori error analysis of the obstacle problem with affine obstacles and Courant finite elements compares five classes of error estimates for accurate guaranteed error control

FOR CONFORMING OBSTACLE PROBLEMS^?

C. CARSTENSEN AND C. MERDON

Abstract. This paper on the a posteriori error analysis of the obstacle problem with affine obstacles and Courant finite elements compares five classes of error estimates for accurate guaranteed error control. In order to treat interesting computational benchmarks, the first part extends the Braess methodology from 2005 of the resulting a posteriori error control to mixed inhomogeneous boundary conditions. The resulting guaranteed global upper bound (GUB) involves some auxiliary partial differential equation and leads to four contributions with explicit constants. Their efficiency is examined affirmatively for five benchmark examples.

1. Introduction

The a posteriori error analysis is well developed for variational equations of second-order elliptic partial differential equations with explicit constants and at least five different types of error estimators. Considerably less is known about the a posteriori error control for vari- ational inequalities, in particular, for its model obstacle problem [CHL, BCH07, BC04a, NSV05,NSV03a,Bra05,BHS08,Vee01a]. Braess’ error estimator split for an obstacle prob- lem [Bra05] leads to some computable term with an extended discrete Lagrange multiplier plus the Galerkin discretisation error of an auxiliary variational equation.

Given a bounded Lipschitz domain Ω⊂R² with boundary∂Ω and its closed subset Γ_D of positive surface measure, the data of the obstacle problem are the right-hand sidesf ∈H¹(Ω) and g∈L²(Γ_N) on Γ_N :=∂Ω\Γ_D plus the Dirichlet datau_D ∈C⁰(Γ_D) edgewise inH² and the obstacleχinH¹(Ω)∩C⁰(Ω) withχ≤u_D along Γ_D. The weak form relies on the closed and convex set

K:={v∈H¹(Ω)

v=u_D on Γ_D and χ≤v a.e. in Ω} 6=∅. The obstacle problem reads: Seeku∈K such that

a(u, u−v)≤F(u−v) for all v∈K.

(1.1)

Here and throughout the paper, we use standard notation for Lebesgue and Sobolev spaces and denote the energy norm

||| · |||:=a(·,·)^1/2 with respect to the bilinear formain H¹(Ω),

a(u, v) :=

Z

Ω

∇u· ∇vdx (1.2)

Key words and phrases. adaptive finite element methods, a posteriori error estimators, elliptic variational inequalities, obstacle problems.

?This work was supported by DFG Research Center MATHEON and by the World Class University (WCU) program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education, Science and Technology R31-2008-000-10049-0.

1

for allu∈u_D+V andv∈V :={v∈H¹(Ω)

v|_ΓD = 0}. The right-hand sideF in the dual V^? of V reads

F(v) :=

Z

Ω

f vdx+ Z

ΓN

gvds for all v∈V.

(1.3)

It is well known that (1.1) has a unique solution [KS80,Zei97,Rod87].

Given a regular triangulation T of Ω in the sense of Ciarlet [BS94,Cia78] with nodes N and nodal basis ϕ_z

z∈ N

of the first-order finite element method, the nodal interpolation of v∈C⁰(Ω) reads

Iv= X

z∈N

v(z)ϕ_z ∈P₁(T)∩C⁰(Ω).

The discrete version of problem (1.1) employs the discrete set K(T) :={v_h ∈P₁(T)∩C⁰(Ω)

v_h =Iu_D on Γ_D and Iχ≤v_h in Ω}

(with the piecewise affine functionsP₁(T)) and reads: Seek u_h ∈K(T) such that a(u_h, u_h−v_h)≤F(u_h−v_h) for all v_h ∈K(T).

(1.4)

The finite element method is calledconforming ifK(T)⊂K, for instance when χ =Iχ equals its nodal interpolation Iχ, and is called non-conforming otherwise.

After [Bra05] and for a particular choice of Λh, the discrete solution of the obstacle problem u_h solves the discrete version of the Poisson problem forw∈ Iu_D+V with

a(w, v) =F(v)− Z

Ω

Λhvdx for all v∈V.

(1.5)

The energy norm difference |||w−uh||| between uh and the exact solution w of the Poisson problem (1.5) can be estimated by a large collection of error estimators [AO00,BR96,BS01, EEHJ96,Ver96,Rep08,RSS03,LW04,Bra07], in particular the ones compared in [BCK01, CM10].

This paper extends the result in [Bra05] to problems with inhomogeneous Dirichlet bound- ary data. Suppose wD ∈ H¹(Ω) satisfies wD|_ΓD = uD − Iu_D and χ−uh ≤ wD. In the conforming case χ ≤ Iχ, our main result (Theorem 3.2) on w from (1.5) implies some computable global upper bound (GUB) slightly sharper than

|||e||| ≤ |||w−uh|||+|||Λ_h−JΛh|||∗+ Z

Ω

(χ−uh−wD)JΛhdx 1/2

+ 2|||w_D|||.

(1.6)

This bounds the errore:=u−u_h in the energy norm by the error |||w−u_h|||. The first extra term

|||Λ_h−JΛ_h|||_∗:= sup

v∈V\{0}

Z

Ω

(Λ_h−JΛ_h)vdx/|||v|||

relates to a non-positive approximationJΛhof Λh, whileR

Ω(χ−uh−wD)JΛhdx contributes only on the transition between the contact zone and the non-contact zone. The term w_D refers to inhomogeneous Dirichlet data approximation ands its minimal H¹(Ω) extension.

For affine obstacles, Theorem4.1proves efficiency in the sense of GUB.|||e|||+|||Λ−Λ_h|||_∗

up to pertubations that are of higher order, at least for benchmarks with ΓN := ∅. The Lagrange multiplier Λ is defined below in (3.10).

The remaining parts of this paper are organized as follows. The interpolation operator J and a suitable choice for wD are described in Section 2. Section 3 presents the details of our construction of Λ_h and discusses reliability of the above estimate. Section 4 discusses its efficiency in case of affine obstacles. Section 5 describes an adaptive mesh refinement algo- rithm and six known a posteriori error estimators from Poisson problems for the estimation of |||w−u_h||| in the global upper bound. Section 6 discusses five computational benchmark examples. Two affine problems verify the theoretical findings and suggest that the overhead apart from the estimator contributions is not dominant.

2. Preliminaries

Consider a regular triangulationT of Ω⊆R²with nodesN, free nodesM:=N \Γ_D, fixed nodes N(ΓD) :=N \ M, edges E, Dirichlet edges E(Γ_D) := {E ∈ E

E ⊆ ΓD}, Neumann edges E(Γ_N) := {E ∈ E

E ⊆ΓN} and interior edges E(Ω) :=E \(E(Γ_D)∪ E(Γ_N)). Each node z inN is associated with its nodal basis functions ϕ_z and node patch ω_z :={ϕ_z >0}

with diameter h_z := diam(ω_z). The quantity hE ∈ P₀(E) denotes the local edge length, i.e., hE|_E := hE := |E| for E ∈ E. Similarly, hT denotes the local triangle diameter, hT|_T :=h_T := diam(T) forT ∈ T. Each element T ∈ T is the convex hull of the setN(T) of its three vertices and associated to its element patch ωT := S

z∈N(T)ωz. The set E(T) denotes the edges on the boundary of T. For two expressions A and B, we write ”A. B”

if B ≤ CA for some generic constant C that depends only on the shape regularity of the triangulation.

Given anyv∈H¹(Ω), setJ v∈P₁(T)∩C⁰(Ω) similar to [BC04a] by J v:= X

z∈N

vzϕz with vz :=

Z

Ω

vϕzdx/ Z

Ω

ϕzdx ∈R.

Givenf ∈L²(Ω), let osc(f,N) :=

P

z∈Nh²_zmin_fz∈Rkf−f_zk²_L2(ω

z)

1/2

. Lemma 2.1. Forf ∈L²(Ω) and v∈H¹(Ω), it holds

Z

Ω

(J f)vdx = Z

Ω

f(J v)dx and Z

Ω

f(v−J v)dx .osc(f,N)|||v|||.

Proof. The first assertion follows from a direct calculation as in [BC04a]. For the proof of

the second assertion, we refer to [Car99,CV99].

Theorem 2.2 ([BCD04]). Assume that uD ∈ H¹(ΓD)∩C⁰(ΓD) satisfies uD ∈ H²(E) for all E ∈ E(Γ_D) with the edgewise second partial derivative ∂_E²u_D/∂s² of u_D along Γ_D and let Iu_D =P

z∈N(Γ_D)uD(z)ϕz denote the nodal interpolation of uD. Then there exists wD ∈H¹(Ω) with

wD|_ΓD =uD|_ΓD − Iu_D|_ΓD, supp(wD)⊂[

{T ∈ T

T∩ΓD 6=∅}, kw_Dk_L∞(Ω) =ku_D − Iu_Dk_L∞(ΓD),

|||w_D|||.kh^3/2_E ∂_E²uD/∂s²k_L2(Γ_D). Furthermore, if χ≤ Iχ≤u_D it holds

χ−u_h ≤w_D.

Proof. The inequalityχ−u_h ≤w_D is an immediate consequence of the design in [BCD04]

withχ−u_h≤ Iχ−uh ≤uD−IuD =wD on E⊂ΓD and (Iχ) (mid(T))−uh(mid(T))≤0 = w_D(mid(T)) plus the linearity ofw_D andIχ−u_h along the lines that connect the boundary points on E with mid(T). The remaining details are contained in [BCD04].

3. Reliable Error Estimation

This section is devoted to the design of reliable error estimators for the Poisson problem (1.5) with some function Λh. Subsections 3.1-3.3 introduce Λh in three steps from a post- processed Riesz representation of a residual similar to [Bra05,BC04a].

3.1. Residuals. The exact solution u of (1.1) defines the residual functional σ :=F−a(u,·)∈V^∗ := dual of V.

(3.1)

The discrete solution u_h of (1.4) defines the discrete residual functional σh :=F−a(uh,·)∈V(T)^∗

(3.2)

with the test function space V(T) := P1(T)∩V. The two properties (3.3) and (3.4) are important in the sequel. The discrete complementary conditions read

(u_h(z)− Iχ(z))σ_h(ϕ_z) = 0 and σ_h(ϕ_z)≤0 for allz∈ M.

(3.3)

ForwD ∈H¹(Ω) withwD =uD − Iu_D on ΓD and χ−uh≤wD in Ω, it holds 0≤σ(e−wD).

(3.4)

3.2. Boundary modifications. The complementary conditions (3.3) state thatσ_h ∈V(T)^∗ contains information about the contact zone {u_h =Iχ}. That information will be used to define some contact force on the entire domain. However, we have to define σh also at Dirichlet nodes. Veeser mentioned thatσ_h= 0 on Γ_D may lead to big refinement indicators [Vee01a, p. 153 line 25]. Our own numerical experiments confirm this observation. Here we avoid to enforceσ_h|_ΓD = 0 by extrapolation of the known contact information as follows.

For each fixed node z ∈ N \M, we choose a neighbouring free node ζ(z) ∈ M and set ζ(z) :=z forz ∈ M. In Section 4 we assume thatζ(y)∈ ωT for all y ∈ N and that every triangle has a free node. The proof of efficiency will make use of this further assumption, but it is not needed for the proof of reliability at the end of Section 3.

The mapping ζ defines a partition ofN into card(M) many classes ζ⁻¹(z) :={y∈ N

ζ(y) =z} for each z∈ M.

For eachz∈ M set

ψ_z := X

y∈ζ⁻¹(z)

ϕ_y ∈P₁(T)∩C⁰(Ω).

(3.5)

Notice that ψz

z∈ M

is a partition of unity [Car99]. For each z∈ N set

bσ_h(ϕz) :=

(0 ifχ(z)< u_h(z), σ_h(ϕ_ζ(z))

R

Ωϕzdx R

Ωϕ_ζ(z)dx else.

(3.6)

A direct consequence of (3.3) is

X

z∈N

σb_h(ϕz) R

Ωϕ_zdxϕz ≤0.

(3.7)

Remark 3.1. Veeser [Vee01a] suggests the alternative boundary modification with

bσ_h(ϕz) :=

(0 ifχ(z)< uh(z), min

n

0, F(ϕz)−a(uh, ϕz) +R

ΓDϕz∇u_h·νds o

else.

This resembles the original definition of the residual with an additional boundary term where ν denotes the outer unit normal along∂Ω.

3.3. Riesz representation and auxiliary problem. The Riesz representation Λ_h ∈P₁(T)

∩C⁰(Ω) of the extendedσbh∈ P1(T)∩C⁰(Ω)∗

from (3.6) in the Hilbert spaceL²(Ω) sat- isfies

Z

Ω

Λ_hϕ_zdx =bσ_h(ϕ_z) for all z∈ N. (3.8)

This defines the auxiliary problem (1.5) from the introduction: Seekw∈ Iu_D+V with a(w, v) =F(v)−

Z

Ω

Λ_hvdx for all v∈V.

(3.9)

In view of (3.8), the discrete solution u_h of (1.4) equals the finite element approximation of the exact solution w of (3.9). The residual (3.1) of the exact problem also defines some Riesz representation Λ∈L²(Ω) with

Z

Ω

Λvdx =σ(v) for allv∈V.

(3.10)

3.4. Main result. The view of Braess onto the obstacle problem at hand leads to the following reliable error estimate

|||u−u_h||| ≤ |||w−u_h|||+|||Λ_h−JΛ_h|||_∗+ Z

Ω

(χ−u_h−w_D)JΛ_hdx 1/2

+ 2|||w_D|||

The following result is slightly sharper and employs

a:=|||w−uh|||+|||Λ_h−JΛh|||_∗+|||w_D||| and b:=

Z

Ω

(χ−uh−wD)JΛhdx . Theorem 3.2. Let u denote the exact solution of (1.1) and let u_h denote the discrete solution of (1.4) with the associated Λh from (3.8) and Λ from (3.10). For χ ≤ Iχ and w_D ∈H¹(Ω) with w_D =u_D− Iu_D on Γ_D and χ−u_h ≤w_D in Ω, it holds

|||u−u_h||| ≤a/2 +p

a²/4 +b+|||w_D|||,

|||Λ−Λh|||∗ ≤ |||w−uh|||+|||u−uh|||.

Some remarks are in order before the proof of Theorem3.2 concludes this section.

Remark 3.3. The energy norm|||w−u_h|||of the discretisation error in the auxiliary problem equals the dual norm

|||Res|||_∗ := sup

v∈V

|||v|||=1

Res(v) of the residual Res∈V^∗ defined, for v∈V, by

Res(v) :=

Z

Ω

(f−Λh)vdx + Z

Γ_N

gvds−a(u_h, v).

(The Riesz representation theorem in the Hilbert space (V, a) yields direct proof, c.f. [BC04b]

for details.) The substitute of |||w−uh|||=|||Res|||∗ and Vh ⊆ ker(Res) allows an interpre- tation of the a posteriori error control in Theorem 3.2in the spirit of the unified approach [Car05]. In fact, all error estimators of [Car05,CEHL] apply to the control of|||Res|||∗. Remark 3.4. The termR

Ω(χ−uh−w_D)JΛhdx can be evaluated exactly (up to quadrature errors). In case χ=Iχ, this term only contributes on a layer between the discrete contact zone and the discrete non-contact zone{T ∈ T

∃z, y∈ N(T), χ(z)< uh(z) &bσh(ϕy)<0}.

Remark 3.5. The properties of wD from Theorem 2.2 guarantee that the term |||w_D||| is bounded by the higher-order term

|||w_D|||.kh^3/2_E ∂_E²u_D/∂s²k_L2(ΓD).

For triangulations that consists only of right isosceles triangles along the Dirichlet boundary, the constant in this estimate equals one. This will be proven in [CMx1]. Furthermore wD

contributes to R

Ω(χ−u_h−w_D)JΛ_hdx only on elements with contact near the boundary.

For homogeneous Dirichlet data, wD ≡0.

Remark 3.6. Theorem3.2 implies the reliable upper bound

|||Λ−Λ_h|||_∗ ≤2|||w−u_h|||+|||Λ_h−JΛ_h|||_∗+ Z

Ω

(χ−u_h−w_D)JΛ_hdx 1/2

+ 2|||w_D|||.

Remark 3.7. As a conclusion, for the estimation of|||u−u_h||| or|||Λ−Λ_h|||_∗ it remains to estimate |||w−u_h|||. This can be done with any a posteriori estimator known for Poisson problems, collected in [BCK01,CM10].

Proof of Theorem 3.2. The proof is similar to the proof of Lemma 3.1 in [Bra05] for σ_h⁺ :=

−JΛ_h. Indeed, forv ∈V, (3.1) and (3.9) imply a(u−w, v) =

Z

Ω

vΛ_hdx−σ(v)

= Z

Ω

vJΛhdx−σ(v) + Z

Ω

v(Λh−JΛh)dx. Forv:=u−u_h−w_D =e−w_D ∈V, it holds

Z

Ω

(u−u_h−wD)JΛ_hdx−σ(e−wD)

= Z

Ω

(χ−u_h−w_D)JΛ_hdx− Z

Ω

(χ−u)JΛ_hdx −σ(e−w_D).

The properties (3.4),χ−u≤0, andJΛh≤0 from (3.7) yield 0≤

Z

Ω

(χ−u)JΛ_hdx and 0≤σ(e−w_D).

Hence

Z

Ω

(u−u_h−w_D)JΛ_hdx−σ(e−w_D)≤ Z

Ω

(χ−u_h−w_D)JΛ_hdx. The combination of the previous equality and inequality with some algebra leads to

|||u−u_h−w_D|||²

=a(u−w, u−u_h−w_D) +a(w−u_h, u−u_h−w_D)−a(w_D, u−u_h−w_D)

≤

|||w−u_h|||+|||Λ_h−JΛ_h|||_∗+|||w_D|||

|||u−u_h−w_D|||+ Z

Ω

(χ−u_h−w_D)JΛ_hdx.

This is an inequality of the formx² ≤ ax + band Braess reasonably concludesx≤a+b^1/2. Since we are interested in sharp estimates, we deduce

0≤x≤a/2 +p

a²/4 +b.

This and the triangle inequality

|||u−u_h||| ≤ |||u−u_h−w_D|||+|||w_D|||

prove the first assertion. Furthermore, (3.1) and (3.9) yield Z

Ω

(Λ−Λ_h)vdx =a(u−w, v)≤ |||u−w||||||v|||for all v∈V.

Hence,

|||Λ−Λ_h|||_∗ ≤ |||u−w|||.

The triangle inequality concludes the proof for the second assertion.

4. Efficiency

This section discusses the efficiency of the global upper bound GUB and involves some notation T_DC, T_C, T_i, Ω₁, Ω₂ as follows. The set of all triangles T ∈ T along the Dirichlet boundary ΓD with contact of the discrete solution on the neighbourhood ωT := {K ∈ T

T ∩K6=∅} ⊆ {u_h=Iχ} is denoted by T_DC :=

T ∈ T

E(T)∩ E(Γ_D)6=∅ andu_h=Iχ onω_T .

The set of all trianglesT with contact of the discrete solution {u_h=Iχ} is denoted by T_C :=

T ∈ T

u_h =Iχ onT .

The set of all triangles T in some layer betweenT_DC∪ T_C and the set{Iχ≤u_h}is denoted by

T_i :={T ∈ T

∃x, y∈ N(ωT), Iχ(x) =u_h(x) &Iχ(y)< u_h(y)}.

Here,N(ω_T) :={z∈ N

z∈ω_T} denotes all nodes in the element patchω_T of the triangle T. With each T ∈ T_i we associate some (preferably interior) node zT ∈ N ∩ωT such that χ(z_T) = u_h(z_T) and set Ωb_zT := {x ∈Ω

x ∈ω_T orψ_zT(x) >0}. All elements T ∈ T_i with zT ∈ΓN form the set

T_N :={T ∈ T_i

zT ∈ΓN} and theL²-norm over all triangles in T_N reads k · k_L2(T_N).

The next theorem establishes efficiency for the global upper bound GUB := 3|||w−u_h|||+ 2|||Λ_h−JΛ_h|||_∗+ 2

Z

Ω

(χ−u_h−w_D)JΛ_hdx 1/2

+ 4|||w_D|||

from Theorem 3.2. This is the converse of

|||u−u_h|||+|||Λ−Λ_h|||_∗ ≤GUB

up to pertubation terms like oscillations with respect to node patches osc(Λ,N) and elemen- twise or edgewise oscillations defined by

osc(f,T) :=kh_T(f−fT)k_L2(Ω), osc(g,E(Γ_N)) :=kh_E(g−gE)k_L2(Γ_N)

with elementwise integral mean fT|_T := R

T fdx/|T| and edgewise integral mean gE|_E :=

R

Egds /|E|.

Theorem 4.1. For an affine obstacleχ=Iχ∈P₁(Ω) and f ∈H¹(Ω), it holds GUB.|||u−u_h|||+|||Λ−Λ_h|||_∗+ osc(Λ,N) + osc(f,T) + osc(g,E(Γ_N)) +|||w_D|||

+



 X

T∈Ti\T_C

kh²_T∇fk²_L2(ωT)





1/2

+



 X

T∈T_DC

k∇w_Dk_L2(T)kh_Tfk_L2(T)





1/2

+k∇(u−χ)k_L2(TN). Proof. The proof of Theorem4.1follows by a combination of Claim 1 - Claim 6 below.

Claim 1. It holds |||w−u_h||| ≤ |||u−u_h|||+|||Λ−Λ_h|||_∗.

Proof of Claim 1. This is already known from [Bra05].

Claim 2. It holds |||JΛ_h−Λ_h|||_∗.|||Λ_h−Λ|||_∗+ osc(Λ,N).

Proof of Claim 2. For anyv∈V, Lemma2.1shows that Z

Ω

(JΛ_h)vdx = Z

Ω

Λ_hJ vdx. This and the second assertion of Lemma 2.1lead to

Z

Ω

(Λh−JΛh)vdx = Z

Ω

(Λh−Λ)(v−J v)dx+ Z

Ω

Λ(v−J v)dx

.|||v|||(|||Λ_h−Λ|||_∗+ osc(Λ,N)). Claim 3. It holds

Z

Ω

(χ−u_h−w_D)(JΛ_h)dx . X

T∈T_DC

h_Tk∇w_Dk_L2(T)kJΛ_hk_L2(T)+ X

T∈Ti\TDC

h²_Tk∇w_Dk_L2(T)k∇(JΛ_h)k_L2(ωT)

+ X

T∈T_i\T_C

h²_Tk∇(JΛ_h)k_L2(ωT) min

qz∈(^P1(T(bΩ_zT))∩C(bΩ_zT))²

k∇(χ−u_h)−q_zk_L2(bΩ

zT). The proof of Claim 3 employs Lemma 8 of [BC04a] which is recalled here for convenient reading.

Lemma 4.2([BC04a]). Letz∈ N be either an interior point of Ω or a nonconvex boundary point (so convex corner, in particular points on straight line segments are excluded). Suppose T ∈ T, ω_T := {P

z∈N(T)ϕz >0} withz∈ω_T and set Ωbz :={x∈Ω

ψz(x)>0} ∪ω_T. Let w_h∈P₁(T)∩C(Ω) satisfyw_h(z) = 0 and 0≤w_h onΩb_z. Then, it holds

kw_hk_L2(bΩ

z) .h_z min

qz∈(^{P1(T(Ωbz))∩C(Ωbz)})²

k∇w_h−q_zk_L2(bΩ

z). Proof of Claim 3. The integralR

Ω(χ−u_h−w_D)(JΛ_h)dx is analysed for eachT ∈ T. In case that χ < u_h on ωT, (3.3) yields

JΛ_h = X

z∈N(T)

ϕ_zσb_h(ϕ_z)/

Z

Ω

ϕ_zdx = 0 on T.

For T ∈ T_i\ T_C with |∂T ∩ΓD|= 0, it holds wD = 0 on T and (u_h−χ)(zT) = 0 for some z_T ∈ N(ω_T). Furthermore, there exists somey_T ∈ N(T) with χ(y_T)< u_h(y_T). Since (3.6)

and (3.3) yield

JΛh(yT) =σbh(yT)/

Z

ϕyTdx = 0, a discrete Friedrichs’ inequality shows

kJΛ_hk_L2(T).h_Tk∇(JΛ_h)k_L2(ωT). (4.1)

This and Lemma 4.2yield Z

T

(χ−u_h)(JΛ_h)dx ≤ kχ−u_hk_L2(T)kJΛ_hk_L2(T)

.h²_Tk∇(JΛ_h)k_L2(ωT) min

qz∈(^P1(T(ΩbzT))∩C(Ωb_zT))²

k∇(χ−u_h)−q_zk_L2(bΩ

zT). IfzT ∈ΓD, Lemma4.2is not applicable. However, this case is insignificant for the following reason. Since z_T is chosen preferably as an inner node, z_T ∈ Γ_D implies N(ω_z)∩ {u_h = Iχ} ⊆ N(Γ_D). Hence σb_h(y) = 0 for all y ∈ N(T). Consequently, JΛ_h = 0 on T. In case that the isolated contact node zT belongs to a convex corner or a straight-line segment of Γ_N, it follows that

Z

T

(χ−uh)(JΛh)dx ≤h²_Tk∇(χ−uh)k_L2(T)kJΛ_hk_L2(T).

In fact, free nodes on convex corners lead to exceptional situations in some second-order positive approximation [NW02].

ForT ∈ T with|∂T ∩ΓD|>0 the integral equals Z

T

(χ−u_h−wD)(JΛ_h)dx = Z

T

(χ−u_h)(JΛ_h)dx− Z

T

wD(JΛ_h)dx

≤ kχ−u_hk_L2(T)kJΛ_hk_L2(T)+kw_Dk_L2(T)kJΛ_hk_L2(T). Since w_D = 0 on∂T \Γ_D, a Friedrichs inequality shows

Z

T

(χ−u_h−w_D)(JΛ_h)dx .kχ−u_hk_L2(T)kJΛ_hk_L2(T)+h_Tk∇w_Dk_L2(T)kJΛ_hk_L2(T). The first summand vanishes if u_h =χ on T or χ < u_h on ω_T. Otherwise, it holds T ∈ T_i and Lemma4.2 leads forz=zT ∈ N(Ω) to

kχ−uhk_L2(Ωbz).hzT min

q∈(^P1(T(bΩ_zT^))∩C(bΩ_zT⁾)²

k∇(χ−uh)−qk_L2(Ωb_zT).

The factorkJΛ_hk_L2(T)can be treated as in (4.1), except in caseu_h =χon ω_T which implies

T ∈ T_DC.

Claim 4. For anyT ∈ T, it holds (4.2) hTkJΛhk_L2(T) .hTkfk_L2(ωT)

+ min

qT∈(P₁(T(ωT))∩C(ω_T))²

k∇u_h−q_Tk_L2(ωT)+h^1/2_T k(g−q_T ·ν)k_L2(ΓN∩∂ω_T)

,

(4.3) h²_Tk∇(JΛ_h)k_L2(T).h²_Tk∇fk_L2(ωT)

+ min

qT∈(P₁(T(ωT))∩C(ω_T))²

k∇u_h−q_Tk_L2(ωT)+h^1/2_T k(g−q_T ·ν)k_L2(ΓN∩∂ω_T)

.

Proof of Claim 4. SinceJΛ_h =%_h, this is Lemma 7 in [BC04a].

Claim 5. It holds min

qT∈(P1(T(ωT))∩C(ωT))²

k∇u_h−q_Tk²_L2(ωT)+h_Tkg−q_E·νk²_L2(ΓN∩∂ω_T)

. X

E∈E(ω_T)

min

qE∈(P₁(T(ωE))∩C(ω_T))²

k∇u_h−q_Ek²_L2(ωE)+h_Tkg−q_E·νk²_L2(ΓN∩E)

. X

E∈E(ω_T)

hEk[∇u_h·ν]k²_L2(E)+ X

E∈E(ΓN)

hEkg− ∇u_h·νk²_L2(ΓN∩E),

min

qz∈(^P1(T(ΩbzT))∩C(Ωb_zT))²

k∇(χ−u_h)−qzk_L2(bΩ

zT). X

E∈E(Ω_zT)

h_Ek[∇(χ−u_h)·ν]k²_L2(E). Proof of Claim 5. The first estimate follows from (3.4) in [CB02, p. 951]. Consider an inner edgeE ∈ E(Ω) and set q_E := (∇u_h|_T1− ∇u_h|_T2)/2. This yields

k∇u_h−qEk²_L2(ωE)= 1/4k[∇u_h·ν]k²_L2(ωE)=|ω_E|/(4|E|)k[∇u_h·ν]k²_L2(E)

For any Neumann edge E ∈ E(Γ_N), ω_E consists of only one element T and we can set qE :=∇u_h|_T. This proves the second asserted estimate.

Claim 6. It holds

h_Tkf_T + Λ_hk_L2(T).k∇(w−u_h)k_L2(T)+ osc(f, T) for all T ∈ T,

h^1/2_E k[∇u_h·ν]k_L2(E).k∇(w−uh)k_L2(ωE)+ osc(f,T(ωE)) for all E∈ E(Ω),

h^1/2_E k∇u_h·ν−gk_L2(E).k∇(w−u_h)k_L2(TE)+ osc(f, T_E) + osc(g, E) for all E ∈ E(Γ_N).

Proof of Claim 6. Those estimates are well-known and follow from an error analysis for the residual estimator for Poisson problems with bubble functions [Bra07,Ver96].

This section concludes with some remarks in order to support our claim that in model examples —such as the benchmarks below— GUB is equivalent to |||u−u_h|||+|||Λ−Λ_h|||∗

up to higher-order terms.

Remark 4.3 (Comment on osc(Λ,N)). Since Λ is merely anL²(Ω) function, it is not clear a priori that the term osc(Λ,N) is of higher order. Consider the maximal open setC with u=χon C and the set U :=S

ε>0B_ε with the maximal open setB_ε withχ+ε≤uon B_ε. Then,

Λ = 0 onU and Λ =f on C.

The set C is regarded as the set of contact and the set U is the set of noncontact. Hence, the oscillations of Λ on C ∪ U are bounded by the oscillations of f, while the contributions of Λ within the free boundary F := Ω\(C ∪ U) may be not because of discontinuities.

Remark 4.4(Heuristic analysis on osc(Λ,N) = higher-order term). In simple model scenar- ios, the free boundaryF is indeed a one-dimensional submanifold, cf. Examples 1-3 below, where Λ is piecewise smooth and bounded. Therefore, we expect

Λminz∈R

kΛ−Λzk²_L2(ωz)≈ |ω_z|

close to the free boundary while it is of higher order elsewhere. For the local mesh-sizeh(s) along the parametrisation of the curve F by arc-length 0 ≤ s ≤ L := |γ| . 1, it holds

|ω_z| ≈h(s)² near any pointγ(s)∈ωz. The node patches N(F) :={z∈ N

ω_z∩ F 6=∅}

along F are coupled with J := |N(F)|+ 1 many points γ(t₀), . . . , γ(t_J), along γ with 0 = t0 < t1 < . . . < tJ = L. The nonsmooth contributions osc(Λ,N(F)) of osc(Λ,N) in the neighbourhood of γ sum up to

X

z∈N(F)

h²_z min

Λz∈R

kΛ−Λ_zk²_L2(ωz).

J

X

j=0

h(t_j)⁴.

J

X

j=0

h(t_j)³(t_j+1−t_j) (witht_J+1 :=L+h(L))

= Z L

0

h(s)³ds ≤L h³_max.h³_max.

Here hmax denotes the maximal mesh-size along γ which is relatively small compared with the maximal mesh-size max_T∈T h_T of the triangulationT for all the adaptive meshes in the numerical examples of Section 6. Therefore,

osc(Λ,N(F)).h^3/2_max is of higher order compared to |||u−u_h|||.

Remark 4.5. The remaining critical contributionk∇w_Dk_L2(T)kh_Tfk_L2(T) arises only for a relatively small number of triangles along the Dirichlet boundary. It vanishes for boundary triangles T /∈ T_DC without contact at the Dirichlet boundary. It also vanishes for piecewise affine Dirichlet data u_D.

Remark 4.6. The contribution k∇(u −χ)k_L2(T_N) vanishes for pure Dirichlet boundary problems with ΓN = ∅. This is valid for all our benchmark examples in Section 6. It also vanishes for triangles without contact at the Neumann boundary in its neighbourhoodω_T.

5. Error Estimation and Adaptive Mesh-refinement Algorithm

This section studies the five classes of a posteriori estimators from Table 5.1and explains our adaptive mesh-refinement algorithm.

Table 5.1. Classes of a posteriori error estimators used in this paper.

No Classes of error estimators Class representatives 1 explicit residual-based ηR

2 averaging ηMP1

3 equilibration ηB,ηLW,ηEQL

4 least-square ηLS

5 localisation ηCF

5.1. Five types of a posteriori error estimators.

(i) Standard residual estimator. The standard residual estimator ηR:=khT(f −Λh)k_L2(Ω)+ X

E∈E

hEk[∇u_h·νE]k²_L2(E)

!1/2

is a guaranteed upper bound of|||w−u_h|||. In all our examples,T consists of right isosceles triangles, hence |||w−u_h||| ≤ η_R [CF99]. Here, [∇u_h·ν_E] denotes the jump of [∇u_h ·ν_E] acrossE ∈ E, which is set to zero along any Dirichlet edgeE ∈ E(∂Ω).

(ii) Minimal P₁(T;R²) averaging [Car04]. The error estimator η_MP1:= min

q∈P1(T;R²)∩C(Ω;R²)

k∇u_h−qk_L2(Ω)

shows very accurate results for the Laplace equation, but only yields an upper bound for

|||w−u_h||| up to some not-displayed reliability constantC_rel.

(iii) Least-square estimator. An integration by parts yields, for any q ∈H(div,Ω) and fb=f−Λ_h with elementwise integral mean fbT ∈P₀(T), that

Z

Ω

∇(w−uh)· ∇vdx = Z

Ω

(fb−fbT)vdx+ Z

Ω

(fbT + divq)vdx + Z

Ω

(∇u_h−q)· ∇vdx. After [Rep08,RSS03,CM10], this results in the error estimator

η_LS:= min

q∈RT₀(T)C_FkfbT + divqk_L2(Ω)+k∇u_h−qk_L2(Ω)+ osc(f ,bT)/π

with (upper bounds of the) Friedrichs’ constant C_F := sup_v∈V\{0}kvk_L2(Ω)/|||v|||. Our inter- pretation of Repin’s variant (without the oscillation split) reads

η_REPIN:= min

q∈RT₀(T)C_Fkfb+ divqk_L2(Ω)+k∇u_h−qk_L2(Ω)+ osc(f ,bT)/π.

This paper studies the least-square variant η_LS rather than Repin’s majorant η_REPIN for reasons discussed in [CM10, Subsection 4.2].

(iv) Luce-Wohlmuth error estimator[LW04]. Luce and Wohlmuth suggest to solve local problems around each node on the dual triangulationT^? ofT and compute some equilibrated quantityq_LW. The dual triangulationT^? connects each triangle center mid(T),T ∈ T, with the edge midpoints mid(E(T)) and nodes N(T) and so divides each triangle T ∈ T into 6 subtriangles of area|T|/6.

Consider some node z∈ N(T) and its nodal basis function ϕ^?_z with the fine patchω_z^? :=

{ϕ^?_z > 0} of the dual triangulation T^? and its neighbouring triangles T^?(z) := {T^? ∈ T^?

z ∈ N^?(T)}. Since σ_h ∈ P₀(T) is continuous along ∂ω^?_z ∩T for any T ∈ T, q ·ν = σ_h ·ν ∈ P0(E^?(∂ω_z^?)) is well-defined on the boundary edges E^?(∂ω_z^?) of ω_z^?. With fT,z :=

−R

T(f−Λ_h)ϕ_zdx/|T^?|and the local spaces Q(T^?(z)) :=

τ_h ∈RT₀(T^?(z))

divτ_h|_T^?+f_{T ,z} = 0 onT^? ∈ T^? with

N^?(T^?)∩ N(T) ={z} andq·ν =σ_h·ν along∂ω_z^?\∂Ω , the mixed finite element method solves

q|_ω^?

z := argmin

τh∈Q(T^?(z))

kq_h−τ_hk_L2(ω_z^?).

This choice of the divergence [CMx2] differs from the original one of [LW04] for an im- proved bound for |||f−Λ_h+ divq_LW|||_? with explicitly known constants, namely

|||f−Λh+ divqLW|||_? ≤ khT(f−Λh+ divqLW)k_L2(Ω)/π.

For details cf. [CMx2]. The remaining degrees of freedom permit proper boundary fluxes and

Z

Ω

q_LW·Curlϕ^?_zdx = Z

Ω

∇u_h·Curlϕ^?_zdx for all z∈ N.

Here, Curl denotes the rotated gradient Curlv := (−∂v/∂x₂, ∂v/∂x₁). Then, the Luce- Wohlmuth error estimator reads

η_LW:=k∇u_h−q_LWk_L2(Ω)+kh_T(f −Λ_h+ divq_LW)k_L2(Ω)/π.

(v) Equilibration error estimator by Braess. Braess [BS08, Bra07] designs patchwise broken Raviart-Thomas functionsr_z ∈RT−1(T(z)) that satisfy

divrz|_T =− Z

T

(f −Λh)ϕzdx/|T| forT ∈ T(z)

[r_z·ν_E]_E =−[σ_h·ν_E]_E/2 on E∈ E(z)∩ E(∂Ω)

r_z·ν= 0 along∂ω_z\ E(∂Ω).

Eventually, the quantity qB := σ_h +P

z∈N rz ∈ RT0(T) satisfies divqB|_T = −R

T(f − Λ_h)dx/|T|. The resulting error estimator reads

η_B :=k∇u_h−q_Bk_L2(Ω)+ osc(f−Λ_h,T)/π.

(vi) Equilibration error estimator by Ladeveze. The fluxes q_L designed by Ladeveze- Leguillon [LL83] act as Neumann boundary conditions for local problems on each triangle, cf. also [AO00] for details. Given the local function spaceH_D¹(T) :=H¹(T)/Rif|T∩Γ_D|= 0 and H_D¹(T) :={v∈H¹(T)

v= 0 on∂T ∩Γ_D}otherwise, seek φ_T ∈H_D¹(T) with Z

T

φ_T · ∇vdx = Z

T

(f−Λ_h)vdx− Z

T

∇u_h· ∇vdx + Z

∂T

q_L·ν_T vds

for all v∈H_D¹(T). Then the error estimate reads

|||w−u_h||| ≤ηEQL := X

T∈T

k∇φ_Tk²_L2(T)

!1/2

.

(vii) Carstensen-Funken error estimator. The partition of unity property of the nodal basis functions leads in [CF99] to the solution of local problems on node patches: For every z∈ N seek

w_z ∈W_z:=

({v∈H_loc¹ (ω_z)

kϕ^1/2_z ∇vk_L2(ωz)<∞, v= 0 on Γ_D∩∂ω_z} ifz∈Γ_D, {v∈H_loc¹ (ωz)

kϕ^1/2_z ∇vk_L2(ωz)<∞}/R otherwise with

Z

ωz

ϕ_z∇w_z· ∇vdx = Z

ωz

ϕ_z(f −Λ_h)vdx− Z

ωz

∇u_h· ∇(ϕ_zv)dx for all v∈W_z. Then the error estimator reads

|||w−u_h||| ≤η_CF := X

z∈N

kϕ^1/2_z ∇w_zk²_L2(ωz)

!1/2

.

In the computations for ηCF and ηEQL, all the local problems are solved with fourth-order polynomials for simplicity. The computed values are regarded as very good approximations.

However, strictly speaking the values displayed for ηEQL or ηCF are lower bounds of the guaranteed upper bounds.

E1 E2

E3 E

3 E

2

E1

E2

E1 E3

Figure 5.1. Red-, blue-, and green-refinement of a triangle.

5.2. Adaptive mesh refinement algorithm and notation. This section explains our adaptive mesh refinement algorithm and notation for the global upper bound GUB(η_xyz).

Throughout this section,

|||w_D|||.kh^3/2_E ∂²_EuD/∂s²k_L2(Γ_D)

denotes the computable upper bound for|||w_D|||from Theorem 2.2. The constant hidden in .is lower than 1 as proven in [CMx1]. The quantity|||Λ_h−JΛ_h|||_∗ is estimated by its upper bound from Lemma 2.1

|||Λ_h−JΛ_h|||_∗.osc(Λ_h,N) := X

z∈N

h²_zmin

fz∈R

kΛ_h−fzk²_L2(ωz)

!1/2

.

In our computations, the constants hidden in.were set to 1 for simplicity. Moreover, there is the contact-related contribution

µ_h:=

Z

Ω

(χ−u_h−wD)JΛ_hdx 1/2

.

Automatic mesh refinement generates a sequence of meshesT₀,T₁,T₂... by succesive mesh refinement according to a bulk criterion with parameter 0<Θ≤1.

Algorithm. INPUT coarse mesh T₀, 0< θ ≤1. For level `= 0,1,2, . . . until termination do

COMPUTE discrete solution u_h on T_` (e.g. with the MATLAB routinequadprog) and Λ_h from Section 3.

ESTIMATE by

GUB(η_xyz) = (η_xyz+ osc(Λ_h,N) + 3|||w_D|||)/2 + q

µ²_h+ (η_xyz+ osc(Λ_h,N) +|||w_D|||)² withη_xyzreplaced by any of the estimatorsη_R,η_MP1,η_LS,η_LW,η_EQL, orη_CF from Section5.1.

MARK minimal set M_{`</sub> ⊆ T<sub>`} of elements such that the refinement indicators η(T)² =|T|kf −Λ_hk²_L2(T)+ X

E∈E(T)

|T|^1/2k[∇u_h]_E ·ν_Ek²_L2(E)

+ Z

T

(χ−u_h−w_D)JΛ_hdx+1 3

X

z∈N(T)

h²_zmin

fz∈R

kΛ_h−JΛ_h(z)k²_L2(ωz)

+ X

E∈E(T)∩E(Γ_D)

h³_Ek∂²_EuD/∂s²k²_L2(E) for all T ∈ T_`,

satisfy

Θ X

T∈T_`

η(T)²≤ X

T∈M_`

η(T)².

REFINE byred-refinement of elements inM_{`</sub>andred-green-blue-refinement (see Figure5.1) of further elements to avoid hanging nodes and compute T<sub>`+1}.

OUTPUT efficiency index GUB(η_xyz)/|||e||| and relative contribution η_xyz/GUB(η_xyz) for ηxyz.

6. Numerical Examples

This section studies the five benchmark examples from Table6.1.

Table 6.1. Benchmark examples and corresponding section numbers.

Sec Short name Problem data Feature

6.1 Square domain f 6=uD 6= 0, χ≡0 smooth solution 6.2 L-shaped domain f 6= 0, χ≡uD ≡0 corner singularity 6.3 Square domain f 6=χ6=uD 6= 0 cusp obstacle 6.4 Square domain f ≡1, χ= dist(x, ∂Ω), uD ≡0 1d contact zone 6.5 Square domain f =−∆χ, χ= (1−x²)(y²−1) non-affine obstacle 6.1. Benchmark Example 1. The first benchmark from [NSV03b] concerns the constant obstacle χ= Iχ ≡0 on the square domain Ω = (−1,1)² subject to smooth Dirichlet data u_D(r, ϕ) =r²−0.49 and right-hand side

f(r, ϕ) =

(−16r²+ 3.92 forr >0.7

−5.8408 + 3.92r² forr≤0.7.

The exact solution of (1.1) reads

u(r, ϕ) = max{0, r²−0.49}².

The adaptive algorithm of Section5.2ran on uniform and adaptive meshes, one is displayed in Figure 6.1 on the right-hand side. The contact zone {r < 0.7} is less refined while its boundary {r = 1} is much more refined than the remaining part of the domain caused by contributions of the extra terms µh and osc(Λh,N). Since there is no contact along the boundary ∂Ω, the critical boundary term of Remark 4.5 does not arise and there holds efficiency as discussed in part one of this paper.

Figure 6.1 displays the convergence history of the exact error for uniform and adaptive mesh refinement. Since the solution is smooth, we observe the optimal empirical conver- gence rate 1/2 for uniform mesh refinement and marginal improvement by adaptive mesh refinement. The concentration on the boundary of the contact zone indeed improves the convergence rate slightly and supports the heuristic argument of Remark 4.4. In a neigh- bourhood of the boundary between the contact zone {u = χ} and the non-contact zone {u > χ} the mesh is quasi-uniformly contributed with a relatively small local maximal mesh-size hmax.

Figure 6.2 compares the efficiency indices I_xyz := GUB(η_xzy)/|||e||| of the global upper bounds GUB(ηxzy) for the six choices of ηxyz and the relative contribution ηxyz/GUB(ηxzy).

The efficiency index is around 10 but decreases slowly to values between 1 and 2 except for