**A PRIORI AND A POSTERIORI PSEUDOSTRESS-VELOCITY**
**MIXED FINITE ELEMENT ERROR ANALYSIS FOR THE STOKES**

**PROBLEM**^{∗}

CARSTEN CARSTENSEN* ^{†}*, DONGHO KIM

*,*

^{‡}_{AND}EUN-JAE PARK

^{§}**Abstract.** The pseudostress-velocity formulation of the stationary Stokes problem allows a
Raviart–Thomas mixed ﬁnite element formulation with quasi-optimal convergence and some super-
convergent reconstruction of the velocity. This local postprocessing gives rise to some averaging a
posteriori error estimator with explicit constants for reliable error control. Standard residual-based
explicit a posteriori error estimation is shown to be reliable and eﬃcient and motivates adaptive
mesh-reﬁning algorithms. Numerical experiments conﬁrm our theoretical ﬁndings and illustrate the
accuracy of the guaranteed upper error bounds even with reduced regularity.

**Key words.** mixed ﬁnite element approximations, a posteriori error estimates, Stokes problem,
pseudostress-velocity formulation

**AMS subject classifications.**65K10, 65M12, 65M60
**DOI.**10.1137/100816237

**1. Introduction.** Adaptive mesh-reﬁning plays an important practical role in
accurate calculation of the numerical solutions of partial diﬀerential equations, espe-
cially when the continuous solutions have local singularities or sharp layers. Adaptive
mesh-reﬁning algorithms consist of successive loops of SOLVE, ESTIMATE, MARK,
and REFINE. A posteriori error estimators provide quantitative estimates for the ac-
tual error and motivate local mesh-reﬁnement. Those are computed from the known
values such as the given data of the problem and the computed numerical solutions.

Various a posteriori error estimators for ﬁnite element methods or mixed ﬁnite element methods for second order elliptic problems have already been studied in [2, 4, 12, 22, 23] and for the Stokes problem in [16, 28]. These estimators are of implicit or explicit type and based on the velocity-pressure formulation of the Stokes problem. On the other hand there is a growing interest in a posteriori error estimators which are completely free of unknown constants and lead to guaranteed upper bounds on the numerical error. See, for example, [15, 1, 29] in this direction.

The stress-velocity-pressure formulation [20] is the original physical equations for incompressible Newtonian ﬂows induced by the conservation of momentum and the constitutive law. Arnold and Falk [3] proposed the pseudostress formulation for the equations of linear elasticity which does not require symmetric stress tensors. This allows for an easy discretization via mixed ﬁnite elements developed for scalar second order elliptic equations. Cai, Lee, and Wang [9] exploited the pseudostress-velocity formulation to study least-squares methods for the Stokes system.

*∗*Received by the editors November 29, 2010; accepted for publication (in revised form) August 31,
2011; published electronically December 15, 2011. This work was supported by the WCU program
through NRF (R31-2008-000-10049-0).

http://www.siam.org/journals/sinum/49-6/81623.html

*†*Department of Computational Science and Engineering, Yonsei University, Seoul 120-749, Ko-
rea and Department of Mathematics, Humboldt-Universit¨at zu Berlin, D-10099 Berlin, Germany
(cc@math.hu-berlin.de).

*‡*University College, Yonsei University, Seoul 120-749, Korea (dhkimm@yonsei.ac.kr).

*§*Department of Computational Science and Engineering, Yonsei University, Seoul 120-749, Korea
(ejpark@yonsei.ac.kr). The research of this author was supported in part by NRF-2010-0013374.

2501

Recently, Cai and Wang [11] used the pseudostress-velocity formulation for the mixed discretization of the Stokes system. The pseudostress is nonsymmetric and the approximation of the pressure, the velocity gradient, or even the stress can be algebraically obtained from the approximate value of the pseudostress. In this paper, we establish a priori and a posteriori error estimates for Raviart–Thomas mixed ﬁnite element methods for the pseudostress-velocity formulation of the Stokes problem. We design explicit residual-based reliable and eﬃcient a posteriori error estimators with a possible application to adaptive mesh-reﬁning algorithms. Explicit constants for some averaging estimator make this asymptotically exact.

The motivation for the postprocessing to improve the velocity ﬁeld is its usage in sharp guaranteed error control; this is the subtle point in the a posteriori error analysis of nonconforming or mixed ﬁnite element technologies. The numerical experiments of this paper conﬁrm eﬃciency indices between one and three. To the best of the authors’ knowledge, this is the ﬁrst result on a posteriori error analysis of the mixed FEMs for the pseudostress-velocity formulation of the Stokes problem.

The remainder of this paper is organized as follows: Section 2 starts with the pseudostress-velocity formulation for the Stokes problem while section 3 introduces its mixed ﬁnite element approximation. Section 4 establishes some superconvergent local postprocessing of the velocity for all ﬁxed polynomial degrees. Sections 5 and 6 present an a posteriori error analysis for the stress and the velocity errors. In the ﬁnal section we present numerical experiments to validate the theoretical results of the previous sections and to explore the accuracy of the guaranteed upper error bounds.

The a posteriori stress error estimator requires a ﬁnite dimensional approximation for actual computation. Numerical realization is discussed with various adaptive strategies for mesh reﬁnement. Then, numerical examples are given with concluding remarks.

We ﬁnish this section with notation and function spaces used in this paper. Let
V^{2} be the set of two-dimensional vectors and M^{2×2} the set of 2*×*2 matrices. For
* v*= (

*v*1

*, v*2)

*,*

^{t}*= (*

**τ***τ*

*ij*)

_{2×2}

*,*and

*= (*

**σ***σ*

*ij*)

_{2×2}, deﬁne

**curlv**:=

*−*^{∂v}_{∂y}^{1} ^{∂v}_{∂x}^{1}

*−*^{∂v}_{∂y}^{2} ^{∂v}_{∂x}^{2}

*,* *∇ v* :=

*∂v*1

*∂x* *∂v*1

*∂v*2 *∂y*

*∂x* *∂v*2

*∂y*

*,*

**curlτ**:=

*∂τ*12

*∂x* *−*^{∂τ}_{∂y}^{11}

*∂τ*22

*∂x* *−*^{∂τ}_{∂y}^{21}

*,* **divτ** :=

*∂τ*11

*∂x* +^{∂τ}_{∂y}^{12}

*∂τ*21

*∂x* +^{∂τ}_{∂y}^{22}

*,*
tr* τ*:=

*τ*11+

*τ*22

*,*

*:=*

**τ v***τ*11*v*1+*τ*12*v*2

*τ*21*v*1+*τ*22*v*2

*,*
* τ* :

*:=*

**σ***i,j*

*τ**ij**σ**ij**,* * δ*:= 2

*×*2 unit matrix

*.*

Standard Sobolev spaces*H** ^{s}*(

*ω*) and

*H*

_{0}

*(*

^{s}*ω*) for

*s≥*0 with associated norm

*·*

*are employed throughout this paper and (*

_{s,ω}*·,·*)

*denotes the*

_{ω}*L*

^{2}(

*ω*) inner product. In the case

*ω*= Ω the lower index is dropped, e.g.,

*·*

*=*

_{s,Ω}*·*

*and (*

_{s}*·,·*)

_{Ω}= (

*·,·*). We deﬁne

*H*

*(*

^{−s}*ω*) := (

*H*

_{0}

*(*

^{s}*ω*))

*as the dual space of*

^{∗}*H*

_{0}

*(*

^{s}*ω*). Extending the deﬁnitions to vector- and matrix-valued functions, we let

**H***(*

^{s}*ω,*V

^{2}) (simply

**H***(*

^{s}*ω*)) and

**H***(*

^{s}*ω,*M

^{2×2}) denote the Sobolev spaces over the set of two-dimensional vector- and 2

*×*2 matrix- valued functions, respectively. Finally, we deﬁne the space

* H(div,*Ω

*,*M

^{2×2}) :=

**τ***∈ L*

^{2}(Ω

*,*M

^{2×2})

**divτ**

*∈*

**L**^{2}(Ω)

with the norm **τ**^{2}_{H}_{(div,Ω,M}2×2) := (τ*, τ*) + (divτ

*,*

**divτ**)

*.*Here, (

*·,·*) denotes the

**L**^{2}(Ω

*,*M

^{2×2}) inner product

Ω* τ* :

*and the*

**τ**dx

**L**^{2}(Ω) inner product

Ω**divτ***·***divτ***dx*
in turn. The extended*L*^{2}(Γ) product along the boundary Γ is denoted by the duality
brackets*·,·*.

For short notation on generic constants*C*, for any two real numbers or functions
or expressions*A* and*B*,

*AB* abbreviates *A≤C B.*

The point is that this multiplicative constant *C* does not depend on the local or
global mesh-sizes but may solely depend on domain Ω. Similarly,*A≈B* abbreviates
*ABA*.

**2. Pseudostress-velocity formulation.** Given a bounded simply connected
polygonal domain Ω*⊂*R^{2} with (connected) Lipschitz boundary Γ ﬁlled with a ﬂuid
of viscosity*ν >*0 and given data **f***∈ L*

^{2}(Ω) and

**g***∈*

**H**^{1}(Ω), the stationary Stokes problem for the unknown velocity

*reads*

**u, and the pressure**p(1)

*−ν u*+

*∇p*=

*−*in Ω

**f***,*div

*= 0 in Ω*

**u***,*

*=*

**u***on Γ*

**g***.*The compatibility conditions read

*∂Ω*

* g·nds*= 0 and

Ω*p dx*= 0*.*
Let ˜* σ*= (˜

*σ*

*ij*)

_{2×2}be the stress tensor and let

* (u) :=* 1
2

*∇ u*+ (

*∇*

**u)**

^{t}be the deformation rate tensor. The aforementioned Stokes problem is derived from the stress-velocity-pressure formulation, the original physical equations for incom- pressible Newtonian ﬂow, i.e.,

(2)

**div*** σ*˜ =

*in Ω*

**f***,*

˜

* σ*+

*p*2

**δ**−*ν*

**(u) =****0**in Ω

*,*div

*= 0 in Ω*

**u***,*

*=*

**u***on Γ*

**g***.*

The elimination of the stress in the above system yields the problem (1). To avoid
the diﬃculties caused by the symmetry constraint of the stress tensor we use the
(nonsymmetric) pseudostress* σ* :=

*−p*+

**δ***ν∇*of [9]. Direct algebraic calculations recover the velocity gradient, stress, and pressure; the two formulations are equivalent.

**u**For simplicity, we assume that*ν*= 1 in the Stokes problem (1).

The framework of Cai and Wang [11] enables the design of the pseudostress-
velocity formulation as follows. Let*A*:M^{2×2}*→*M^{2×2} be the deviatoric operator

*A τ* :=

*1*

**τ**−2(tr* τ*)δ for all

**τ***∈*M

^{2×2}

*.*

Note that*Ker*(*A*) =*{q δ|q*is a scalar function

*}*and

*A*is a trace-free tensor called deviatoric part. The following properties of the operator

**τ***A*are immediate:

(*A τ,σ) = (τ,Aσ),*

(*A τ,Aτ*) = (

*A*) = (τ

**τ**,**τ***,*)

**τ***−*1

2(tr* τ,*tr

*)*

**τ***,*

*A*

**τ**_{0}

*≤*

**τ**_{0}

*.*

The pseudostress

* σ*:=

*−p*+

**δ***∇*

**u**allows for the pseudostress-velocity formulation for the Stokes problem (1),

(3)

**divσ**=* f* in Ω

*,*

*A*=

**σ**− ∇**u****0**in Ω

*,*

*=*

**u***on Γ*

**g***.*The second equation of (3) is obtained from

tr (*∇ u) = divu*= 0 and tr

*=*

**σ***−*2

*p.*

The compatibility condition

Ω*pdx*= 0 implies

Ω

tr* σdx*= 0

*.*

We have the following well-known regularity results (see [11, 21]) for suﬃciently
smooth boundary Γ or for a convex domain. For**f***∈ L*

^{2}(Ω), the solutions to problems (1) and (3) satisfy

**u**∈**H**^{2}(Ω)

*∩*

**H**^{1}

_{0}(Ω),

*p∈H*

^{1}(Ω)

*/R*,

**σ**∈**H**^{1}(Ω

*,*M

^{2×2}), and (4)

**u**_{2}+

*p*

_{H}^{1}

_{(Ω)/R}+

**σ**_{1}

**f**_{0}+

**g**_{2}

*.*

Hereafter, the Stokes problem is said to be *H*^{2}-regular if estimate (4) holds;*H** ^{k+2}*-
regularity is similarly deﬁned via the shift. With

*:=*

**V**

**L**^{2}(Ω) and

**Φ**:=* H(div,*Ω

*,*M

^{2×2})

*/span{*=

**δ**}**τ***∈ H(div,*Ω

*,*M

^{2×2})

Ω

trτ*dx*= 0

*,*
the weak form for the problem (3) reads: *Find* **σ**∈**Φ** *and u∈V*

*such that*

(*A σ,τ*) + (divτ

*,*for all

**u) =**<**g**,**τ n**>

**τ***∈*

**Φ**

*,*(5)

(divσ*, v) = (f,v)* for all

*(6)*

**v**∈**V**.This problem has a unique solution from the well-known inf-sup condition in the mixed formulation and the following lemma [7, 8].

Lemma 2.1. *For all τ*

*∈*

**Φ, we have**

**τ**^{2}_{0}*A*^{1/2}**τ**^{2}_{0}+**divτ**^{2}_{−1}*.*

**3. Mixed ﬁnite element method.** Let*{T*_{h}*}* be a family of shape-regular tri-
angulations of Ω into triangles*T* of diameter*h**T*. For each*T** _{h}*, denote

*E*

*to be the set of all edges of*

_{h}*T*

*. Given*

_{h}*T*

*∈ T*

*, we let*

_{h}*E*(

*T*) be the set of its edges. Further, for an edge

*E*

*∈ E*(

*T*), we let

**t***= (*

_{E}*−n*

_{2}

*, n*1)

*be the unit tangential vector along*

^{t}*E*for the unit outward normal

**n***= (*

_{E}*n*1

*, n*2)

*to*

^{t}*E*. In what follows,

*h*

*E*stands for the length of the edge

*E∈ E*

*h*. Moreover, we deﬁne the jump [w] of

*by*

**w**[w]

*E*:= (w

*T*+)

*E**−*(w

*T**−*)

*E* if *E*=*T*+*∩T**−**,*

where**n*** _{E}*points from

*T*+into its neighbor element

*T*

*−*. For an edge

*E*=

*T*+

*∩*Γ on the boundary, the jump reﬂects boundary conditions with

*and hence [w]*

**g***E*:=* g−w.*
We deﬁne the ﬁnite element spaces associated with the decomposition of Ω

*,*

**Φ*** _{T}* :=

**τ***∈*M^{2×2}(*τ**i1**, τ**i2*)*∈RT**k*(*T*) for*i*= 1*,*2 *,*
**Φ*** _{h}*:=

**τ***∈***Φ**for all*T* *∈ T**h**,* **τ**|*T* *∈***Φ**_{T}*,*
**V*** _{h}*:=

* v*= (

*v*1

*, v*2)

^{t}*∈*

**L**^{2}(Ω)for all

*T*

*∈ T*

_{h}*, v*

*i*

*|*

_{T}*∈P*

*k*(

*T*) for

*i*= 1

*,*2

*,*where

*RT*

*k*(

*T*) is the Raviart–Thomas element of index

*k*introduced in [25], and

*P*

*k*(

*T*) is the set of polynomials of total degree

*≤k*on domain

*T*. We notice that

**Φ**

_{h}*⊂*

**Φ**and hence if

**τ**

_{h}*∈*

**Φ**

*, then*

_{h}

**τ***has continuous normal components and the constraint*

_{h}Ωtr**τ**_{h}*dx*= 0 holds.

The mixed ﬁnite element methods reads: *Find σ*

_{h}*∈*

**Φ**

_{h}*and*

**u**

_{h}*∈*

**V**

_{h}*such that*(

*A*

**σ**

_{h}*,*

**τ***) + (divτ*

_{h}

_{h}*,*

**u***) =*

_{h}*<*

**g**,**τ**

_{h}*for all*

**n**>

**τ**

_{h}*∈*

**Φ**

_{h}*,*

(7)

(divσ_{h}*, v*

*) = (f*

_{h}*,*

**v***) for all*

_{h}

**v**

_{h}*∈*

**V**

_{h}*.*(8)

By Lemma 2.1 and the discrete inf-sup condition of the*RT**k* element space (cf., [7]),
the discrete problem is well posed and has a unique solution.

We consider an interpolation operator over the space**Φ. Let ˜**Π* _{h}*denote Raviart–

Thomas interpolation operator [7] associated with the degrees of freedom onto**Φ*** _{h}*+

*span{*. We deﬁne Π

**δ**}*:*

_{h}**Φ**

*→*

**Φ**

*by*

_{h}(9) Π_{h}* τ* = ˜Π

_{h}

**τ**−

Ω(tr ˜Π_{h}* τ*)

*dx*

2*|*Ω*|* * δ* for all

**τ***∈*

**Φ**

*,*where

*|*Ω

*|*=

Ω*dx*. We notice that

Ω(trΠ_{h}* τ*)

*dx*= 0. Let

**P***be the*

_{h}

**L**^{2}projection onto

**V***with the well-known approximation property*

_{h}(10) **P**_{h}**v**−**v**_{0}*Ch*^{k+1}*| v|*

*for all*

_{k+1}

**v***∈*

**H***(Ω)*

^{k+1}*.*

Then the following two lemmas hold: *k*= 0 was treated in [11] and*k≥*1 in [10].

Lemma 3.1. *The commutative property***divΠ*** _{h}*=

**P**

_{h}**div**

*holds. Furthermore, for*

**τ***∈*

**Φ**

*∩*

**H***(Ω*

^{k+1}*,*M

^{2×2}), we have

* τ−*Π

_{h}*0*

**τ***h*

^{k+1}*|*

**τ**|*k+1*

*,*(11)

**divτ***−***div(Π**_{h}**τ)**_{0}*h*^{k+1}**divτ**_{k+1}*.*
(12)

Lemma 3.2. *For the exact solution* (σ*, u)∈* (Φ

*∩*

**H***(Ω*

^{k+1}*,*M

^{2×2}))

*×*

**H***(Ω)*

^{k+1}*of problem*(1)

*and the approximate solution*(σ

_{h}*,*

**u***)*

_{h}*of problem*(7)–(8),

*we have*

**σ**−**σ**_{h}_{0}*h*^{k+1}*| σ|*

_{k+1}*,*

**u**−**u*** _{h}*0

*h*

*(*

^{k+1}*|*

**u**|*k+1*+

*|*

**σ**|*k+1*)

*.*

*Remark* 3.3. We note that physical quantities such as pressure, gradient velocity,
and stress can be expressed as

*p*=*−*1

2tr**σ**,*∇ u*=

*A*

**σ**,*˜ =*

**σ***+ (*

**σ***∇*

**u)***=*

^{t}*+ (*

**σ***A*

**σ)***=*

^{t}*A*+

**σ**

**σ**

^{t}*.*

From these identities, the approximation of the pressure, gradient velocity, and stress can be deﬁned by

*p**h*:=*−*1

2trσ_{h}*,* (*∇ u)*

*:=*

_{h}*A*

**σ**

_{h}*,*

*˜*

**σ***:=*

_{h}*A*

**σ***+*

_{h}

**σ**

^{t}

_{h}*.*Then the following relations hold:

*p−p**h*0= 1

2trσ*−*trσ* _{h}*0

*≤*

**σ**−**σ***0*

_{h}*,*

*∇ u−*(

*∇*

**u)***0=*

_{h}*A*(σ

*−*

**σ***)0*

_{h}*≤*

**σ**−**σ***0*

_{h}*,*

* σ*˜

*−*˜

**σ***0*

_{h}*≤ A*(σ

*−*

**σ***)0+(σ*

_{h}*−*

**σ***)*

_{h}*0*

^{t}*≤*2

**σ**−**σ***0*

_{h}*.*

The estimate for **P**_{h}**u**−**u*** _{h}*0 in the following theorem allows for the error
estimates of the postprocessed velocity below.

Theorem 3.4. *Any σ∈H*

*(Ω*

^{k+1}*,*M

^{2×2})

*and*=

**f****divσ**

*∈*

**H***(Ω)*

^{k+1}*satisfy*(13)

**P**

_{h}

**u**−**u***0*

_{h}*h*

*(*

^{k+2}*|*

**σ**|*k+1*+

*|*

**divσ**

*|*

*k+1*)

*.*

*Proof.* We start with a duality argument. Let (η*, z)∈*

**Φ**

*×*be the solution to (

**V***A*) + (divτ

**η**,**τ***,*for all

**z) = 0**

**τ***∈*

**Φ**

*,*

(14)

(divη*, v) = (P*

_{h}

**u**−**u**

_{h}*,*for all

**v)***(15)*

**v**∈**V**.The convexity of Ω and a priori estimates (4) imply

**z**_{2}**P**_{h}**u**−**u**_{h}_{0} and **η**_{1}**P**_{h}**u**−**u**_{h}_{0}*.*
(16)

Since (15) and**divΠ*** _{h}*=

**P**

_{h}**div, we deduce**

**P**_{h}**u**−**u**_{h}^{2}_{0}= (P_{h}**u**−**u**_{h}*,***divη)**

= (P_{h}**u**−**u**_{h}*,***divΠ**_{h}**η)**

= (u*− u*

_{h}*,*

**divΠ**

_{h}*Subtracting (7)–(8) from (5)–(6) we obtain the error identities*

**η)**.(*A*(σ*− σ*

*)*

_{h}*,*

**τ***) + (divτ*

_{h}

_{h}*,*

**u**−**u***) = 0 for all*

_{h}

**τ**

_{h}*∈*

**Φ**

_{h}*,*(17)

(div(σ*− σ*

*)*

_{h}*,*

**v***) = 0 for all*

_{h}

**v**

_{h}*∈*

**V**

_{h}*.*(18)

The identities (17), (14), and (18) and the estimates (10)–(11) yield
**P**_{h}**u**−**u**_{h}^{2}_{0}=*−*(*A*(σ*− σ*

*)*

_{h}*,*Π

_{h}*(σ*

**η**−**η)**−*−*

**σ**

_{h}*,A*

**η)**=*−*(*A*(σ*− σ*

*)*

_{h}*,*Π

_{h}

**η**−**η) + (div(σ**−**σ***)*

_{h}*,*

**z**−**P**

_{h}

**z)***h*(

**σ**−**σ**

_{h}_{0}

**η**_{1}+

**div(σ**

*−*

**σ***)*

_{h}_{0}

**z**_{1})

*.*

(19)

To analyze the term **div(σ***− σ*

*)*

_{h}_{0}, consider

**divΠ**

*=*

_{h}

**P**

_{h}**div**and (18). For any

**w**

_{h}*∈*

**V***, this leads to*

_{h}(div(Π_{h}**σ**−**σ*** _{h}*)

*,*

**w***) = (div(σ*

_{h}*−*

**σ***)*

_{h}*,*

**w***) = 0*

_{h}*.*

Let**w*** _{h}*=

**div(Π**

_{h}

**σ**−**σ***) and deduce*

_{h}**div(Π**

_{h}

**σ**−**σ***) = 0. Hence*

_{h}**div(σ***− σ*

*)0=*

_{h}**div(σ**

*−*Π

_{h}*0=*

**σ)****divσ**

*−*

**P**

_{h}**divσ**0

*h*

^{k+1}*|*

**divσ**

*|*

*k+1*

*.*(20)

Lemma 3.2 and the inequalities (16), (19)–(20) lead to
**P**_{h}**u**−**u*** _{h}*0

*h*

*(*

^{k+2}*|*

**σ**|*k+1*+

*|*

**divσ**

*|*

*k+1*)

*.*

**4. Postprocessing.** From Remark 3.3, we note that*A σ*

*=*

_{h}

**σ***+*

_{h}*p*

*h*

*is a good approximation of*

**δ***∇*

**u. This implies an improved approximate solution of the velocity***through local postprocessing in the spirit of Stenberg [26]. Let*

**u****W**^{∗}* _{h}*=

*{*

**v**∈**L**^{2}(Ω)

**v**|

_{T}*∈P*

*k+1*(

*T*)

^{2}for all

*T*

*∈ T*

_{h}*}.*

We deﬁne**u**^{∗}_{h}*∈ W*

^{∗}*on each*

_{h}*T*with

**P***=*

_{T}

**P**

_{h}*|*

*as the solution to the system*

_{T}

**P**

_{T}

**u**

^{∗}*=*

_{h}

**u**

_{h}*,*

(21)

(*∇ u*

^{∗}

_{h}*,∇*

**v***)*

_{h}*= (σ*

_{T}

_{h}*,∇*

**v***)*

_{h}*+ (*

_{T}*p*

*h*

*,*divv

*)*

_{h}*for all*

_{T}

**v**

_{h}*∈*(

*I−*

**P***)W*

_{T}

^{∗}

_{h}*|*

_{T}*.*(22)

Note that**u**^{∗}_{h}*|** _{T}* is the Riesz representation of the linear functional (σ

_{h}*,∇*(

*·*))

*+ (*

_{T}*p*

*h*

*,*div(

*·*))

*in the Hilbert space (*

_{T}*I*

*−*

**P***)W*

_{T}

^{∗}

_{h}*|*

_{T}*≡*(

*P*

*k+1*(

*T*)

*/P*

*k*(

*T*))

^{2}with scalar product (

*∇*(

*·*)

*,∇*(

*·*))

*. The Poincar´e inequality yields positive deﬁniteness on each triangle and the well-posedness of (21)–(22) follows. Observe from (22) in case*

_{T}*k*= 0 that div

**u**

^{∗}*= 0 on each*

_{h}*T*. A generalization of this property div

_{h}

**u**

^{∗}*= 0 to*

_{h}*k*

*≥*1 appears nonobvious.

From the identity* σ*=

*−p*+

**δ***∇*and (22), we obtain the error identity (

**u***∇*(u

*−*

**u**

^{∗}*)*

_{h}*,∇*

**v)***= (σ*

_{T}*−*

**σ**

_{h}*,∇*

**v)***+ (*

_{T}*p−p*

*h*

*,*divv)

_{T}*∀*(

**v**∈*I−*

**P***)W*

_{T}

^{∗}

_{h}*|*

*T*

*.*(23)

Theorem 4.1. *Let* **u***∈* **H*** ^{k+2}*(Ω),

**σ***∈*

**H***(Ω*

^{k+1}*,*M

^{2×2}), and

*=*

**f****divσ**

*∈*

**H***(Ω)*

^{k+1}*solve*(1)

*and let*

**u**

^{∗}

_{h}*∈*

**W**

^{∗}

_{h}*satisfy*(21)–(22).

*Then we have*

**u**−**u**^{∗}* _{h}*0

*h*

*(*

^{k+2}*|*

**u**|*k+2*+

*|*

**σ**|*k+1*+

*|*divσ

*|*

*k+1*)

*.*

*Proof.* Let ˆ* u* be the

*L*

^{2}-projection of

*onto*

**u**

**W**

^{∗}*. Deﬁne*

_{h}

**v***∈*

**W**

^{∗}*by*

_{h}

**v**|*T*= (δ

*−*

**P***)( ˆ*

_{T}

**u**−**u**

^{∗}*) for each*

_{h}*T*. Then we have from (23) and Schwarz inequality

*| v|*

^{2}

_{1,T}= (

*∇*( ˆ

**u**−**u**

^{∗}*)*

_{h}*,∇*

**v)**

_{T}*−*(

*∇*

**P***( ˆ*

_{T}

**u**−**u**

^{∗}*)*

_{h}*,∇*

**v)**

_{T}= (*∇*( ˆ**u**−**u)**,∇**v)*** _{T}* + (σ

*−*

**σ**

_{h}*,∇*

**v)***+ (*

_{T}*p−p*

*h*

*,*divv)

_{T}*−*(*∇ P*

*( ˆ*

_{T}

**u**−**u**

^{∗}*)*

_{h}*,∇*

**v)**

_{T}*≤ | u*ˆ

*−*1,T

**u**|*|*1,T+

**v**|

**σ**−**σ***0,T*

_{h}*|*1,T+

**v**|*p−p*

*h*0,Tdivv0,T

+*| P*

*( ˆ*

_{T}

**u**−**u**

^{∗}*)*

_{h}*|*

_{1,T}

*|*

**v**|_{1,T}

*.*(24)

Sincedivv^{2}_{0,T} *≤*2*∇ v*

^{2}

_{0,T}, we obtain

*| v|*1,T

*≤ |*ˆ

**u***−*1,T +

**u**|

**σ**−**σ***0,T+*

_{h}*√*

2*p−p**h*0,T+*| P*

*( ˆ*

_{T}

**u**−**u**

^{∗}*)*

_{h}*|*1,T

*.*Since (δ

*−*

**P***)w= 0 if*

_{T}*0(*

**w**∈P*T*)

^{2}, the Poincar´e inequality (with Payne–Weinberger constant 1

*/π*) reads

**v**_{0,T} *≤h**T**/π| v|*

_{1,T}

*.*

This inequality and the inverse estimate

*| P*

*( ˆ*

_{T}

**u**−**u**

^{∗}*)*

_{h}*|*1,T

*h*

^{−1}

_{T}

**P***( ˆ*

_{T}

**u**−**u**

^{∗}*)0,T*

_{h}yield

(δ*− P*

*)( ˆ*

_{T}

**u**−**u**

^{∗}*)*

_{h}_{0,T}

*h*

*T*

*|*

**v**|_{1,T}

**P***( ˆ*

_{T}

**u**−**u**

^{∗}*)*

_{h}_{0,T}

+*h**T*(*| u*ˆ

*−*1,T+

**u**|

**σ**−**σ***0,T +*

_{h}*p−p*

*h*0,T)

*.*(25)

After squaring and summing over all*T* *∈ T**h*, we conclude the proof from Theorem 3.4
and the following estimates:

**P*** _{h}*( ˆ

**u**−**u**

^{∗}*)*

_{h}_{0}=

**P**

_{h}

**u**−**P**

_{h}

**u**

^{∗}

_{h}_{0}=

**P**

_{h}

**u**−**u**

_{h}_{0}

*h*

*(*

^{k+2}*|*

**σ**|*+*

_{k+1}*|*divσ

*|*

*)*

_{k+1}*,*

*p−p*

*h*0

*≤*

**σ**−**σ***0*

_{h}*h*

^{k+1}*|*

**σ**|*k+1*

*,*

**u**−**u**^{∗}_{h}_{0}*≤ u−u*ˆ

_{0}+(δ

*−*

**P***)( ˆ*

_{h}

**u**−**u**

^{∗}*)*

_{h}_{0}+

**P***( ˆ*

_{h}

**u**−**u**

^{∗}*)*

_{h}_{0}

*.*

**5. A posteriori stress error estimation.** In its ﬁrst part, this section follows
the uniﬁed approach in [13] to obtain reliable and eﬃcient error estimators. The
second part analyzes the constants explicitly and leads to asymptotic exactness in
Theorem 5.3.

Recall that * ε* :=

**σ**−**σ***and*

_{h}*:=*

**e**

**u**−**u***for the unique approximate solution (σ*

_{h}

_{h}*,*

**u***)*

_{h}*∈*

**Φ**

_{h}*×*

**V***of the mixed ﬁnite element methods (7)–(8). The well posedness of the Stokes system leads to equivalence of errors and residuals. The generic constants, hidden in the equivalence*

_{h}*≈*below, represent the norms of some operators on the continuous level (cf. section 3 of [13] for details and proofs) and are independent of

**ε**,**u**,**v**,**σ**

_{h}*,*, and

**f**

**f***:=*

_{h}

**P**

_{h}*, etc.*

**f**Theorem 5.1. *Given the exact solution σ∈*

**Φ**

*and*+

**u**∈**g**

**H**^{1}

_{0}(Ω)

*from*(5)–(6)

*and the discrete solution*

**σ**

_{h}*∈*

**Φ**

_{h}*and*

**u**

_{h}*∈*

**V**

_{h}*from*(7)–(8), any

**v***∈*

*+*

**g**

**H**^{1}

_{0}(Ω)

*satisﬁes*

**ε**_{0}+**u**−**v**_{1}*≈ A σ*

_{h}*− ∇*

**v**_{0}+

**f***−*

**f**

_{h}

_{−1}*.*

*Proof. The ideas of the proof are contained in [13] and recalled here for the*
particular situation at hand for convenient reading. The bilinear form from (5)–(6)
will be recast into the primal form with the bilinear form

*B*:

**L**^{2}(Ω;M^{2×2})*/span{ δ} ×H*

^{1}

_{0}(Ω)

*×*

**L**^{2}(Ω;M^{2×2})*/span{ δ} ×H*

^{1}

_{0}(Ω)

*→*R
deﬁned for any (σ*, u),*(τ

*,*

**v)**∈**L**^{2}(Ω;M

^{2×2})

*/span{*

**δ**} ×**H**^{1}

_{0}(Ω) by

*B*((σ*, u),*(τ

*,*)

**v)) := (**A**σ**,**τ***−*(τ

*,∇*(σ

**u)**−*,∇*

**v)**.This bilinear form is bounded and symmetric and deﬁnes an isomorphism. The Brezzi splitting theorem for the proof of the global inf-sup condition [5, 6, 7] requires the following inf-sup condition on the bilinear form:

*b*:**L**^{2}(Ω;M^{2×2})*/span{ δ} ×H*

^{1}

_{0}(Ω)

*→*R,(τ

*,*(τ

**v)**→*,∇*

**v)**.Indeed, any**v**∈**H**^{1}_{0}(Ω) with gradient* τ* :=

*∇*

**v**∈**L**^{2}(Ω;M

^{2×2})

*/span{*satisﬁes

**δ**}*b*(τ

*,*(

**v) =**b*∇*

**v**,**v) =**∇**v**^{2}

_{0}

*≈*

**v**^{2}

_{1}

**τ**_{0}

**v**_{1}

*.*

The second condition in the Brezzi splitting theorem is the ellipticity of the bilinear form

*a*:**L**^{2}(Ω;M^{2×2})*× L*

^{2}(Ω;M

^{2×2})

*→*R,(σ

*,*(

**τ)**→*A*on the subspace

**σ**,**τ)***Z*:=*{ τ∈L*

^{2}(Ω;M

^{2×2})

*/span{*for all

**δ**}|

**v**∈**H**^{1}

_{0}(Ω)

*,*(τ

*,∇*

**v) = 0**}.An integration by parts shows that **τ***∈Z* is divergence free, and hence Lemma 2.1
leads to

**τ**^{2}_{0}*A*^{1/2}**τ**^{2}_{0}=*a*(τ*, τ).*

This completes the proof of the global inf-sup condition on the bilinear form*B*. Con-
sequently, given the exact solution **σ***∈* **Φ** and **u***∈* **H**^{1}_{0}(Ω) and the discrete solu-
tion **σ**_{h}*∈* **Φ*** _{h}* and

**u**

_{h}*∈*

**V***plus any*

_{h}

**v***∈*

**H**^{1}

_{0}(Ω), the norm of (σ

*−*

**σ**

_{h}*,*in

**u**−**v)**

**L**^{2}(Ω;M

^{2×2})

*×*

**H**^{1}

_{0}(Ω) is equivalent to the dual norm of

*B*((σ

*−*

**σ**

_{h}*,*) in the dual space of

**u**−**v)**,·

**L**^{2}(Ω;M

^{2×2})

*×*

**H**^{1}

_{0}(Ω). Since

*A*=

**σ***∇*, it follows for all

**u**

**τ***∈*

**L**^{2}(Ω;M

^{2×2})

*/span{*and

**δ**}

**w**∈**H**^{1}

_{0}(Ω) that

*B*((σ*− σ*

_{h}*,*(τ

**u**−**v)**,*,*

**w)) = (**∇**v**− A**σ**

_{h}*,*) + (σ

**τ***−*

**σ**

_{h}*,∇*Since

**w)**.**divσ**=

*and*

**f****divσ**

*=*

_{h}

**f***=*

_{h}

**P**

_{h}*, an integration by parts leads to*

**f**(σ*− σ*

_{h}*,∇*(div(σ

**w) =**−*−*

**σ***)*

_{h}*,*(f

**w) =**−*−*

**f**

_{h}*,*Hence, the square of the dual norm of

**w)**.*B*((σ

*−*

**σ**

_{h}*,*) equals

**u**−**v)**,·*A σ*

_{h}*− ∇*

**v**^{2}

_{0}+

**f**−**f**

_{h}^{2}

_{−1}*.*

This is equivalent to the error norm**ε**^{2}_{0}+**u**−**v**^{2}_{1} and concludes the proof.

The optimal choice of * v* on the right-hand side leads to the deﬁnition of the
nonconformity error estimator

*μ**h*:= min

**v***∈** g*+

**H**^{1}

_{0}(Ω)

*A*

**σ**

_{h}*− ∇*0

**v***.*Then, the theorem yields

(26) * ε*0

*μ*

*h*+

**f**−**f**

_{h}*−1*

*.*

The exact computation of the optimal **v***∈ V* in

*μ*

*h*is certainly too costly, but any approximation of it will lead to a guaranteed upper error bound. Notice that a Poincar´e inequality (with the factor 1

*/π*for convex domains after Payne–Weinberger [24])

**f**−**f**_{h}* _{−1}*:= sup

**v***∈***H**^{1}_{0}(Ω)\{0}

(f*− f*

*)*

_{h}*·*

**v**dx/∇**v**_{0}

*≤*osc

*(f*

_{k}*,T*

*)*

_{h}*/π*

for the higher-order data oscillation (of order*k*+ 2 for piecewise*H** ^{k+1}*data

*) leads to*

**f**osc* _{k}*(f

*,T*

*h*)

^{2}:=

*T**∈T**h*

*h*^{2}_{T}**f**−**f**_{h}^{2}_{0,T}*.*

The choice* v*=

*and the relation*

**u***∇*=

**u***A*imply that

**σ***μ*

*h*

*≤ A*(σ

_{h}*−*0

**σ)***≤*0

**ε***.*

Therefore, in view of (26), the estimator*μ**h* is reliable and eﬃcient in the sense of
*μ**h**≤ ε*0

*μ*

*h*+

**f***−*

**f**

_{h}*−1*

*.*

In the discussion so far, the generic constants are supressed in the notation behind
the signs*≈*or. To exploit those constants in the second part of this section, consider
the spaces

*X* :=*{ v∈H*

^{1}

_{0}(Ω) : div

*= 0*

**v***}*and

*Y* :=

**τ***∈ L*

^{2}(Ω;M

^{2×2}) : for all

**v**∈X,Ω

* τ* :

*∇*= 0

**v**dx*.*
The orthogonal split

**L**^{2}(Ω;M^{2×2}) =*∇X⊕Y*

is known from another reformulation of the Stokes equations [19].

Lemma 5.2 (see [19, Lemma 2]). *The positive constant*
*c*0:= inf

*q∈L*^{2}_{0}(Ω)\{0} sup

**v***∈***H**^{1}_{0}(Ω)\{0}

Ω*q*div**v**dx

*∇ v*

_{0}

*q*

_{0}

*satisﬁes*

*for all τ*

*∈Y*

*∃ω∈L*

^{2}

_{0}(Ω)

*∇ω*= div

**τ***andc*0

*ω*0

*≤*0

**τ***.*

Numerical values for the inf-sup constant*c*0can be found in the literature [17, 27].

Recall that * ε* :=

**σ***−*

**σ***and that*

_{h}*A*=

**ε***∇*

**u**− A**σ***is its deviatoric part. The second main result of this section is the following estimate which leads to asymptotic exactness of the estimator:*

_{h}(27) *η**h*:= min

**v***∈** g*+

**H**^{1}

_{0}(Ω)(

*A*(σ

_{h}*− ∇*+

**v)***c*

^{−1}_{0}div

*)*

**v***.*

Indeed, the following estimate implies asymptotic exactness in the sense of:

*A ε*

^{2}

_{0}

*−η*

_{h}^{2}

*≤*

**f**−**f**

_{h}^{2}

_{−1}*≤*osc

*(f*

_{k}*,T*

*)*

_{h}^{2}

*/π*

^{2}

*,*which is of high order for any suﬃciently smooth right-hand side

*.*

**f**Theorem 5.3. *For any* * v∈g*+

**H**^{1}

_{0}(Ω), it holds

*η*

_{h}^{2}

*≤ A*

**ε**^{2}

_{0}

*≤η*

_{h}^{2}+

**f**−**f**

_{h}^{2}

_{−1}*.*

*Proof. GivenA*

**ε**∈L^{2}(Ω;M

^{2×2}), the linear functional

(*A ε,∇·*) :

*X*

*→*R,

*(*

**v**→*A*

**ε**,∇**v)**with the scalar product (*·,·*) in **L**^{2}(Ω;M^{2×2}) allows a unique Riesz representation
**a***∈X* in the Hilbert space*X* with respect to the scalar product (*∇·,∇·*) (which is
equivalent to the scalar product in**H**^{1}_{0}(Ω)). In other words,* a∈X* satisﬁes

(*∇ a,∇v) = (Aε,∇v)* for all

**v**∈X.The remainder*A ε− ∇a*belongs to

*Y*. Hence Lemma 5.2 leads to the representation of its distributional divergence as a distributional gradient of some

*ω∈L*

^{2}

_{0}(Ω),

div(*A ε− ∇a) =∇ω* and

*c*0

*ω*0

*≤ A*0

**ε**− ∇**a***.*The orthogonal decomposition

**L**^{2}(Ω;M

^{2×2}) =

*∇X⊕Y*, motivates the split

*A ε*

^{2}= (

*A*

**ε**,∇**a) + (**A**ε**,A**ε**− ∇**a)**.The analysis of the ﬁrst term (*A ε,∇a) involves diva*= tr

*∇*= 0 and an integration by parts,

**a**(*A ε,∇a) = (ε,∇a) = (−f*+ div

**σ**

_{h}*,*

**a)**.Since div**σ*** _{h}*=

**f***is the piecewise polynomial best approximation of*

_{h}*, (*

**f***A*(f

**ε**,∇**a) =**−*−*

**f**

_{h}*,*

**a)**≤**f**−**f**

_{h}

_{−1}*∇*

**a**_{0}

*.*

The analysis of the second term involves a test function * v∈g*+

**H**^{1}

_{0}(Ω) and utilizes

*A*=

**ε***∇*

**u**− A**σ***. Then,*

_{h}*and so*

**u**−**g**∈X*A*is orthogonal onto

**ε**− ∇**a***∇*(u

*−*i.e.,

**g),**(*A ε,Aε− ∇a*) = (

*A*(

*∇*

**v**−**σ***)*

_{h}*,A*(g

**ε**− ∇**a) + (**∇*−*

**v)**,A**ε**− ∇**a)**.The last term involves the distributional divergence of *A ε− ∇a* which equals the
gradient of the aforementioned function

*ω*. Therefore and with div

*= 0,*

**g**(*∇*(g*− v),Aε− ∇a*) = (

*ω,*div

*0div*

**v)**≤ ω*0*

**v***.*

The aforementioned estimate of*ω*0 with the inf-sup constant*c*0 allows for a con-
trol of the last term. The combination of the resulting estimate with the previous
inequalities leads to

(*A ε,Aε− ∇a)≤ Aε− ∇a*0

*A*(*∇ v−σ*

*)0+div*

_{h}*0*

**v***/c*0

*.*
The combination of the two estimates with the orthogonal sum

*A ε*

^{2}

_{0}=

*∇*

**a**^{2}

_{0}+

*A*

**ε**− ∇**a**^{2}

_{0}concludes the proof of the theorem.

Remark that *η**h* *≈* *μ**h*. This follows from the orthogonality of deviatoric and
isotropic part of matrices and Pythagoras’ theorem (*| · |* denotes the Frobenius norm
for matrices)

*|*(*A σ*

_{h}*− ∇*)

**v)(**x*|*

^{2}=

*|A*(σ

_{h}*− ∇*)

**v)(**x*|*

^{2}+

*|*tr(

*∇*

**v)**|^{2}

*/*2

*.*

The explicit residual-based error estimators can be derived following arguments in the literature [18, 13, 16, 19, 28].

Theorem 5.4. *With the tangential unit vector t*

_{E}*and the jump of the tangential*

*components of the discrete stress deviator*[t

_{E}*· A*

**σ***], it holds reliability in the sense*

_{h}**ε**_{0} *η** ^{Res}*:=

*T*

*h** _{T}*(f

*−*

**f***)*

_{h}^{2}

_{0,T}+

*T*

*h*_{T}**curl(***A σ*

*)*

_{h}^{2}

_{0,T}

+

*E*

*h*^{1/2}* _{E}* [t

_{E}*· A*

**σ***]*

_{h}^{2}

_{0,E}

_{1/2}

*.*
*Local eﬃciency holds in the sense that, for eachT* *∈ T*_{h}*,*

*h**T***curl(***A σ*

*)0,T*

_{h}*0,T;*

**ε***for each interior*

*E∈ E*

_{h}*with edge-patchω*

*E*

*,*

*h*^{1/2}* _{E}* [t

_{E}*· A*

**σ***]0,E*

_{h}*0,ω*

**ε***E*;

*for any edgeE⊂*Γ *on the boundary (ω**E* *is one element domain),*

*h*^{1/2}* _{E}* [t

_{E}*· A*

**σ***]*

_{h}_{0,E}=

*h*

^{1/2}

_{E}

**t**

_{E}*·*(

*∇*

**g**− A**σ***)*

_{h}_{0,E}

**ε**_{0,ω}

*+*

_{E}*osc(t*

_{E}*· ∇*)

**g**, E*with edge oscillations osc(γ, E*) :=

*h*

^{1/2}

_{E}*γ−γ*

_{E}_{0,E}

*forγ*:=

**t**

_{E}*·∇*

**g**and its polynomial*L*

^{2}(

*E*)

*best approximation*

*γ*

*E*

*inP*

*k*(

*E*)

^{2}

*.*

*Proof.* For the reliability in view of (26), ﬁrst note that the Stokes problem
with right-hand side (div*A*(σ*− σ*

*)*

_{h}*,*0) instead of (f

*,*(a

**g) in (1) leads to some solution***, b*)

*∈*

**H**^{1}

_{0}(Ω)

*×L*

^{2}

_{0}(Ω) with diva= 0 and

**a****H**^{1}(Ω)+*b*_{L}^{2}_{(Ω)}div*A*(σ*− σ*

*)*

_{h}*||*

**H***(Ω)*

^{−1}*A*(σ

*−*

**σ***)*

_{h}

_{L}^{2}

_{(Ω)}

*.*

The Stokes equations imply that *A*(σ *− σ*

*)*

_{h}*− ∇*+

**a***b*is divergence free in the simply connected domain Ω and hence equals some rotation

**δ****curlr**of some function

**r**∈**H**^{1}(Ω). The traces in this Helmholtz decomposition

*A*(σ*− σ*

*) =*

_{h}*∇*+

**a**−b**δ****curlr**plus diva= 0 = divu and

*A*=

**σ***∇*imply

**u****v***∈** g*min+

**H**^{1}

_{0}(Ω)

*A*

**σ**

_{h}*− ∇*0

**v***≤ A*

**σ**

_{h}*− ∇*(u

*−*0=

**a)***b*

**δ**−**curlr**0

*.*

On the other hand, the orthogonality in the Helmholtz decomposition up to boundary terms lead to

*b δ−*

**curlr**

^{2}

_{0}= (

*b*

**δ**−**curlr**

*,A*(σ

*−*

**σ***)) = (curlr*

_{h}*,A*

**σ***)*

_{h}*−*

**g**,**curlrn**

*.*Given any weak interpolation

**r***of*

_{h}*in piecewise aﬃne and globally continuous functions, the divergence of*

**r****curlr**

*exists in*

_{h}

**L**^{2}(Ω) and vanishes. Hence

**τ***:=*

_{h}**curlr**_{h}*−s δ∈*

**Φ**

*with the real number*

_{h}*s*:=

Ωtr**curlr**_{h}*dx/*2, leads to
(*A ε,τ*

*) = (e*

_{h}*,*divτ

*) = 0*

_{h}*.*

This implies*< g,*

**curlr**

_{h}*= (*

**n**>*A*

**σ**

_{h}*,*

**curlr**

*). Consequently*

_{h}*b δ−*

**curlr**

^{2}

_{0}= (curl(r

*−*

**r***)*

_{h}*,A*

**σ***)*

_{h}*−*

**g**,**curl(r**

*−*

**r***)n*

_{h}*.*

An elementwise integration by parts leads to the jumps of tangential components
**t**· A**σ*** _{h}* across interior edges of triangles and to the terms

*(*

**t**·*∇*

**g**− A**σ***) along the outer boundary. Standard stability and approximation results then conclude the proof of the asserted reliability estimate.*

_{h}The eﬃciency follows from the discrete test-function technology due to Verf¨urth
and is standard amongst the experts. The function* d*:= [t

_{E}*· A*

**σ***] =*

_{h}*−*[t

_{E}*· A*polynomial along an edge

**ε] is a***E∈ E*

*h*with a quadratic function

*ψ*

*E*deﬁned as the product of the two nodal basis functions of a conforming linear ﬁnite element scheme. With

**curl(**

*A*

**σ) =****0, one deduces**

**d**^{2}_{0,E}**ψ**^{1/2}_{E}**d**^{2}_{0,E} =*− ψ*

_{E}*[t*

**d**,

_{E}*· A*

**ε]***E*

*.*

Some extension operator*P*extof the polynomial and an integration by parts separately
on the two neighboring triangles of the edge patch *ω**E* = *{ψ*_{E}*>* 0*}* show that the
previous integral over the edge*E*equals

(ψ_{E}*P*ext**d**,**curl(***A σ*

*))*

_{h}*2(ω*

_{L}*E*)

*−*(

*A*

**ε**,**curl(ψ**

_{E}*P*ext

**d))***2(ω*

_{L}*E*)

*≤ ε*

_{L}^{2}

_{(ω}

_{E}_{)}

*|*

**ψ**

_{E}*P*ext

*1,ω*

**d**|*E*+

**ψ**

_{E}*P*ext

**d**

_{L}^{2}

_{(ω}

_{E}_{)}

**curl(**

*A*

**σ***)*

_{h}

_{L}^{2}

_{(ω}

_{E}_{)}

*.*(28)

Some stability analysis of the edge-bubble functions as well as of the extension oper- ator leads to

**d**_{0,E} **ε**_{L}^{2}_{(ω}_{E}_{)}+**curl(***A σ*

*)*

_{h}

_{L}^{2}

_{(ω}

_{E}_{)}

*.*

Together with the ﬁrst asserted eﬃciency estimate for**curl(***A σ*

*)*

_{h}

_{L}^{2}

_{(T}

_{)}, this proves the second. Since the remaining details are standard and well known for elliptic PDEs, further details are omitted.

**6. A posteriori velocity error estimation.** A duality technique allows the
control of the velocity error * e* =

**u**−**u***via the regularity results (4). Recall that*

_{h}

**u**

^{∗}*denotes the postprocessed solution of section 4 and consider*

_{h}

**e***=*

^{∗}

**u**−**u**

^{∗}*. The exact solution*

_{h}*and its piecewise polynomial*

**u***L*

^{2}best approximation

**P**

_{h}*in*

**u***P*

*k*(

*T*

*h*)

^{2}satisfy

**e**^{2}_{0}=**u**−**P**_{h}**u**^{2}_{0}+**P**_{h}**e**^{2}_{0}*.*

While the ﬁrst term on the right-hand side is expected to be of order*k*+ 1, the second
may be of higher order. This is seen from the a priori error estimates of Theorem 4.1
and**P**_{h}**u**^{∗}* _{h}*=

**P**

_{h}

**u***=*

_{h}

**u***(provided that the Stokes problem is*

_{h}

**H***-regular),*

^{k+2}**P**_{h}**e**^{∗}_{0}=**P**_{h}**e**_{0}*≤ u−u*

^{∗}

_{h}_{0}

*h*

^{k+2}*.*

The following list of a posteriori error estimates reﬂects this higher order behavior
and the leading ﬁrst-order term. The smooth right-hand side enters in terms of
oscillations of order*k*++ 1*≥k*+ 2 with= 1 for *k*= 0 while= 2 for*k≥*1,

osc_{2,k}(f*,T**h*)^{2}:=

*T**∈T**h*

**f**_{h}*∈P*min*k*(T)^{2}*h*^{2}_{T}**f***− f*

_{h}^{2}

_{0,T}

*.*

Theorem 6.1. *Provided that the Stokes problem is H*

^{2}

*-regular, it holds that*(a)

**P**

_{h}

**u**−**u**

_{h}_{0}

*osc*

_{2,k}(f

*,T*

*) + min*

_{h}**v***∈** g*+

**H**^{1}

_{0}(Ω)

*h*(

*A*

**σ**

_{h}*− ∇*

**v)**_{0}

*.*