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^{2} with boundary∂Ω and its closed subset Γ_{D} of
positive surface measure, the data of the obstacle problem are the right-hand sidesf ∈H^{1}(Ω)
and g∈L^{2}(Γ_{N}) on Γ_{N} :=∂Ω\Γ_{D} plus the Dirichlet datau_{D} ∈C^{0}(Γ_{D}) edgewise inH^{2} and
the obstacleχinH^{1}(Ω)∩C^{0}(Ω) withχ≤u_{D} along Γ_{D}. The weak form relies on the closed
and convex set

K:={v∈H^{1}(Ω)

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^{1}(Ω),

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^{1}(Ω)

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^{0}(Ω) reads

Iv= X

z∈N

v(z)ϕ_{z} ∈P_{1}(T)∩C^{0}(Ω).

The discrete version of problem (1.1) employs the discrete set
K(T) :={v_{h} ∈P_{1}(T)∩C^{0}(Ω)

v_{h} =Iu_{D} on Γ_{D} and Iχ≤v_{h} in Ω}

(with the piecewise affine functionsP_{1}(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^{1}(Ω) 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^{1}(Ω) 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^{2}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_{0}(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^{1}(Ω), setJ v∈P_{1}(T)∩C^{0}(Ω) similar to [BC04a] by
J v:= X

z∈N

vzϕz with vz :=

Z

Ω

vϕzdx/ Z

Ω

ϕzdx ∈R.

Givenf ∈L^{2}(Ω), let osc(f,N) :=

P

z∈Nh^{2}_{z}min_{f}_{z}∈Rkf−f_{z}k^{2}_{L}_{2}_{(ω}

z)

1/2

.
Lemma 2.1. Forf ∈L^{2}(Ω) and v∈H^{1}(Ω), 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^{1}(ΓD)∩C^{0}(ΓD) satisfies uD ∈ H^{2}(E)
for all E ∈ E(Γ_{D}) with the edgewise second partial derivative ∂_{E}^{2}u_{D}/∂s^{2} 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^{1}(Ω) with

wD|_{Γ}_{D} =uD|_{Γ}_{D} − Iu_{D}|_{Γ}_{D},
supp(wD)⊂[

{T ∈ T

T∩ΓD 6=∅},
kw_{D}k_{L}∞(Ω) =ku_{D} − Iu_{D}k_{L}∞(ΓD),

|||w_{D}|||.kh^{3/2}_{E} ∂_{E}^{2}uD/∂s^{2}k_{L}2(Γ_{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^{1}(Ω) 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
ζ^{−1}(z) :={y∈ N

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

For eachz∈ M set

ψ_{z} := X

y∈ζ^{−1}(z)

ϕ_{y} ∈P_{1}(T)∩C^{0}(Ω).

(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

Ωϕ_{z}dxϕ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_{1}(T)

∩C^{0}(Ω) of the extendedσbh∈ P1(T)∩C^{0}(Ω)∗

from (3.6) in the Hilbert spaceL^{2}(Ω) sat-
isfies

Z

Ω

Λ_{h}ϕ_{z}dx =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

Ω

Λ_{h}vdx 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^{2}(Ω) 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Λ_{h}dx
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^{1}(Ω) with w_{D} =u_{D}− Iu_{D} on Γ_{D} and χ−u_{h} ≤w_{D} in Ω, it holds

|||u−u_{h}||| ≤a/2 +p

a^{2}/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}^{2}u_{D}/∂s^{2}k_{L}2(Γ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Λ_{h}dx 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Λ_{h}dx
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Λ_{h}dx−σ(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Λ_{h}dx−σ(e−wD)

= Z

Ω

(χ−u_{h}−w_{D})JΛ_{h}dx−
Z

Ω

(χ−u)JΛ_{h}dx −σ(e−w_{D}).

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

Z

Ω

(χ−u)JΛ_{h}dx and 0≤σ(e−w_{D}).

Hence

Z

Ω

(u−u_{h}−w_{D})JΛ_{h}dx−σ(e−w_{D})≤
Z

Ω

(χ−u_{h}−w_{D})JΛ_{h}dx.
The combination of the previous equality and inequality with some algebra leads to

|||u−u_{h}−w_{D}|||^{2}

=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Λ_{h}dx.

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

0≤x≤a/2 +p

a^{2}/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}, Ω_{1}, Ω_{2} 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_{z}_{T} := {x ∈Ω

x ∈ω_{T} orψ_{z}_{T}(x) >0}. All elements T ∈ T_{i} with
zT ∈ΓN form the set

T_{N} :={T ∈ T_{i}

zT ∈ΓN}
and theL^{2}-norm over all triangles in T_{N} reads k · k_{L}2(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Λ_{h}dx
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_{L}2(Ω), osc(g,E(Γ_{N})) :=kh_{E}(g−gE)k_{L}2(Γ_{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_{1}(Ω) and f ∈H^{1}(Ω), it holds
GUB.|||u−u_{h}|||+|||Λ−Λ_{h}|||_{∗}+ osc(Λ,N) + osc(f,T) + osc(g,E(Γ_{N})) +|||w_{D}|||

+

X

T∈Ti\T_{C}

kh^{2}_{T}∇fk^{2}_{L}2(ωT)

1/2

+

X

T∈T_{DC}

k∇w_{D}k_{L}2(T)kh_{T}fk_{L}2(T)

1/2

+k∇(u−χ)k_{L}2(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

Ω

Λ_{h}J 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_{T}k∇w_{D}k_{L}2(T)kJΛ_{h}k_{L}2(T)+ X

T∈Ti\TDC

h^{2}_{T}k∇w_{D}k_{L}2(T)k∇(JΛ_{h})k_{L}2(ωT)

+ X

T∈T_{i}\T_{C}

h^{2}_{T}k∇(JΛ_{h})k_{L}2(ωT) min

qz∈(^{P}1(T(bΩ_{zT}))∩C(bΩ_{zT}))^{2}

k∇(χ−u_{h})−q_{z}k_{L}_{2}_{(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_{1}(T)∩C(Ω) satisfyw_{h}(z) = 0 and 0≤w_{h} onΩb_{z}. Then, it holds

kw_{h}k_{L}_{2}_{(b}_{Ω}

z) .h_{z} min

qz∈(^{P}^{1}^{(T}^{(}^{Ω}^{b}^{z}^{))∩C(}^{Ω}^{b}^{z}^{)})^{2}

k∇w_{h}−q_{z}k_{L}_{2}_{(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

Ω

ϕ_{z}dx = 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Λ_{h}k_{L}2(T).h_{T}k∇(JΛ_{h})k_{L}2(ωT).
(4.1)

This and Lemma 4.2yield Z

T

(χ−u_{h})(JΛ_{h})dx ≤ kχ−u_{h}k_{L}2(T)kJΛ_{h}k_{L}2(T)

.h^{2}_{T}k∇(JΛ_{h})k_{L}2(ωT) min

qz∈(^{P}^{1}^{(T}^{(}^{Ω}^{b}zT))∩C(Ωb_{zT}))^{2}

k∇(χ−u_{h})−q_{z}k_{L}_{2}_{(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^{2}_{T}k∇(χ−uh)k_{L}2(T)kJΛ_{h}k_{L}2(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_{h}k_{L}2(T)kJΛ_{h}k_{L}2(T)+kw_{D}k_{L}2(T)kJΛ_{h}k_{L}2(T).
Since w_{D} = 0 on∂T \Γ_{D}, a Friedrichs inequality shows

Z

T

(χ−u_{h}−w_{D})(JΛ_{h})dx .kχ−u_{h}k_{L}2(T)kJΛ_{h}k_{L}2(T)+h_{T}k∇w_{D}k_{L}2(T)kJΛ_{h}k_{L}2(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_{L}2(Ωbz).hzT min

q∈(^{P}^{1}^{(T}^{(b}^{Ω}_{zT}^{))∩C(b}^{Ω}_{zT}^{)})^{2}

k∇(χ−uh)−qk_{L}2(Ωb_{zT}).

The factorkJΛ_{h}k_{L}2(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_{L}2(T) .hTkfk_{L}2(ωT)

+ min

qT∈(P_{1}(T(ωT))∩C(ω_{T}))^{2}

k∇u_{h}−q_{T}k_{L}2(ωT)+h^{1/2}_{T} k(g−q_{T} ·ν)k_{L}2(ΓN∩∂ω_{T})

,

(4.3) h^{2}_{T}k∇(JΛ_{h})k_{L}2(T).h^{2}_{T}k∇fk_{L}2(ωT)

+ min

qT∈(P_{1}(T(ωT))∩C(ω_{T}))^{2}

k∇u_{h}−q_{T}k_{L}2(ωT)+h^{1/2}_{T} k(g−q_{T} ·ν)k_{L}2(Γ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))^{2}

k∇u_{h}−q_{T}k^{2}_{L}2(ωT)+h_{T}kg−q_{E}·νk^{2}_{L}2(ΓN∩∂ω_{T})

. X

E∈E(ω_{T})

min

qE∈(P_{1}(T(ωE))∩C(ω_{T}))^{2}

k∇u_{h}−q_{E}k^{2}_{L}2(ωE)+h_{T}kg−q_{E}·νk^{2}_{L}2(ΓN∩E)

. X

E∈E(ω_{T})

hEk[∇u_{h}·ν]k^{2}_{L}2(E)+ X

E∈E(ΓN)

hEkg− ∇u_{h}·νk^{2}_{L}2(ΓN∩E),

min

qz∈(^{P}^{1}^{(T}^{(}^{Ω}^{b}zT))∩C(Ωb_{zT}))^{2}

k∇(χ−u_{h})−qzk_{L}_{2}_{(b}_{Ω}

zT). X

E∈E(Ω_{zT})

h_{E}k[∇(χ−u_{h})·ν]k^{2}_{L}2(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}|_{T}_{1}− ∇u_{h}|_{T}_{2})/2. This yields

k∇u_{h}−qEk^{2}_{L}2(ωE)= 1/4k[∇u_{h}·ν]k^{2}_{L}2(ωE)=|ω_{E}|/(4|E|)k[∇u_{h}·ν]k^{2}_{L}2(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_{T}kf_{T} + Λ_{h}k_{L}2(T).k∇(w−u_{h})k_{L}2(T)+ osc(f, T) for all T ∈ T,

h^{1/2}_{E} k[∇u_{h}·ν]k_{L}2(E).k∇(w−uh)k_{L}2(ωE)+ osc(f,T(ωE)) for all E∈ E(Ω),

h^{1/2}_{E} k∇u_{h}·ν−gk_{L}2(E).k∇(w−u_{h})k_{L}2(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^{2}(Ω) 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^{2}_{L}2(ω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)^{2} 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_{0}), . . . , γ(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^{2}_{z} min

Λz∈R

kΛ−Λ_{z}k^{2}_{L}2(ωz).

J

X

j=0

h(t_{j})^{4}.

J

X

j=0

h(t_{j})^{3}(t_{j+1}−t_{j}) (witht_{J+1} :=L+h(L))

= Z L

0

h(s)^{3}ds ≤L h^{3}_{max}.h^{3}_{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_{D}k_{L}2(T)kh_{T}fk_{L}2(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_{L}2(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_{L}2(Ω)+ X

E∈E

hEk[∇u_{h}·νE]k^{2}_{L}2(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_{1}(T;R^{2}) averaging [Car04]. The error estimator
η_{MP1}:= min

q∈P1(T;R^{2})∩C(Ω;R^{2})

k∇u_{h}−qk_{L}2(Ω)

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_{0}(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_{0}(T)C_{F}kfbT + divqk_{L}2(Ω)+k∇u_{h}−qk_{L}2(Ω)+ osc(f ,bT)/π

with (upper bounds of the) Friedrichs’ constant C_{F} := sup_{v∈V}_{\{0}}kvk_{L}2(Ω)/|||v|||. Our inter-
pretation of Repin’s variant (without the oscillation split) reads

η_{REPIN}:= min

q∈RT_{0}(T)C_{F}kfb+ divqk_{L}2(Ω)+k∇u_{h}−qk_{L}2(Ω)+ 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_{0}(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})ϕ_{z}dx/|T^{?}|and the local spaces
Q(T^{?}(z)) :=

τ_{h} ∈RT_{0}(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}−τ_{h}k_{L}2(ω_{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_{L}2(Ω)/π.

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

Z

Ω

q_{LW}·Curlϕ^{?}_{z}dx =
Z

Ω

∇u_{h}·Curlϕ^{?}_{z}dx for all z∈ N.

Here, Curl denotes the rotated gradient Curlv := (−∂v/∂x_{2}, ∂v/∂x_{1}). Then, the Luce-
Wohlmuth error estimator reads

η_{LW}:=k∇u_{h}−q_{LW}k_{L}2(Ω)+kh_{T}(f −Λ_{h}+ divq_{LW})k_{L}2(Ω)/π.

(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_{B}k_{L}2(Ω)+ 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}^{1}(T) :=H^{1}(T)/Rif|T∩Γ_{D}|= 0
and H_{D}^{1}(T) :={v∈H^{1}(T)

v= 0 on∂T ∩Γ_{D}}otherwise, seek φ_{T} ∈H_{D}^{1}(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}^{1}(T). Then the error estimate reads

|||w−u_{h}||| ≤ηEQL := X

T∈T

k∇φ_{T}k^{2}_{L}2(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}^{1} (ω_{z})

kϕ^{1/2}_{z} ∇vk_{L}2(ωz)<∞, v= 0 on Γ_{D}∩∂ω_{z}} ifz∈Γ_{D},
{v∈H_{loc}^{1} (ωz)

kϕ^{1/2}_{z} ∇vk_{L}2(ωz)<∞}/R otherwise
with

Z

ωz

ϕ_{z}∇w_{z}· ∇vdx =
Z

ωz

ϕ_{z}(f −Λ_{h})vdx−
Z

ωz

∇u_{h}· ∇(ϕ_{z}v)dx for all v∈W_{z}.
Then the error estimator reads

|||w−u_{h}||| ≤η_{CF} := X

z∈N

kϕ^{1/2}_{z} ∇w_{z}k^{2}_{L}2(ω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} ∂^{2}_{E}uD/∂s^{2}k_{L}2(Γ_{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^{2}_{z}min

fz∈R

kΛ_{h}−fzk^{2}_{L}2(ω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Λ_{h}dx
1/2

.

Automatic mesh refinement generates a sequence of meshesT_{0},T_{1},T_{2}... by succesive mesh
refinement according to a bulk criterion with parameter 0<Θ≤1.

Algorithm. INPUT coarse mesh T_{0}, 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

µ^{2}_{h}+ (η_{xyz}+ osc(Λ_{h},N) +|||w_{D}|||)^{2}
withη_{xyz}replaced by any of the estimatorsη_{R},η_{MP1},η_{LS},η_{LW},η_{EQL}, orη_{CF} from Section5.1.

MARK minimal set M_{`} ⊆ T_{`} of elements such that the refinement indicators
η(T)^{2} =|T|kf −Λ_{h}k^{2}_{L}2(T)+ X

E∈E(T)

|T|^{1/2}k[∇u_{h}]_{E} ·ν_{E}k^{2}_{L}2(E)

+ Z

T

(χ−u_{h}−w_{D})JΛ_{h}dx+1
3

X

z∈N(T)

h^{2}_{z}min

fz∈R

kΛ_{h}−JΛ_{h}(z)k^{2}_{L}2(ωz)

+ X

E∈E(T)∩E(Γ_{D})

h^{3}_{E}k∂^{2}_{E}uD/∂s^{2}k^{2}_{L}2(E) for all T ∈ T_{`},

satisfy

Θ X

T∈T_{`}

η(T)^{2}≤ X

T∈M_{`}

η(T)^{2}.

REFINE byred-refinement of elements inM_{`}andred-green-blue-refinement (see Figure5.1)
of further elements to avoid hanging nodes and compute T_{`+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^{2})(y^{2}−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)^{2} subject to smooth Dirichlet data
u_{D}(r, ϕ) =r^{2}−0.49 and right-hand side

f(r, ϕ) =

(−16r^{2}+ 3.92 forr >0.7

−5.8408 + 3.92r^{2} forr≤0.7.

The exact solution of (1.1) reads

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

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