∎

¹¹institutetext: X. Chen ²²institutetext: Department of Mechanical Engineering, Pennsylvania State University, University Park, PA 16802, USA
²²email: [email protected] ³³institutetext: Y. Li ⁴⁴institutetext: Department of Mathematics, The Pennsylvania State University, University Park, PA 16802, USA
⁴⁴email: [email protected]

Superconvergent pseudostress-velocity finite element methods for the Oseen and Navier–Stokes equations

Xi Chen Yuwen Li

(Received: / Accepted: date)

Abstract

We present a priori and superconvergence error estimates of mixed finite element methods for the pseudostress-velocity formulation of the Oseen equation. In particular, we derive superconvergence estimates for the velocity and a priori error estimates under unstructured grids, and obtain superconvergence results for the pseudostress under certain structured grids. A variety of numerical experiments validate the theoretical results and illustrate the effectiveness of the superconvergent recovery-based adaptive mesh refinement. It is also numerically shown that the proposed postprocessing yields apparent superconvergence in a benchmark problem for the incompressible Navier–Stokes equation.

Keywords:

pseudostresssuperconvergencesuperclosenesspostprocessingOseen equationNavier–Stokes equation

MSC:

65N1265N1565N30

1 Introduction

In fluid mechanics, the Oseen equations describe the flow of a viscous and incompressible fluid at small Reynolds numbers. Instead of dropping the advective term completely, the Oseen approximation linearizes the advective acceleration term $\bm{u}\cdot\nabla\bm{u}$ to $\bm{b}\cdot\nabla\bm{u}$ , where $\bm{b}$ represents the velocity at large distance. As a result it provides a lowest-order solution that is uniformly valid everywhere in the flow field. Since the Oseen equations partly account for the inertia terms (at large distance), they have better approximation in the far field while keeping the same order of accuracy as Stokes approximation near the body, see, e.g., kundu2015fluid . Mathematically speaking, the Oseen equation can be viewed as a linearized Navier–Stokes equation arising from fixed point iteration.

Let $\Omega\subset\mathbb{R}^{d}$ be a bounded Lipschiz domain with $d\in\{2,3\}$ , $\bm{b},\bm{f}:\Omega\rightarrow\mathbb{R}^{d}$ , $c:\Omega\rightarrow\mathbb{R}$ and $\bm{g}:\partial\Omega\rightarrow\mathbb{R}^{d}$ be given data. Let $\bm{u}:\Omega\rightarrow\mathbb{R}^{d}$ denote the velocity field and $p:\Omega\rightarrow\mathbb{R}$ denote the pressure subject to the zero-mean constraint

\int_{\Omega}pdx=0.

(1.1)

The Oseen equation under the Dirichlet boundary condition is


$\displaystyle-\nu\Delta\bm{u}+\bm{b}\cdot\nabla\bm{u}+c\bm{u}+\nabla p$	$\displaystyle=\bm{f}\text{ in }\Omega,$	(1.2a)
$\displaystyle\nabla\cdot\bm{u}$	$\displaystyle=0\text{ in }\Omega,$	(1.2b)
$\displaystyle\bm{u}$	$\displaystyle=\bm{g}\text{ on }\partial\Omega,$	(1.2c)

where $\nu>0$ is the viscosity constant. As a default, vectors such as $\bm{b},\bm{f},\bm{u},\bm{g}$ are always arranged as columns unless confusion arises. In (1.2a), $\nabla\bm{u}$ is the Jacobian matrix of $\bm{u},$ and we adopt the convention $\bm{b}\cdot\nabla\bm{u}:=(\nabla\bm{u})\bm{b}$ . When $\bm{b}=\bm{0}$ (resp. $\bm{b}=\bm{0},c=0$ ), (1.2) reduces to the Brinkman equation (resp. Stokes equation). Replacing $\bm{b}$ with $\bm{u}$ in (1.2) yields the Navier–Stokes equation


$\displaystyle-\nu\Delta\bm{u}+\bm{u}\cdot\nabla\bm{u}+\nabla p$	$\displaystyle=\bm{f}\text{ in }\Omega,$	(1.3a)
$\displaystyle\nabla\cdot\bm{u}$	$\displaystyle=0\text{ in }\Omega,$	(1.3b)
$\displaystyle\bm{u}$	$\displaystyle=\bm{g}\text{ on }\partial\Omega.$	(1.3c)

Numerical analysis of Oseen, Brinkman, Stokes, or Navier–Stokes equations based on the velocity-pressure formulation is extensive, see, e.g., GiraultRaviart1986 ; BoffiBrezziFortin2013 ; GuzmanNeilan2014 ; FalkNeilan2013 for conforming mixed methods, ccs2006 ; CCNP2013 ; CockburnSayas2014 ; CesmeliogluCockburnQiu2017 ; FuQiu2019 for hybridized discontinuous Galerkin methods, Mulin2014 ; WangYe2016 ; LLC2016 for weak Galerkin methods, and ABMV2014 for virtual element methods.

Let $\bm{I}\in\mathbb{R}^{d\times d}$ denote the identity matrix. The nonsymmetric pseudostress and symmetric stress are

	$\displaystyle\bm{\sigma}$	$\displaystyle:=\nu\nabla\bm{u}-p\bm{I},$
	$\displaystyle\bm{\sigma}^{S}$	$\displaystyle:=\frac{\nu}{2}(\nabla\bm{u}+(\nabla\bm{u})^{\top})-p\bm{I},$

respectively. In a series of papers Caiwang2007 ; CaiWangZhang2010a ; CaiTongVassilevskiWang2010b , Cai et al. proposed the pseudostress-velocity formulation of the Stokes and Navier–Stokes equations, which could be discretized by classical Raviart–Thomas or Brezzi–Douglas–Marini elements without any stabilization term, see (1.12). Let $\text{Tr}\bm{\sigma}$ be the trace of the matrix $\bm{\sigma}$ and

\mathcal{A}\bm{\sigma}:=\bm{\sigma}-\frac{1}{d}(\text{Tr}\bm{\sigma})\bm{I}

denote the deviatoric part of $\bm{\sigma}$ . Direct calculation shows that the Oseen equation (1.2) equation is equivalent to the following pseudostress-velocity formulation

	$\displaystyle\mathcal{A}\bm{\sigma}$	$\displaystyle=\nu\nabla\bm{u},$		(1.4)
	$\displaystyle-\nabla\cdot\bm{\sigma}+\nu^{-1}\bm{b}\cdot\mathcal{A}\bm{\sigma}+c\bm{u}$	$\displaystyle=\bm{f}$		(1.4)

subject to the constraints

\bm{u}|_{\partial\Omega}=\bm{g},\quad\int_{\Omega}\text{Tr}\bm{\sigma}dx=0.

Here $\nabla\cdot\bm{\sigma}$ is the row-wise divergence of the matrix-valued function $\bm{\sigma}$ .

In contrast to extensive numerical results on velocity-pressure formulation of Stokes-related equations, numerical analysis of pseudostress-based methods is restricted to classical mixed methods CaiWangZhang2010a ; Gatica2014 , virtual element methods Caceres2017 ; Caceres2017b , discontinuous Galerkin methods Qian2019 , and adaptive mixed methods CKP2011 ; CGS2013 ; HuYu2018 ; Li2020MCOM for the Stokes or Brinkman equation. There are a few works devoted to the mixed formulation of the Oseen equation (1.4), see, e.g., Park2014 for a lowest order upwinded mixed method on rectangular meshes and BarriosCasconGonzalez2017 ; BarriosCasconGonzalez2020 for least-squares mixed methods. It is the purpose of this paper to shed light on a priori and superconvergence analysis of pseudostress-velocity mixed methods for the Oseen and Navier–Stokes equations.

In this paper, we shall develop a priori and superconvergence error estimates for (1.4). To the best of our knowledge, even a priori error estimates of (1.4) is not available in literature. We emphasize that the presence of lower order terms is a major challenge in error analysis of mixed methods, see, e.g., DouglasRoberts1985 ; Demlow2002 for elliptic equations with lower order terms and ArnoldLi2017 for the perturbed Hodge–Laplace equation. For instance, due to the convection term $\bm{b}\cdot\mathcal{A}\bm{\sigma}$ and variable coefficient $c$ , the construction of the discrete inf-sup condition is far from obvious. Hence we shall use the duality argument to derive a priori error estimates, which also yields improved decoupled error estimates. It is noted that mixed methods for the scalar elliptic equation

-\nabla\cdot(\bm{A}\nabla u+\bm{b}u)+cu=f

(1.5)

have been analyzed in DouglasRoberts1985 ; Demlow2002 , where $\bm{A}$ is a uniformly elliptic matrix-valued coefficient. In contrast, the deviatoric operator $\mathcal{A}$ in (1.2) is singular, which is a major difficulty in our analysis.

Furthermore, we shall derive superconvergence results for $\bm{\sigma}$ based on postprocessing and a new estimate on $\|\Pi_{h}\bm{\sigma}-\bm{\sigma}_{h}\|$ for several families of RT elements. Even for the Stokes equation, such estimates are not known in literature. The superconvergence postprocessing procedure for $\bm{\sigma}_{h}$ and $\bm{u}_{h}$ is element-wise and flexible, see Section 4 and see, e.g., ZZ1992a ; CKP2011 ; BankXu2003a ; ZhangNaga2005 ; HuangLiWu2013 ; BankLi2019 ; Li2021JSC for related postprocessing schemes in the literature. Since the postprocessing scheme is independent of the underlying physical model, we also test the proposed postprocessing in the benchmark problem for steady incompressible Navier–Stokes equation, where apparent superconvergence phenomena in $\bm{\sigma}$ and $\bm{u}$ are observed. As for superconvergence of Stokes-type equations in velocity-pressure form, readers are referred to e.g., WangYe2001 ; Ye2001 ; CuiYe2009 for superconvergence by two-grid $L^{2}$ -projections, Pan1997 ; Eichel2011 ; LN2008 for superconvergent recovery of lowest-order methods, and CockburnSayas2014 ; FuQiu2019 for superconvergent hybridized discontinuous Galerkin methods.

1.1 Preliminary notation

Given a vector space $X$ , let $[X]^{n}$ denote the Cartesian product of $n$ copies of $X$ and $[X]^{n\times n}$ the space of $n\times n$ matrices whose components are contained in $X$ . Let $H(\operatorname{div},\Omega)=\{\bm{v}\in[L^{2}(\Omega)]^{d}:\nabla\cdot\bm{v}\in L^{2}(\Omega)\}$ be the usual $H(\operatorname{div})$ space. Let

	$\displaystyle\bm{V}$	$\displaystyle=[L^{2}(\Omega)]^{d},$
	$\displaystyle\bm{\Sigma}$	$\displaystyle=\{\bm{\tau}\in[L^{2}(\Omega)]^{d\times d}:\bm{\tau}_{i}\in H(\operatorname{div},\Omega),~{}1\leq i\leq d,~{}\int_{\Omega}\text{Tr}\bm{\tau}dx=0\},$

where $\bm{\tau}_{i}$ is the $i$ -th row of $\bm{\tau}.$ For scalar-, vector-, or matrix-valued functions, we use $(\cdot,\cdot)$ and $\langle\cdot,\cdot\rangle$ to denote the $L^{2}(\Omega)$ - and $L^{2}(\partial\Omega)$ -inner product, respectively. The variational formulation of (1.4) seeks $\bm{\sigma}\in\bm{\Sigma}$ and $\bm{u}\in\bm{V}$ such that

	$\displaystyle(\mathcal{A}\bm{\sigma},\bm{\tau})+\nu(\nabla\cdot\bm{\tau},\bm{u})$	$\displaystyle=\nu\langle\bm{g},\bm{\tau}\bm{n}\rangle,\quad\bm{\tau}\in\bm{\Sigma},$		(1.6)
	$\displaystyle-\nu(\nabla\cdot\bm{\sigma},\bm{v})+(\bm{b}\cdot\mathcal{A}\bm{\sigma},\bm{v})+\nu(c\bm{u},\bm{v})$	$\displaystyle=\nu(\bm{f},\bm{v}),\quad\bm{v}\in\bm{V},$		(1.6)

where $\bm{n}$ is the outward normal to $\partial\Omega.$

Let $\mathcal{T}_{h}$ be a conforming simplicial or cubical partition of $\Omega$ . Given an element $K\in\mathcal{T}_{h},$ $r_{K}$ and $\rho_{K}$ denote radii of circumscribed and inscribed spheres of $K$ , $h_{K}$ is the diameter of $K$ , and $h:=\max_{K\in\mathcal{T}_{h}}h_{K}<1$ is the mesh size of $\mathcal{T}_{h}.$ We assume that the aspect ratio of elements in $\mathcal{T}_{h}$ is uniformly bounded, i.e.,

\max_{K\in\mathcal{T}_{h}}r_{K}/\rho_{K}:={C}_{\text{mesh}}<\infty

for some absolute constant $C_{\text{mesh}}$ . Given an integer $r\geq 0,$ let $\widetilde{\bm{\Sigma}}_{r}(K)$ , $V_{r}(K)$ be suitable finite-dimensional vector spaces defined on $K$ and

	$\displaystyle\widetilde{\bm{\Sigma}}_{h}$	$\displaystyle:=\{\bm{\tau}\in H(\operatorname{div},\Omega):\bm{\tau}\|_{K}\in\widetilde{\bm{\Sigma}}_{r}(K)~{}\text{for}~{}K\in\mathcal{T}_{h}\},$
	$\displaystyle V_{h}$	$\displaystyle:=\{v\in L^{2}(\Omega):v\|_{K}\in V_{r}(K)~{}\text{for}~{}K\in\mathcal{T}_{h}\}.$

For instance, when $\mathcal{T}_{h}$ is a simplicial partition, one could take

\widetilde{\bm{\Sigma}}_{r}(K):=[\mathcal{P}_{r}(K)]^{d}+\mathcal{P}_{r}(K)\bm{x},\quad V_{r}(K):=\mathcal{P}_{r}(K),

(1.7)

where $\mathcal{P}_{r}(K)$ is the space of polynomials of degree at most $r$ on $K$ , and $\bm{x}=(x_{1},\ldots,x_{d})^{\top}$ denotes the position vector. In this case, the corresponding $\widetilde{\bm{\Sigma}}_{h}\times V_{h}$ is the classical Raviart–Thomas (RT) RaviartThomas1977 space. Another possible choice is

\widetilde{\bm{\Sigma}}_{r}(K):=[\mathcal{P}_{r+1}(K)]^{d},\quad V_{r}(K):=\mathcal{P}_{r}(K),

(1.8)

which corresponds to the Brezzi–Douglas–Marini (BDM) BrezziDouglasMarini1985 element.

Given a hypercube $K\subset\mathbb{R}^{d}$ , let $\mathcal{Q}_{r}(K)$ denote the space of polynomials on $K$ of degree at most $r$ in $x_{i}$ with $1\leq i\leq d$ . For a cubical mesh $\mathcal{T}_{h}$ , the pair $\widetilde{\bm{\Sigma}}_{h}\times V_{h}$ can be the cubical RT element using the shape functions

\widetilde{\bm{\Sigma}}_{r}(K):=[\mathcal{Q}_{r}(K)]^{d}+\bm{x}\mathcal{Q}_{r}(K),\quad V_{r}(K):=\mathcal{Q}_{r}(K).

(1.9)

Let $\widetilde{\Pi}_{h}:[H^{1}(\Omega)]^{d}\rightarrow\widetilde{\bm{\Sigma}}_{h}$ be the canonical interpolation onto $\widetilde{\bm{\Sigma}}_{h}$ determined by the degrees of freedom of $\widetilde{\bm{\Sigma}}_{h}$ . Let $P_{h}:L^{2}(\Omega)\rightarrow V_{h}$ be the $L^{2}$ -projection onto $V_{h},$ i.e.,

(P_{h}v,w)=(v,w),\quad v\in L^{2}(\Omega),~{}w\in V_{h}.

(1.10)

Throughout the rest of this paper, we assume the commutativity property

\nabla\cdot\widetilde{\Pi}_{h}=P_{h}\nabla\cdot,

(1.11)

which is satisfied by at least the aforementioned RT, BDM, and cubical RT elements. Readers are referred to e.g., BrezziDouglasMarini1985 ; BDDF1987 ; BDFM1987 for other cubical or rectangular mixed elements that satisfy (1.11). Let $\bm{\Sigma}_{h}$ be the matrix version of $\widetilde{\bm{\Sigma}}_{h}$ , i.e.,

\displaystyle\bm{\Sigma}_{h}:=\{\bm{\tau}\in\bm{\Sigma}:\bm{\tau}_{i}\in\widetilde{\bm{\Sigma}}_{h}\text{ for }1\leq i\leq d\}.

Let $\bm{V}_{h}$ be the vector-valued broken piecewise polynomial space

\displaystyle\bm{V}_{h}:=\{\bm{v}=(v_{1},v_{2},\ldots,v_{d})^{\top}\in\bm{V}:{v}_{i}\in{V}_{h}\text{ for }1\leq i\leq d\}.

The mixed method for (1.4) seeks $(\bm{\sigma}_{h},\bm{u}_{h})\in\bm{\Sigma}_{h}\times\bm{V}_{h}$ satisfying

	$\displaystyle(\mathcal{A}{\bm{\sigma}_{h}},\bm{\tau})+\nu(\nabla\cdot\bm{\tau},\bm{u}_{h})$	$\displaystyle=\nu\langle\bm{g},\bm{\tau}\bm{n}\rangle,\quad\bm{\tau}\in\bm{\Sigma}_{h},$		(1.12)
	$\displaystyle-\nu(\nabla\cdot{\bm{\sigma}_{h}},\bm{v})+(\bm{b}\cdot\mathcal{A}{\bm{\sigma}_{h}},\bm{v})+\nu(c\bm{u}_{h},\bm{v})$	$\displaystyle=\nu(\bm{f},\bm{v}),\quad\bm{v}\in\bm{V}_{h}.$		(1.12)

The outline of this work is as follows. In Section 2, we develop a apriori error estimates of the mixed method (1.12) and supercloseness estimate on $\|P_{h}\bm{u}-\bm{u}_{h}\|$ . In Section 3, we develop supercloseness estimate on $\|\Pi_{h}\bm{\sigma}-\bm{\sigma}_{h}\|$ . Section 4 is devoted to superconvergent postprocessing procedures. Numerical results for smooth, singular, convection-dominated, and nonlinear problems are presented in Section 5. Throughout the rest of this paper, we say $A\lesssim B$ provided $A\leq CB,$ where $C>0$ is a generic constant dependent solely on $\bm{b}$ , $c$ , $\Omega$ , $C_{\text{mesh}}$ , $\nu$ . In the error analysis, we assume $\nu=1$ without loss of generality.

2 A priori error estimates

Let $\|\cdot\|_{r}$ denote the $\|\cdot\|_{H^{r}(\Omega)}$ -norm, $|\cdot|_{r}$ be the $|\cdot|_{H^{r}(\Omega)}$ -semi norm, $\|\cdot\|:=\|\cdot\|_{0}$ be the $\|\cdot\|_{L^{2}(\Omega)}$ -norm, and ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\cdot\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}$ be the $H(\operatorname{div})$ -norm

{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\bm{\tau}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{2}:=\|\bm{\tau}\|^{2}+\|\nabla\cdot\bm{\tau}\|^{2}.

In this section, we derive a priori error estimates of the mixed method (1.12) in Theorem 2.1. The same analysis implies that (1.12) admits a unique solution provided $h$ is sufficiently small.

The operator $\mathcal{A}$ is singular and satisfies

(\mathcal{A}\bm{\tau}_{1},\bm{\tau}_{2})=(\bm{\tau}_{1},\mathcal{A}\bm{\tau}_{2})=(\mathcal{A}\bm{\tau}_{1},\mathcal{A}\bm{\tau}_{2}),\quad\bm{\tau}_{1},\bm{\tau}_{2}\in[L^{2}(\Omega)]^{d\times d}.

(2.1)

In addition, it holds that (see Caiwang2007 ; ArnoldDouglasGupta1984 )

\|\bm{\tau}\|\lesssim\|\bm{\tau}\|_{\mathcal{A}}+\|\nabla\cdot\bm{\tau}\|_{-1},\quad\bm{\tau}\in\bm{\Sigma},

(2.2)

where $\|\bm{\tau}\|_{\mathcal{A}}:=(\mathcal{A}\bm{\tau},\bm{\tau})^{\frac{1}{2}}.$ The inequality (2.2) is a key ingredient in the analysis of pseudostress-based methods. Define $\Pi_{h}:[H^{1}(\Omega)]^{d\times d}\rightarrow\bm{\Sigma}_{h}$ by

\Pi_{h}\bm{\tau}:=\widetilde{\Pi}_{h}\bm{\tau}-\frac{1}{d|\Omega|}\left(\int_{\Omega}\text{Tr}\widetilde{\Pi}_{h}\bm{\tau}dx\right)\bm{I},

(2.3)

where $\widetilde{\Pi}_{h}$ is applied to each row of $\bm{\tau}.$ By abuse of notation, we may also use $P_{h}:[L^{2}(\Omega)]^{d}\rightarrow\bm{V}_{h}$ to denote the component-wise $L^{2}$ -projection. It follows from (1.11) that

\nabla\cdot\Pi_{h}\bm{\tau}=P_{h}\nabla\cdot\bm{\tau},\quad\forall\bm{\tau}\in[H^{1}(\Omega)]^{d\times d}.

(2.4)

For convenience, we introduce the interpolation errors

\bm{\rho}_{h}:=\bm{u}-P_{h}\bm{u},\quad\bm{\zeta}_{h}:=\bm{\sigma}-\Pi_{h}\bm{\sigma},

which can be easily estimated by (see CaiTongVassilevskiWang2010b )


$\displaystyle\\|\bm{v}-P_{h}\bm{v}\\|$	$\displaystyle\lesssim h^{s}\|\bm{v}\|_{s},$	(2.5a)
$\displaystyle\\|\bm{\tau}-\Pi_{h}\bm{\tau}\\|$	$\displaystyle\lesssim h^{s}\|\bm{\tau}\|_{s},$	(2.5b)
$\displaystyle\\|\nabla\cdot(\bm{\tau}-{\Pi}_{h}\bm{\tau})\\|$	$\displaystyle\lesssim h^{s}\|\operatorname{div}\bm{\tau}\|_{s},$	(2.5c)

where $0\leq s\leq r+1$ , $\int_{\Omega}\text{Tr}\bm{\tau}dx=0$ , and $\bm{v},\bm{\tau}$ satisfy the regularity indicated by the right hand sides. For the BDM element (1.8), it additionally holds that

\|\bm{\tau}-\Pi_{h}\bm{\tau}\|\lesssim h^{s}|\bm{\tau}|_{s},\quad 0\leq s\leq r+2.

(2.6)

The essential errors to be estimated are

\bm{e}_{h}:=P_{h}\bm{u}-\bm{u}_{h},\quad\bm{\xi}_{h}:=\Pi_{h}\bm{\sigma}-\bm{\sigma}_{h}.

Subtracting (1.12) from (1.6) (with $\nu=1$ ), we obtain the error equation


$\displaystyle(\mathcal{A}(\bm{\sigma}-{\bm{\sigma}_{h}}),\bm{\tau}_{h})+(\nabla\cdot\bm{\tau}_{h},\bm{u}-\bm{u}_{h})$	$\displaystyle=0,\quad\bm{\tau}_{h}\in\bm{\Sigma}_{h},$	(2.7a)
$\displaystyle-(\nabla\cdot(\bm{\sigma}-{\bm{\sigma}_{h}}),\bm{v}_{h})+(\bm{b}\cdot\mathcal{A}(\bm{\sigma}-{\bm{\sigma}_{h}}),\bm{v}_{h})+(c(\bm{u}-\bm{u}_{h}),\bm{v}_{h})$	$\displaystyle=0,\quad\bm{v}_{h}\in\bm{V}_{h}.$	(2.7b)

To estimate $\|\bm{e}_{h}\|$ , we consider the adjoint problem of (1.6): Find $\widetilde{\bm{\sigma}}\in\bm{\Sigma}$ and $\widetilde{\bm{u}}\in\bm{V}$ such that


$\displaystyle(\mathcal{A}\widetilde{\bm{\sigma}},\bm{\tau})+(\nabla\cdot\bm{\tau},\widetilde{\bm{u}})-(\bm{b}\cdot\mathcal{A}\bm{\tau},\widetilde{\bm{u}})$	$\displaystyle=0,\quad\bm{\tau}\in\bm{\Sigma},$	(2.8a)
$\displaystyle-(\nabla\cdot\widetilde{\bm{\sigma}},\bm{v})+(c\widetilde{\bm{u}},\bm{v})$	$\displaystyle=(\bm{e}_{h},\bm{v})\quad\bm{v}\in\bm{V}.$	(2.8b)

In the analysis, it is assumed that (2.8) admits the elliptic regularity

\|\widetilde{\bm{\sigma}}\|_{1}+\|\widetilde{\bm{u}}\|_{1}\lesssim\|\bm{e}_{h}\|.

(2.9)

The next lemma is a supercloseness estimate which is crucial for both a priori and superconvergence error analysis.

Lemma 2.1

It holds that

\|\bm{e}_{h}\|\lesssim h\big{(}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\bm{\sigma}-{\bm{\sigma}_{h}}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}+\|\bm{u}-\bm{u}_{h}\|+\|\bm{\rho}_{h}\|\big{)}.

Proof

Taking $\bm{v}=\bm{e}_{h}$ in (2.8b), we obtain

\|\bm{e}_{h}\|^{2}=-(\nabla\cdot\widetilde{\bm{\sigma}},\bm{e}_{h})+(c\widetilde{\bm{u}},\bm{e}_{h}).

(2.10)

Using the definition (1.10) and the property (2.4), we have

		$\displaystyle(\nabla\cdot\widetilde{\bm{\sigma}},\bm{e}_{h})=(P_{h}\nabla\cdot\widetilde{\bm{\sigma}},P_{h}\bm{u}-\bm{u}+\bm{u}-\bm{u}_{h})$		(2.11)
		$\displaystyle\quad=(P_{h}\nabla\cdot\widetilde{\bm{\sigma}},\bm{u}-\bm{u}_{h})=(\nabla\cdot\Pi_{h}\widetilde{\bm{\sigma}},\bm{u}-\bm{u}_{h}).$		(2.11)

Combining (2.7a) with $\bm{\tau}_{h}=\Pi_{h}\widetilde{\bm{\sigma}}$ implies that

		$\displaystyle(\nabla\cdot\Pi_{h}\widetilde{\bm{\sigma}},\bm{u}-\bm{u}_{h})=-(\mathcal{A}(\bm{\sigma}-{\bm{\sigma}_{h}}),\Pi_{h}\widetilde{\bm{\sigma}})$		(2.12)
		$\displaystyle\quad=(\mathcal{A}(\bm{\sigma}-{\bm{\sigma}_{h}}),\widetilde{\bm{\sigma}}-\Pi_{h}\widetilde{\bm{\sigma}})-(\mathcal{A}(\bm{\sigma}-{\bm{\sigma}_{h}}),\widetilde{\bm{\sigma}}).$		(2.12)

Using (2.1) and setting $\widetilde{\bm{\tau}}=\bm{\sigma}-{\bm{\sigma}_{h}}$ in (2.8a), we have

		$\displaystyle(\mathcal{A}(\bm{\sigma}-{\bm{\sigma}_{h}}),\widetilde{\bm{\sigma}})=(\mathcal{A}\widetilde{\bm{\sigma}},\bm{\sigma}-{\bm{\sigma}_{h}})$
		$\displaystyle=-(\nabla\cdot(\bm{\sigma}-{\bm{\sigma}_{h}}),\widetilde{\bm{u}})+(\bm{b}\cdot\mathcal{A}(\bm{\sigma}-{\bm{\sigma}_{h}}),\widetilde{\bm{u}})$
		$\displaystyle=-(\nabla\cdot(\bm{\sigma}-{\bm{\sigma}_{h}}),\widetilde{\bm{u}}-P_{h}\widetilde{\bm{u}})-(\nabla\cdot(\bm{\sigma}-{\bm{\sigma}_{h}}),P_{h}\widetilde{\bm{u}})+(\bm{b}\cdot\mathcal{A}(\bm{\sigma}-{\bm{\sigma}_{h}}),\widetilde{\bm{u}}).$

It then follows from the above equation and (2.7b) with $\bm{v}_{h}=P_{h}\widetilde{\bm{u}}$ that

		$\displaystyle(\mathcal{A}(\bm{\sigma}-{\bm{\sigma}_{h}}),\widetilde{\bm{\sigma}})=-(\nabla\cdot(\bm{\sigma}-{\bm{\sigma}_{h}}),\widetilde{\bm{u}}-P_{h}\widetilde{\bm{u}})+(\bm{b}\cdot\mathcal{A}(\bm{\sigma}-{\bm{\sigma}_{h}}),\widetilde{\bm{u}})$		(2.13)
		$\displaystyle\quad-(\bm{b}\cdot\mathcal{A}(\bm{\sigma}-{\bm{\sigma}_{h}}),P_{h}\widetilde{\bm{u}})-(c(\bm{u}-\bm{u}_{h}),P_{h}\widetilde{\bm{u}}).$		(2.13)

As a result of (2.11), (2.12) and (2.13), we obtain

	$\displaystyle(\nabla\cdot\widetilde{\bm{\sigma}},\bm{e}_{h})$	$\displaystyle=(\mathcal{A}(\bm{\sigma}-{\bm{\sigma}_{h}}),\widetilde{\bm{\sigma}}-\Pi_{h}\widetilde{\bm{\sigma}})+(\nabla\cdot(\bm{\sigma}-{\bm{\sigma}_{h}}),\widetilde{\bm{u}}-P_{h}\widetilde{\bm{u}})$		(2.14)
		$\displaystyle-(\bm{b}\cdot\mathcal{A}(\bm{\sigma}-{\bm{\sigma}_{h}}),\widetilde{\bm{u}}-P_{h}\widetilde{\bm{u}})+(c(\bm{u}-\bm{u}_{h}),P_{h}\widetilde{\bm{u}}).$		(2.14)

On the other hand, the second term on the right hand side of (2.10) is

	$\displaystyle(c\widetilde{\bm{u}},\bm{e}_{h})$	$\displaystyle=(c\widetilde{\bm{u}},P_{h}\bm{u}-\bm{u}+\bm{u}-\bm{u}_{h})$		(2.15)
		$\displaystyle=(c\widetilde{\bm{u}}-P_{h}(c\widetilde{\bm{u}}),P_{h}\bm{u}-\bm{u})+(c(\bm{u}-\bm{u}_{h}),\widetilde{\bm{u}}).$		(2.15)

Finally with (2.14) and (2.15), the error in (2.10) can be written as

		$\displaystyle\\|\bm{e}_{h}\\|^{2}=(\mathcal{A}(\bm{\sigma}-{\bm{\sigma}_{h}}),\Pi_{h}\widetilde{\bm{\sigma}}-\widetilde{\bm{\sigma}})-(\nabla\cdot(\bm{\sigma}-{\bm{\sigma}_{h}}),\widetilde{\bm{u}}-P_{h}\widetilde{\bm{u}})$		(2.16)
		$\displaystyle\quad+(\bm{b}\cdot\mathcal{A}(\bm{\sigma}-{\bm{\sigma}_{h}}),\widetilde{\bm{u}}-P_{h}\widetilde{\bm{u}})+(c(\bm{u}-\bm{u}_{h}),\widetilde{\bm{u}}-P_{h}\widetilde{\bm{u}})$
		$\displaystyle\quad+(c\widetilde{\bm{u}}-P_{h}(c\widetilde{\bm{u}}),P_{h}\bm{u}-\bm{u}).$

Using (2.16), (2.5), and the Cauchy–Schwarz inequality, we obtain

\|\bm{e}_{h}\|^{2}\lesssim h\|\widetilde{\bm{\sigma}}\|_{1}\|\bm{\sigma}-{\bm{\sigma}_{h}}\|+h\|\widetilde{\bm{u}}\|_{1}\big{(}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\bm{\sigma}-{\bm{\sigma}_{h}}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}+\|\bm{u}-\bm{u}_{h}\|+\|\bm{\rho}_{h}\|\big{)}.

Combining the above estimate with (2.9) completes the proof. ∎

Lemma 2.1 is a supercloseness estimate, i.e., $P_{h}\bm{u}$ and $\bm{u}_{h}$ are much closer than the distance predicted by standard a priori error estimates. However, a priori error estimates on $\|\bm{\sigma}-\bm{\sigma}_{h}\|$ and $\|\nabla\cdot(\bm{\sigma}-\bm{\sigma}_{h})\|$ are not known at the moment. When deriving a priori error estimates, we need the $L^{2}$ and negative norm estimates of $\nabla\cdot\bm{\xi}_{h}.$

Lemma 2.2

It holds that


$\displaystyle\\|\nabla\cdot\bm{\xi}_{h}\\|$	$\displaystyle\lesssim\\|\bm{\zeta}_{h}\\|+\\|\bm{\xi}_{h}\\|+\\|\bm{e}_{h}\\|+h\\|\bm{\rho}_{h}\\|,$	(2.17a)
$\displaystyle\\|\nabla\cdot\bm{\xi}_{h}\\|_{-1}$	$\displaystyle\lesssim h\big{(}{\left\|\kern-1.07639pt\left\|\kern-1.07639pt\left\|\bm{\zeta}_{h}\right\|\kern-1.07639pt\right\|\kern-1.07639pt\right\|}+{\left\|\kern-1.07639pt\left\|\kern-1.07639pt\left\|\bm{\xi}_{h}\right\|\kern-1.07639pt\right\|\kern-1.07639pt\right\|}+\\|\bm{\rho}_{h}\\|\big{)}+\\|\bm{e}_{h}\\|.$	(2.17b)

Lemme 2.2 is also essential for proving the supercloseness estimate in Theorem 3.1 and the proof is postponed in the end of this section. In this work, we say $h$ is sufficiently small provided $h\leq h_{0},$ where $h_{0}$ is an absolute constant relying on $\Omega$ , $\bm{b}$ , $c$ , $C_{\text{mesh}}$ , $\nu$ . Now we are in a position to present the first main result in this paper.

Theorem 2.1

For sufficiently small $h$ , it holds that


$\displaystyle\\|\bm{u}-\bm{u}_{h}\\|$	$\displaystyle\lesssim\\|\bm{\rho}_{h}\\|+h{\left\|\kern-1.07639pt\left\|\kern-1.07639pt\left\|\bm{\zeta}_{h}\right\|\kern-1.07639pt\right\|\kern-1.07639pt\right\|},$	(2.18a)
$\displaystyle{\left\|\kern-1.07639pt\left\|\kern-1.07639pt\left\|\bm{\sigma}-{\bm{\sigma}_{h}}\right\|\kern-1.07639pt\right\|\kern-1.07639pt\right\|}$	$\displaystyle\lesssim{\left\|\kern-1.07639pt\left\|\kern-1.07639pt\left\|\bm{\zeta}_{h}\right\|\kern-1.07639pt\right\|\kern-1.07639pt\right\|}+h\\|\bm{\rho}_{h}\\|,$	(2.18b)
$\displaystyle\\|\bm{\sigma}-{\bm{\sigma}_{h}}\\|$	$\displaystyle\lesssim\\|\bm{\zeta}_{h}\\|+h\\|\nabla\cdot\bm{\zeta}_{h}\\|+h\\|\bm{\rho}_{h}\\|.$	(2.18c)

Proof

Using Lemma 2.1 and the triangle inequality

	$\displaystyle\\|\bm{u}-\bm{u}_{h}\\|$	$\displaystyle\leq\\|\bm{\rho}_{h}\\|+\\|\bm{e}_{h}\\|,$
	$\displaystyle{\left\|\kern-1.07639pt\left\|\kern-1.07639pt\left\|\bm{\sigma}-\bm{\sigma}_{h}\right\|\kern-1.07639pt\right\|\kern-1.07639pt\right\|}$	$\displaystyle\leq{\left\|\kern-1.07639pt\left\|\kern-1.07639pt\left\|\bm{\zeta}_{h}\right\|\kern-1.07639pt\right\|\kern-1.07639pt\right\|}+{\left\|\kern-1.07639pt\left\|\kern-1.07639pt\left\|\bm{\xi}_{h}\right\|\kern-1.07639pt\right\|\kern-1.07639pt\right\|},$

we obtain

\|\bm{e}_{h}\|\lesssim h\big{(}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\bm{\zeta}_{h}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}+{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\bm{\xi}_{h}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}+\|\bm{\rho}_{h}\|\big{)}+h\|\bm{e}_{h}\|.

(2.19)

Plugging (2.17a) into (2.19) and kicking $h\|\bm{e}_{h}\|$ (with sufficiently small $h$ ) back to the left hand side yields

\|\bm{e}_{h}\|\lesssim h\big{(}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\bm{\zeta}_{h}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}+\|\bm{\xi}_{h}\|+\|\bm{\rho}_{h}\|\big{)}.

(2.20)

On the other hand, with help of (2.7a) with $\bm{\tau}_{h}=\bm{\xi}_{h}$ , we obtain

	$\displaystyle(\mathcal{A}\bm{\xi}_{h},\bm{\xi}_{h})$	$\displaystyle=(\mathcal{A}(\Pi_{h}\bm{\sigma}-\bm{\sigma}),\bm{\xi}_{h})-(\nabla\cdot\bm{\xi}_{h},\bm{u}-\bm{u}_{h})$		(2.21)
		$\displaystyle=(\mathcal{A}(\Pi_{h}\bm{\sigma}-\bm{\sigma}),\bm{\xi}_{h})-(\nabla\cdot\bm{\xi}_{h},P_{h}\bm{u}-\bm{u}_{h}).$		(2.21)

It follows from (2.2), (2.21), and a Young’s inequality with $\varepsilon>0$ that

	$\displaystyle\\|\bm{\xi}_{h}\\|^{2}$	$\displaystyle\lesssim(\mathcal{A}\bm{\xi}_{h},\bm{\xi}_{h})+\\|\nabla\cdot\bm{\xi}_{h}\\|_{-1}^{2}$
		$\displaystyle\lesssim\\|\bm{\zeta}_{h}\\|\\|\bm{\xi}_{h}\\|+\\|\nabla\cdot\bm{\xi}_{h}\\|\\|\bm{e}_{h}\\|+\\|\nabla\cdot\bm{\xi}_{h}\\|_{-1}^{2}$
		$\displaystyle\leq\frac{\varepsilon^{2}}{2}\\|\bm{\xi}_{h}\\|^{2}+\frac{\varepsilon^{-2}}{2}\\|\bm{\zeta}_{h}\\|^{2}+\frac{\varepsilon^{2}}{2}\\|\nabla\cdot\bm{\xi}_{h}\\|^{2}+\frac{\varepsilon^{-2}}{2}\\|\bm{e}_{h}\\|^{2}+\\|\nabla\cdot\bm{\xi}_{h}\\|_{-1}^{2},$

or equivalently

\|\bm{\xi}_{h}\|\leq C^{*}\big{(}{\varepsilon}\|\bm{\xi}_{h}\|+{\varepsilon^{-1}}\|\bm{\zeta}_{h}\|+{\varepsilon}\|\nabla\cdot\bm{\xi}_{h}\|+{\varepsilon^{-1}}\|\bm{e}_{h}\|+\|\nabla\cdot\bm{\xi}_{h}\|_{-1}\big{)},

(2.22)

where $C^{*}$ is independent of $\varepsilon$ . Using (2.22), Lemma 2.2, and the triangle inequality, we deduce that

	$\displaystyle\\|\bm{\xi}_{h}\\|$	$\displaystyle\leq C^{*}\big{\{}\varepsilon\big{(}\\|\bm{\xi}_{h}\\|+\\|\bm{\zeta}_{h}\\|+\\|\bm{e}_{h}\\|+h\\|\bm{\rho}_{h}\\|\big{)}$
		$\displaystyle+{\varepsilon^{-1}}\big{(}\\|\bm{\zeta}_{h}\\|+\\|\bm{e}_{h}\\|\big{)}+h{\left\|\kern-1.07639pt\left\|\kern-1.07639pt\left\|\bm{\zeta}_{h}\right\|\kern-1.07639pt\right\|\kern-1.07639pt\right\|}+h\\|\bm{\xi}_{h}\\|+h\\|\bm{\rho}_{h}\\|+\\|\bm{e}_{h}\\|\big{\}}.$

In the above estimate, it suffices to choose sufficiently small $\varepsilon$ and $h$ to obtain

\|\bm{\xi}_{h}\|\leq C\big{(}\|\bm{\zeta}_{h}\|+h\|\nabla\cdot\bm{\zeta}_{h}\|+\|\bm{e}_{h}\|+h\|\bm{\rho}_{h}\|\big{)}.

(2.23)

Plugging (2.23) into (2.20) with sufficiently small $h$ then leads to

\|\bm{e}_{h}\|\leq Ch\big{(}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\bm{\zeta}_{h}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}+\|\bm{\rho}_{h}\|\big{)}.

(2.24)

Therefore we close the loop. As a consequence of (2.23), (2.24), (2.17a), we obtain

	$\displaystyle\\|\bm{\xi}_{h}\\|$	$\displaystyle\leq C\big{(}\\|\bm{\zeta}_{h}\\|+h\\|\nabla\cdot\bm{\zeta}_{h}\\|+h\\|\bm{\rho}_{h}\\|\big{)},$		(2.25)
	$\displaystyle\\|\nabla\cdot\bm{\xi}_{h}\\|$	$\displaystyle\leq C\big{(}\\|\bm{\zeta}_{h}\\|+h\\|\nabla\cdot\bm{\zeta}_{h}\\|+h\\|\bm{\rho}_{h}\\|\big{)}.$		(2.25)

Theorem 2.1 eventually follows from (2.24), (2.25) and a triangle inequality. ∎

For sufficiently smooth $(\bm{\sigma},\bm{u})$ , Theorem 2.1 with (2.5) yields


$\displaystyle\\|\bm{u}-\bm{u}_{h}\\|$	$\displaystyle\lesssim h^{r+1}\big{(}\|\bm{u}\|_{r+1}+h\|\bm{\sigma}\|_{r+1}+h\|\nabla\cdot\bm{\sigma}\|_{r+1}\big{)},$	(2.26a)
$\displaystyle{\left\|\kern-1.07639pt\left\|\kern-1.07639pt\left\|\bm{\sigma}-{\bm{\sigma}_{h}}\right\|\kern-1.07639pt\right\|\kern-1.07639pt\right\|}$	$\displaystyle\lesssim h^{r+1}\big{(}\|\bm{\sigma}\|_{r+1}+\|\nabla\cdot\bm{\sigma}\|_{r+1}+h\|\bm{u}\|_{r+1}\big{)}.$	(2.26b)

When $\bm{\Sigma}_{h}$ is based on the BDM element (1.8), the following improved error estimate follows from (2.18c) and (2.6).

\|\bm{\sigma}-\bm{\sigma}_{h}\|\lesssim h^{r+2}\big{(}|\bm{\sigma}|_{r+2}+|\bm{u}|_{r+1}\big{)}.

(2.27)

The well-posedness (1.12) follows from the same analysis in the proof of Theorem 2.1.

Corollary 1

For sufficiently small $h$ , the mixed method (1.12) has a unique solution.

Proof

To establish the existence and uniqueness of the solution to (1.12), it suffices to show the uniqueness because of linearity. Suppose $(\bm{\sigma}_{h},\bm{u}_{h})$ and $(\widehat{\bm{\sigma}}_{h},\widehat{\bm{u}}_{h})$ are both solutions to (1.12), then we have the error equation

	$\displaystyle(\mathcal{A}\bm{e}_{\sigma},\bm{\tau}_{h})+(\nabla\cdot\bm{\tau}_{h},\bm{e}_{u})$	$\displaystyle=0,\quad\bm{\tau}_{h}\in\bm{\Sigma}_{h},$		(2.28)
	$\displaystyle-(\nabla\cdot\bm{e}_{\sigma},\bm{v}_{h})+(\bm{b}\cdot\mathcal{A}\bm{e}_{\sigma},\bm{v}_{h})+(c\bm{e}_{u},\bm{v}_{h})$	$\displaystyle=0,\quad\bm{v}_{h}\in\bm{V}_{h},$		(2.28)

where $\bm{e}_{\sigma}:=\bm{\sigma}_{h}-\widehat{\bm{\sigma}}_{h}$ , $\bm{e}_{u}:=\bm{u}_{h}-\widehat{\bm{u}}_{h}.$ Consider the dual problem

	$\displaystyle(\mathcal{A}\widetilde{\bm{\sigma}},\bm{\tau})+(\nabla\cdot\bm{\tau},\widetilde{\bm{u}})-(\bm{b}\cdot\mathcal{A}\bm{\tau},\widetilde{\bm{u}})$	$\displaystyle=0,\quad\bm{\tau}\in\bm{\Sigma},$		(2.29)
	$\displaystyle-(\nabla\cdot\widetilde{\bm{\sigma}},\bm{v})+(c\widetilde{\bm{u}},\bm{v})$	$\displaystyle=(\bm{e}_{u},v),\quad\bm{v}\in\bm{V}.$		(2.29)

It then follows from (2.28), (2.29) and the same analysis in Lemma 2.1 that

\|\bm{e}_{u}\|\lesssim h{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\bm{e}_{\sigma}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|},

(2.30)

provided $h$ is sufficiently small. The same argument for proving Lemma 2.2 yields

	$\displaystyle\\|\nabla\cdot\bm{e}_{\sigma}\\|$	$\displaystyle\lesssim\\|\bm{e}_{\sigma}\\|+\\|\bm{e}_{u}\\|,$		(2.31)
	$\displaystyle\\|\nabla\cdot\bm{e}_{\sigma}\\|_{-1}$	$\displaystyle\lesssim h{\left\|\kern-1.07639pt\left\|\kern-1.07639pt\left\|\bm{e}_{\sigma}\right\|\kern-1.07639pt\right\|\kern-1.07639pt\right\|}+\\|\bm{e}_{u}\\|.$		(2.31)

Following the proof of (2.23) and using (2.30), (2.31), we obtain

\|\bm{e}_{\sigma}\|\lesssim\|\bm{e}_{u}\|,

(2.32)

when $h$ is small enough. Finally combining (2.30), (2.32), (2.31) yields

\|\bm{e}_{u}\|\lesssim h\|\bm{e}_{u}\|.

Hence $\bm{e}_{u}=\bm{0}$ and $\bm{e}_{\sigma}=\bm{0}$ provided $h$ is sufficiently small. ∎

Proof (Proof of Lemma 2.2)

Let $\bm{v}_{h}=\nabla\cdot\bm{\xi}_{h}/\|\nabla\cdot\bm{\xi}_{h}\|$ . Using (2.4) and (2.7b), we arrive at

		$\displaystyle\\|\nabla\cdot\bm{\xi}_{h}\\|=(\nabla\cdot\bm{\xi}_{h},\bm{v}_{h})=(P_{h}\nabla\cdot\bm{\sigma}-\nabla\cdot{\bm{\sigma}_{h}},\bm{v}_{h})=(\nabla\cdot(\bm{\sigma}-{\bm{\sigma}_{h}}),\bm{v}_{h})$		(2.33)
		$\displaystyle\quad=(\bm{b}\cdot\mathcal{A}(\bm{\sigma}-{\bm{\sigma}_{h}}),\bm{v}_{h})+(c(P_{h}\bm{u}-\bm{u}_{h}),\bm{v}_{h})+((c-c_{h})(\bm{u}-P_{h}\bm{u}),\bm{v}_{h}).$		(2.33)

where $c_{h}$ is the piecewise average of $c$ w.r.t. $\mathcal{T}_{h}.$ The estimate (2.17a) then follows from (2.33), the Cauchy–Schwarz and triangle inequalities.

To prove (2.17b), let $\bm{v}\in[H_{0}^{1}(\Omega)]^{d}$ with $\|\bm{v}\|_{1}=1$ . It follows from (2.4), (1.10), and (2.7b) that

		$\displaystyle(\nabla\cdot\bm{\xi}_{h},\bm{v})=(P_{h}\nabla\cdot\bm{\sigma}-\nabla\cdot\bm{\sigma},\bm{v})+(\nabla\cdot\bm{\sigma}-\nabla\cdot\bm{\sigma}_{h},\bm{v})$		(2.34)
		$\displaystyle\quad=(P_{h}\nabla\cdot\bm{\sigma}-\nabla\cdot\bm{\sigma},\bm{v}-P_{h}\bm{v})+(\nabla\cdot\bm{\sigma}-\nabla\cdot\bm{\sigma}_{h},\bm{v}-P_{h}\bm{v})$
		$\displaystyle\qquad+(\bm{b}\cdot\mathcal{A}(\bm{\sigma}-\bm{\sigma}_{h}),\bm{v})+(\bm{b}\cdot\mathcal{A}(\bm{\sigma}-\bm{\sigma}_{h}),P_{h}\bm{v}-\bm{v})$
		$\displaystyle\qquad+(c(\bm{u}-P_{h}\bm{u}),P_{h}\bm{v})+(c(P_{h}\bm{u}-\bm{u}_{h}),P_{h}\bm{v}).$

Using (2.1) and (2.7a), we can rewrite $(\bm{b}\cdot\mathcal{A}(\bm{\sigma}-\bm{\sigma}_{h}),\bm{v})$ as

		$\displaystyle\quad(\bm{b}\cdot\mathcal{A}(\bm{\sigma}-\bm{\sigma}_{h}),\bm{v})$		(2.35)
		$\displaystyle=(\mathcal{A}(\bm{\sigma}-\bm{\sigma}_{h}),\bm{v}\bm{b}^{\top})=(\mathcal{A}(\bm{\sigma}-\bm{\sigma}_{h}),\mathcal{A}(\bm{v}\bm{b}^{\top}))$
		$\displaystyle=(\mathcal{A}(\bm{\sigma}-\bm{\sigma}_{h}),\Pi_{h}\mathcal{A}(\bm{v}\bm{b}^{\top}))+(\mathcal{A}(\bm{\sigma}-\bm{\sigma}_{h}),(\bm{I}-\Pi_{h})\mathcal{A}(\bm{v}\bm{b}^{\top}))$
		$\displaystyle=-(\nabla\cdot\Pi_{h}\mathcal{A}(\bm{v}\bm{b}^{\top}),\bm{u}-\bm{u}_{h})+(\mathcal{A}(\bm{\sigma}-\bm{\sigma}_{h}),(\bm{I}-\Pi_{h})\mathcal{A}(\bm{v}\bm{b}^{\top}))$
		$\displaystyle=-(\nabla\cdot\Pi_{h}\mathcal{A}(\bm{v}\bm{b}^{\top}),P_{h}\bm{u}-\bm{u}_{h})+(\mathcal{A}(\bm{\sigma}-\bm{\sigma}_{h}),(\bm{I}-\Pi_{h})\mathcal{A}(\bm{v}\bm{b}^{\top})).$

On the other hand,

(c(\bm{u}-P_{h}\bm{u}),P_{h}\bm{v})=((c-c_{h})(\bm{u}-P_{h}\bm{u}),P_{h}\bm{v}).

(2.36)

Combining (2.34) with (2.35) and (2.36) and using (2.5), we obtain

		$\displaystyle\\|\nabla\cdot\bm{\xi}_{h}\\|_{-1}=\sup_{\\|\bm{v}\\|_{1}=1}(\nabla\cdot\bm{\xi}_{h},\bm{v})$		(2.37)
		$\displaystyle\lesssim h\big{(}\\|P_{h}\nabla\cdot\bm{\sigma}-\nabla\cdot\bm{\sigma}\\|+{\left\|\kern-1.07639pt\left\|\kern-1.07639pt\left\|\bm{\sigma}-\bm{\sigma}_{h}\right\|\kern-1.07639pt\right\|\kern-1.07639pt\right\|}+\\|\bm{u}-P_{h}\bm{u}\\|\big{)}$
		$\displaystyle\quad+(1+\\|\nabla\cdot\Pi_{h}\mathcal{A}(\bm{v}\bm{b}^{\top})\\|)\\|\bm{e}_{h}\\|.$

It follows from (2.4) and elementary calculation that

\|\nabla\cdot\Pi_{h}\mathcal{A}(\bm{v}\bm{b}^{\top})\|=\|P_{h}\nabla\cdot\mathcal{A}(\bm{v}\bm{b}^{\top})\|\leq\|\nabla\cdot\mathcal{A}(\bm{v}\bm{b}^{\top})\|\leq C\|\bm{v}\|_{1}\leq C.

Therefore using the previous estimate and (2.37), we obtain

\|\nabla\cdot\bm{\xi}\|_{-1}\lesssim h\|\nabla\cdot\Pi_{h}\bm{\sigma}-\nabla\cdot\bm{\sigma}\|+h{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\bm{\sigma}-\bm{\sigma}_{h}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}+h\|\bm{\rho}_{h}\|+\|\bm{e}_{h}\|.

The proof is complete. ∎

3 Superconvergence on pseudostress

Using Lemma 2.1 and (2.26), we obtain the supercloseness estimate on $\bm{e}_{h}=P_{h}\bm{u}-\bm{u}_{h}$ .

\|\bm{e}_{h}\|\lesssim h^{r+2}\big{(}|\bm{\sigma}|_{r+1}+|\nabla\cdot\bm{\sigma}|_{r+1}+|\bm{u}|_{r+1}\big{)}.

(3.1)

In this section, we shall prove an improved estimate for $\bm{\xi}_{h}=\Pi_{h}\bm{\sigma}-\bm{\sigma}_{h}$ . To this end, $\bm{\Sigma}_{h}$ is endowed with the $\Lambda$ -inner product

\Lambda(\bm{\tau}_{1},\bm{\tau}_{2}):=(\mathcal{A}\bm{\tau}_{1},\bm{\tau}_{2})+(\nabla\cdot\bm{\tau}_{1},\nabla\cdot\bm{\tau}_{2}),\quad\bm{\tau}_{1},\bm{\tau}_{2}\in\bm{\Sigma}_{h}.

Let $\bm{Z}_{h}:=\ker\big{(}\nabla\cdot|_{\bm{\Sigma}_{h}}\big{)}$ and $\bm{Z}^{\perp}_{h}$ be the orthogonal complement of $\bm{Z}_{h}$ w.r.t. $\Lambda(\cdot,\cdot)$ in $\bm{\Sigma}_{h}$ . We shall decompose $\bm{\bm{\xi}}_{h}\in\bm{\Sigma}_{h}$ by the discrete Helmholtz decomposition

\bm{\Sigma}_{h}=\bm{Z}_{h}\oplus_{\Lambda}\bm{Z}_{h}^{\perp},

(3.2)

which can be analyzed using tools in finite element exterior calculus (FEEC), see, e.g., ArnoldFalkWinther2006 ; ArnoldFalkWinther2010 . In the theory of FEEC, an essential ingredient is the bounded projections that commute with exterior derivatives. In particular, there exist projections $\widetilde{\pi}_{h}:[L^{2}(\Omega)]^{d}\rightarrow\widetilde{\bm{\Sigma}}_{h}$ and ${Q}_{h}:L^{2}(\Omega)\rightarrow{V}_{h}$ (see, e.g., ArnoldFalkWinther2006 ; ChristiansenWinther2008 ) such that

		$\displaystyle\widetilde{\pi}_{h}\|_{\widetilde{\bm{\Sigma}}_{h}}=\text{Id},\quad{Q}_{h}\|_{{V}_{h}}=\text{Id},$		(3.3)
		$\displaystyle\nabla\cdot\widetilde{\pi}_{h}\bm{\tau}={Q}_{h}\nabla\cdot\bm{\tau}\text{ for }\bm{\tau}\in H(\text{div},\Omega),$
		$\displaystyle\\|\widetilde{\pi}_{h}\bm{\tau}\\|\lesssim\\|\bm{\tau}\\|\text{ for }\bm{\tau}\in[L^{2}(\Omega)]^{d},$

where Id is the identity mapping. Starting from $\widetilde{\pi}_{h},Q_{h},$ one can easily obtain commuting bounded projections onto $\bm{\Sigma}_{h}$ and $\bm{V}_{h}$ . Let

\pi_{h}\bm{\tau}:=\widetilde{\pi}_{h}\bm{\tau}-\frac{1}{d|\Omega|}\left(\int_{\Omega}\text{Tr}(\widetilde{\pi}_{h}\bm{\tau})dx\right)\bm{I},

where $\widetilde{\pi}_{h}$ is applied to each row of $\bm{\tau}$ . Similarly, $Q_{h}$ can applied to each component of a vector-valued function in $[L^{2}(\Omega)]^{d}$ . Using the property of $\widetilde{\pi}_{h}$ in (3.3), we have


$\displaystyle\pi_{h}\|_{\bm{\Sigma}_{h}}$	$\displaystyle=\text{Id},\quad{Q}_{h}\|_{\bm{V}_{h}}=\text{Id},$	(3.4a)
$\displaystyle\nabla\cdot{\pi}_{h}\bm{\tau}$	$\displaystyle={Q}_{h}\nabla\cdot\bm{\tau}\quad\text{ for }\bm{\tau}\in\bm{\Sigma},$	(3.4b)
$\displaystyle\\|{\pi}_{h}\bm{\tau}\\|$	$\displaystyle\lesssim\\|\bm{\tau}\\|\quad\text{ for }\bm{\tau}\in[L^{2}(\Omega)]^{d},$	(3.4c)

where Id is the identity mapping. Then with the help of ${\pi}_{h}$ , $Q_{h}$ , we obtain the following lemma for estimating the decomposition component living in $\bm{Z}_{h}^{\perp}.$

Lemma 3.1

It holds that

\|\bm{\tau}_{h}\|_{\mathcal{A}}\lesssim\|\nabla\cdot\bm{\tau}_{h}\|_{-1},\quad\bm{\tau}_{h}\in\bm{Z}_{h}^{\perp}.

Proof

Let $\bm{\phi}\in[H_{0}^{1}(\Omega)]^{d}$ be the weak solution to

\Delta\bm{\phi}=\nabla\cdot\bm{\tau}_{h}.

Then $\bm{\tau}:=\nabla\bm{\phi}$ satisfies $\nabla\cdot\bm{\tau}=\nabla\cdot\bm{\tau}_{h}.$ Using it and (3.4b), (3.4a), we obtain

	$\displaystyle\nabla\cdot(\bm{\tau}_{h}-\pi_{h}\bm{\tau})=\nabla\cdot\bm{\tau}_{h}-{Q}_{h}\nabla\cdot\bm{\tau}$
	$\displaystyle\quad=\nabla\cdot\bm{\tau}_{h}-{Q}_{h}\nabla\cdot\bm{\tau}_{h}=\bm{0},$

i.e., $\bm{\tau}_{h}-\pi_{h}\bm{\tau}\in\bm{Z}_{h}$ . Recall that $\|\cdot\|_{-1}$ denotes the $H^{-1}(\Omega)$ -norm. The natural $H^{1}$ -elliptic regularity of $\Delta\bm{\phi}=\nabla\cdot\bm{\tau}_{h}$ yields

\|\bm{\tau}\|\leq\|\bm{\phi}\|_{1}\lesssim\|\nabla\cdot\bm{\tau}_{h}\|_{-1}.

(3.5)

Hence using $\bm{\tau}_{h}-\pi_{h}\bm{\tau}\in\bm{Z}_{h}$ , $\bm{\tau}_{h}\in\bm{Z}_{h}^{\perp}$ , (3.5), and (3.4c), we have

	$\displaystyle\\|\bm{\tau}_{h}\\|_{\mathcal{A}}^{2}$	$\displaystyle=\Lambda(\bm{\tau}_{h}-\pi_{h}\bm{\tau},\bm{\tau}_{h})+(\mathcal{A}\pi_{h}\bm{\tau},\bm{\tau}_{h})$
		$\displaystyle=(\mathcal{A}\pi_{h}\bm{\tau},\bm{\tau}_{h})\leq\\|\pi_{h}\bm{\tau}\\|_{\mathcal{A}}\\|\bm{\tau}_{h}\\|_{\mathcal{A}}$
		$\displaystyle\lesssim\\|\bm{\tau}\\|\\|\bm{\tau}_{h}\\|_{\mathcal{A}}\lesssim\\|\nabla\cdot\bm{\tau}_{h}\\|_{-1}\\|\bm{\tau}_{h}\\|_{\mathcal{A}}.$

The proof is complete.∎

Remark 3.1

It follows from the discrete abstract Poincaré inequality (see ArnoldFalkWinther2010 , Theorem 3.6) that

\|\bm{\tau}_{h}\|_{\mathcal{A}}\lesssim\|\nabla\cdot\bm{\tau}_{h}\|\text{ for }\bm{\tau}_{h}\in\bm{Z}_{h}^{\perp}.

Hence Lemma 3.1 is an improved discrete Poincaré inequality. The improvement is achieved by utilizing the special property of the divergence operator.

We emphasize that the improved estimate of $\bm{\xi}_{h}$ is dependent on mesh structure, type of finite elements, and quite technical. For simplicity of presentation, first let $\bm{\Sigma}_{h}\times\bm{V}_{h}$ be based on the lowest order RT element (1.7) with $r=0$ in $\mathbb{R}^{2}$ although the estimates can be generalized in several ways, see the end of this section. Given a scalar-valued function $v$ and a vector-valued function $\bm{v}=(v_{1},v_{2})^{\top}$ , let

\displaystyle\nabla^{\perp}v:=(-\partial_{x_{2}}{v},\partial_{x_{1}}{v})^{\top},\quad\bm{v}^{\perp}:=(-{v_{2}},v_{1})^{\top}.

The row-wise rotational gradient or curl is defined as

\nabla^{\perp}\bm{v}=\left[\begin{matrix}(\nabla^{\perp}v_{1})^{\top}\\ (\nabla^{\perp}v_{2})^{\top}\end{matrix}\right].

Introducing the following nodal element spaces

	$\displaystyle{W}_{h}$	$\displaystyle=\{w\in{H}^{1}(\Omega):{w}\|_{K}\in\mathcal{P}_{1}(K)\text{ for }K\in\mathcal{T}_{h}\},$
	$\displaystyle\bm{W}_{h}$	$\displaystyle=\{\bm{w}\in[W_{h}]^{2}:~{}\int_{\Omega}\nabla\cdot\bm{w}^{\perp}dx=0\},$

we obtain a two-dimensional discrete sequence

\displaystyle\bm{W}_{h}\xrightarrow{\nabla^{\perp}}\bm{\Sigma}_{h}\xrightarrow{\nabla\cdot}\bm{V}_{h}\rightarrow 0.

(3.6)

The inclusion $\nabla^{\perp}\bm{W}_{h}\subset\bm{\Sigma}_{h}$ follows from $\int_{\Omega}\nabla\cdot\bm{w}_{h}^{\perp}dx=-\int_{\Omega}\text{Tr}\nabla^{\perp}\bm{w}_{h}dx=0$ $\forall\bm{w}_{h}\in\bm{W}_{h}$ .

The supercloseness estimate of $\bm{\bm{\xi}}_{h}$ does not hold on general unstructured grids. For elliptic equations discretized by RT elements, the author Y. Li et al. in Li2018SINUM ; BankLi2019 derived several supercloseness estimates on the vector unknown in mixed methods under certain mildly structured grids. Readers are referred to Duran1990 ; EwingLiuWang1999 ; Brandts1994 ; Brandts2000 for similar results on Poisson’s equation under rectangular, $h^{2}$ -uniform quadrilateral, and uniform triangular grids. To avoid lengthy descriptions of different mesh structures, we focus on the following piecewise uniform grids.

Refer to caption — Figure 1: Piecewise uniform grid on a square

Definition 3.1

We say $\mathcal{T}_{h}$ is a uniform grid provided each pair of adjacent triangles (two triangles sharing a common edge) form an exact parallelogram. Let $\{\Omega_{i}\}_{i=1}^{N}$ be a fixed polygonal partition of $\Omega.$ We say $\mathcal{T}_{h}$ is a piecewise uniform grid provided $\mathcal{T}_{h}$ is aligned with $\{\Omega_{i}\}_{i=1}^{N}$ and $\mathcal{T}_{h}|_{\Omega_{i}}$ is a uniform grid for each $1\leq i\leq N$ .

For instance, let $\Omega=[0,1]^{2}$ be the unit square that is split into the 19 triangles $\{\Omega_{i}\}_{i=1}^{19}$ given in Figure 1 (left). After 3 consecutive uniform quad-refinement, we obtain the triangulation in Figure 1 (right), which is a piecewise uniform grid (w.r.t. $\{\Omega_{i}\}_{i=1}^{19}$ ). The piecewise uniform mesh structure is only used in the next technical lemma.

Lemma 3.2

Let $\bm{D}\in\mathbb{R}^{2\times 2}$ be a constant matrix. Let $\widetilde{\Pi}_{h}$ be the canonical interpolation onto $\widetilde{\bm{\Sigma}}_{h}$ , the lowest order RT element space. For a piecewise uniform $\mathcal{T}_{h}$ and $\bm{v}\in W^{2}_{\infty}(\Omega)$ , $w_{h}\in W_{h},$ it holds that

(\bm{\bm{D}}(\bm{v}-\widetilde{\Pi}_{h}\bm{v}),\nabla^{\perp}w_{h})\lesssim h^{2}|\log h|^{\frac{1}{2}}\|\bm{v}\|_{W^{2}_{\infty}(\Omega)}\|\nabla^{\perp}w_{h}\|.

In Li2018SINUM , one of the author proved Lemma 3.2 on mildly structured grids when $\bm{D}=\bm{I}$ . For the general anisotropic $\bm{D}$ , we could not find a complete proof in literature. In the appendix, we give a detailed proof of Lemma 3.2 using the technique in Theorem 3.2 of Brandts1994 and Lemma 3.4 of HuMaMa2021 . Now we present a supercloseness estimate for $\bm{\xi}_{h}$ , the second main result, in the next theorem.

Theorem 3.1

Let $\Omega$ be simply-connected and $\mathcal{T}_{h}$ be a piecewise uniform grid. For (1.12) using the lowest order RT element and sufficiently small $h$ , it holds that

\|\bm{\bm{\xi}}_{h}\|\lesssim h^{2}|\log h|^{\frac{1}{2}}(\|\bm{\sigma}\|_{W^{2}_{\infty}(\Omega)}+\|\bm{u}\|_{2}).

Proof

When $\Omega$ is simply connected, (3.6) is an exact sequence, i.e., $\nabla^{\perp}\bm{W}_{h}=\bm{Z}_{h}.$ It then follows from the exactness and the discrete Helmholtz decomposition (3.2) that

\displaystyle\bm{\xi}_{h}=\nabla^{\perp}\bm{\phi}_{h}\oplus_{\Lambda}\bm{\eta}_{h},

(3.7)

where $\bm{\phi}_{h}\in\bm{W}_{h},$ $\bm{\eta}_{h}\in\bm{Z}_{h}^{\perp}.$ Using Lemma 3.1 and $\nabla\cdot\nabla^{\perp}\bm{\phi}_{h}=\bm{0}$ , we obtain

\|\bm{\eta}_{h}\|_{\mathcal{A}}\lesssim\|\nabla\cdot\bm{\eta}_{h}\|_{-1}=\|\nabla\cdot\bm{\xi}_{h}\|_{-1}.

(3.8)

It remains to estimate $\|\nabla^{\perp}\bm{\phi}_{h}\|_{\mathcal{A}}$ . Note that

(\mathcal{A}\nabla^{\perp}\bm{\phi}_{h},\bm{\eta}_{h})=\Lambda(\nabla^{\perp}\bm{\phi}_{h},\bm{\eta}_{h})=0.

Then using the orthogonality, (2.7a), and (2.1), we have

		$\displaystyle\\|\nabla^{\perp}\bm{\phi}_{h}\\|_{\mathcal{A}}^{2}=(\mathcal{A}(\Pi_{h}\bm{\sigma}-\bm{\sigma}_{h}),\nabla^{\perp}\bm{\phi}_{h})$
		$\displaystyle=(\Pi_{h}\bm{\sigma}-\bm{\sigma},\mathcal{A}\nabla^{\perp}\bm{\phi}_{h})=(\widetilde{\Pi}_{h}\bm{\sigma}-\bm{\sigma},\mathcal{A}\nabla^{\perp}\bm{\phi}_{h})$
		$\displaystyle=-(\bm{\sigma}-\widetilde{\Pi}_{h}\bm{\sigma},\nabla^{\perp}\bm{\phi}_{h})-\frac{1}{2}(\bm{\sigma}-\widetilde{\Pi}_{h}\bm{\sigma},(\nabla\cdot\bm{\phi}^{\perp}_{h})\bm{I}).$

Recall that $\bm{\sigma}_{i}$ denotes the $i$ -th row of $\bm{\sigma}.$ Let $\bm{e}_{1}=(1,0)$ , $\bm{e}_{2}=(0,1)$ , and $\bm{\phi}_{h}=(\phi_{h,1},\phi_{h,2})^{\top}$ . It follows from the previous equation and direct calculation that

		$\displaystyle\\|\nabla^{\perp}\bm{\phi}_{h}\\|_{\mathcal{A}}^{2}=\sum_{i=1}^{2}\left\{-(\bm{\sigma}_{i}-\widetilde{\Pi}_{h}\bm{\sigma}_{i},\nabla^{\perp}{\phi}_{h,i})\right.$		(3.9)
		$\displaystyle\quad\left.+\frac{1}{2}(\bm{\sigma}_{i}-\widetilde{\Pi}_{h}\bm{\sigma}_{i},(\partial_{x_{1}}{\phi}_{h,2})\bm{e}_{i})-\frac{1}{2}(\bm{\sigma}_{i}-\widetilde{\Pi}_{h}\bm{\sigma}_{i},(\partial_{x_{2}}{\phi}_{h,1})\bm{e}_{i})\right\}$
		$\displaystyle\quad=\sum_{i,j=1}^{2}(\bm{D}_{ij}(\bm{\sigma}_{i}-\widetilde{\Pi}_{h}\bm{\sigma}_{i})^{\top},\nabla^{\perp}\phi_{h,j}),$

where

\bm{D}_{11}=\begin{pmatrix}-\frac{1}{2}&0\\ 0&-1\end{pmatrix},~{}\bm{D}_{12}=\begin{pmatrix}0&0\\ \frac{1}{2}&0\end{pmatrix},~{}\bm{D}_{21}=\begin{pmatrix}0&\frac{1}{2}\\ 0&0\end{pmatrix},~{}\bm{D}_{22}=\begin{pmatrix}-1&0\\ 0&-\frac{1}{2}\end{pmatrix}.

Then using (3.9), Lemma 3.2 and (2.2) with $\bm{\tau}=\nabla^{\perp}\bm{\phi}_{h}$ , we have

	$\displaystyle\\|\nabla^{\perp}\bm{\phi}_{h}\\|_{\mathcal{A}}^{2}$	$\displaystyle\lesssim h^{2}\|\log h\|^{\frac{1}{2}}\\|\bm{\sigma}\\|_{W^{2}_{\infty}(\Omega)}\\|\nabla^{\perp}\bm{\phi}_{h}\\|$		(3.10)
		$\displaystyle\lesssim h^{2}\|\log h\|^{\frac{1}{2}}\\|\bm{\sigma}\\|_{W^{2}_{\infty}(\Omega)}\\|\nabla^{\perp}\bm{\phi}_{h}\\|_{\mathcal{A}}.$		(3.10)

Combining (2.2) with (3.8), (3.10), (2.17b), (3.1), we obtain

	$\displaystyle\\|\bm{\xi}_{h}\\|$	$\displaystyle\lesssim\\|\bm{\xi}_{h}\\|_{\mathcal{A}}+\\|\nabla\cdot\bm{\xi}_{h}\\|_{-1}$
		$\displaystyle\lesssim\\|\nabla^{\perp}\bm{\phi}_{h}\\|_{\mathcal{A}}+\\|\bm{\eta}_{h}\\|_{\mathcal{A}}+\\|\nabla\cdot\bm{\xi}_{h}\\|_{-1}$
		$\displaystyle\lesssim h^{2}\|\log h\|^{\frac{1}{2}}\big{(}\\|\bm{\sigma}\\|_{W^{2}_{\infty}(\Omega)}+\\|\bm{u}\\|_{2}\big{)}.$

The proof is complete. ∎

The supercloseness estimate on $\|\bm{\xi}_{h}\|$ can be generalized on rectangular meshes. Throughout the rest of this section, let $\mathcal{T}_{h}$ be a rectangular mesh,

	$\displaystyle{W}_{h}$	$\displaystyle=\left\{w_{h}\in{H}^{1}(\Omega):{w}_{h}\|_{K}\in\mathcal{Q}_{r+1}(K)\text{ for }K\in\mathcal{T}_{h}\right\},$
	$\displaystyle\bm{W}_{h}$	$\displaystyle=\left\{\bm{w}_{h}\in[W_{h}]^{2}:~{}\int_{\Omega}\nabla^{\perp}\bm{w}_{h}dx=0\right\},$

and $\bm{\Sigma}_{h}\times\bm{V}_{h}$ be based on the rectangular RT element (1.9) with $d=2$ . When $\Omega$ is simply-connected, we still have the discrete exact sequence (3.6). Similarly to the case of triangular grids, we need the uniform mesh structure.

Definition 3.2

We say $\mathcal{T}_{h}$ is a uniformly rectangular mesh if all the rectangles of $\mathcal{T}_{h}$ are of the same shape and size.

For $w_{h}\in{W}_{h}$ with $(\nabla^{\perp}w_{h})\cdot\bm{n}|_{\partial\Omega}=0$ , Theorem 5.1 of EwingLiuWang1999 implies

(\bm{\bm{D}}(\bm{v}-\widetilde{\Pi}_{h}\bm{v}),\nabla^{\perp}w_{h})\lesssim h^{r+2}\|\bm{v}\|_{r+2}\|\nabla^{\perp}w_{h}\|,

where $\bm{v}\in[H^{r+2}(\Omega)]^{2}$ . The previous estimate is not true when $(\nabla^{\perp}w_{h})\cdot\bm{n}\neq 0$ on $\partial\Omega$ . As claimed in the remark in Example 6.2 of EwingLiuWang1999 , it holds that

(\bm{\bm{D}}(\bm{v}-\widetilde{\Pi}_{h}\bm{v}),\nabla^{\perp}w_{h})\leq Ch^{r+1.5}\|\nabla^{\perp}w_{h}\|,\quad w_{h}\in W_{h}

(3.11)

for some absolute constant $C$ independent of $h$ . Following the same argument in the proof of Theorem 3.1 and replacing Lemma 3.2 by (3.11), we actually obtain the supercloseness estimate

\|\Pi_{h}\bm{\sigma}-\bm{\sigma}_{h}\|=\mathcal{O}(h^{r+1.5}).

Similar estimates may hold on $h^{2}$ -uniform quadrilateral grids described in EwingLiuWang1999 .

4 Postprocessing

We have derived supercloseness estimates in Lemma 2.1 and Theorem 3.1. However, those results can not be directly used because $P_{h}\bm{u}$ and $\Pi_{h}\bm{\sigma}$ are not known in practice. To extract superconvergence information from the smallness of $\bm{e}_{h}$ and $\bm{\xi}_{h}$ , one may design easy-to-compute postprocessed solutions $\bm{u}_{h}^{*}$ and $\bm{\sigma}_{h}^{*}$ . For instance, following the idea of Stenberg1991 , a element-by-element postprocessing procedure for $\bm{u}_{h}$ in the pseudostress-velocity formulation of the Stokes equation is proposed in CKP2011 . In particular, let

\mathcal{P}_{r+1}(K)/\mathcal{P}_{r}(K):=\{\bm{v}_{h}\in\mathcal{P}_{r+1}(K):(\bm{v}_{h},\bm{w}_{h})_{K}=0\text{ for all }\bm{w}_{h}\in\mathcal{P}_{r}(K)\}.

where $(\cdot,\cdot)_{K}$ is the $L^{2}(K)$ -inner product. Note that the pressure could be recovered as

p_{h}:=-\frac{1}{d}\text{Tr}\bm{\sigma}_{h}\approx p.

For each $K\in\mathcal{T}_{h},$ the postprocessed solution $\bm{u}_{h}^{*}|_{K}\in\mathcal{P}_{r+1}(K)$ is determined by

	$\displaystyle\nu(\nabla\bm{u}_{h}^{*},\nabla\bm{v}_{h})_{K}$	$\displaystyle=(\bm{\sigma}_{h},\nabla\bm{v}_{h})_{K}+(p_{h},\nabla\cdot\bm{v}_{h})_{K},~{}\bm{v}_{h}\in\mathcal{P}_{r+1}(K)/\mathcal{P}_{r}(K),$		(4.1)
	$\displaystyle(\bm{u}_{h}^{*},\bar{\bm{v}}_{h})_{K}$	$\displaystyle=(\bm{u}_{h},\bar{\bm{v}}_{h})_{K},\quad\bar{\bm{v}}_{h}\in\mathcal{P}_{r}(K).$		(4.1)

Since the analysis of the mapping $\bm{u}_{h}\mapsto\bm{u}_{h}^{*}$ is independent of equations, we can combine the analysis in Theorem 4.1 of CKP2011 for the Stokes equation with the supercloseness on $\|P_{h}\bm{u}-\bm{u}_{h}\|$ in Theorem 2.1 to obtain the following postprocessing superconvergence estimate.

Theorem 4.1

For sufficiently small $h$ , it holds that

\|\bm{u}-\bm{u}^{*}_{h}\|\lesssim h^{r+2}\big{(}\|\bm{u}\|_{r+2}+|\bm{\sigma}|_{r+1}+|\nabla\cdot\bm{\sigma}|_{r+1}\big{)}.

Proof

Let $\bar{\bm{u}}_{h}$ be the $L^{2}$ -projection of $\bm{u}$ onto the space

\left\{\bm{v}_{h}:\bm{v}_{h}|_{K}\in[\mathcal{P}_{r+1}(K)]^{d}\text{ for all }K\in\mathcal{T}_{h}\right\}.

It has been shown in the proof of Theorem 4.1 in CKP2011 that

		$\displaystyle h^{-1}\\|(\bm{I}-P_{h})(\bar{\bm{u}}_{h}-\bm{u}_{h}^{})\\|\lesssim h^{-1}\\|P_{h}(\bar{\bm{u}}_{h}-\bm{u}_{h}^{})\\|$		(4.2)
		$\displaystyle\quad+\\|\bm{\sigma}-\bm{\sigma}_{h}\\|+\left(\sum_{K\in\mathcal{T}_{h}}\\|\nabla(\bar{\bm{u}}_{h}-\bm{u})\\|^{2}_{L^{2}(K)}\right)^{\frac{1}{2}}.$		(4.2)

Using the triangle inequality, $P_{h}\bar{\bm{u}}_{h}=P_{h}\bm{u},$ $P_{h}{\bm{u}}^{*}_{h}=\bm{u}_{h}$ , and (4.2), we obtain

	$\displaystyle\\|\bm{u}-\bm{u}^{*}_{h}\\|$	$\displaystyle\leq\\|\bm{u}-\bar{\bm{u}}_{h}\\|+\\|(\bm{I}-P_{h})(\bar{\bm{u}}_{h}-\bm{u}^{}_{h})\\|+\\|P_{h}(\bar{\bm{u}}_{h}-\bm{u}^{}_{h})\\|$
		$\displaystyle\lesssim h^{r+2}\|\bm{u}\|_{r+2}+\\|P_{h}\bm{u}-\bm{u}_{h}\\|+h\\|\bm{\sigma}-\bm{\sigma}_{h}\\|.$

We finally conclude the proof from the previous estimate with (3.1) and (2.26b). ∎

Postprocessing technique on the scalar variable in mixed methods for Poisson’s equation can be found in e.g., BrezziDouglasMarini1985 ; BrambleXu1989 ; Stenberg1991 ; LovadinaStenberg2006 .

The postprocessing procedure $\bm{\sigma}_{h}\mapsto\bm{\sigma}_{h}^{*}$ can be derived from existing postprocessing operator $\widetilde{R}_{h}$ . In particular, $\widetilde{R}_{h}:\widetilde{\bm{\Sigma}}_{h}\rightarrow Y\subset[L^{2}(\Omega)]^{2}$ is a linear mapping onto the space $Y$ of suitable piecewise polynomials and satisfies


$\displaystyle\widetilde{R}_{h}\bm{c}$	$\displaystyle=\bm{c},\quad\bm{c}\in\mathbb{R}^{2},$	(4.3a)
$\displaystyle\\|\widetilde{R}_{h}\bm{v}\\|$	$\displaystyle\lesssim\\|\bm{v}\\|,\quad\bm{v}\in\widetilde{\bm{\Sigma}}_{h},$	(4.3b)
$\displaystyle\\|\bm{\tau}-\widetilde{R}_{h}\widetilde{\Pi}_{h}\bm{\tau}\\|$	$\displaystyle\lesssim h^{1+\alpha}\\|\bm{\tau}\\|_{W^{2}_{\infty}(\Omega)},$	(4.3c)

for some positive constant $\alpha$ and sufficiently smooth $\bm{\tau}.$ When $\widetilde{\bm{\Sigma}}_{h}$ is the lowest order RT element space, $\widetilde{R}_{h}$ that satisfies (4.3) is given in e.g., Brandts1994 ; BankLi2019 . The simple nodal or edge averaging ZZ1987 ; Brandts1994 and superconvergent patch recovery ZZ1992a ; XuZhang2004 ; BankLi2019 are also possible choices.

For $\bm{\tau}\in\bm{\Sigma}_{h}$ , let

R_{h}\bm{\tau}:=\widetilde{R}_{h}\bm{\tau}-\frac{1}{2|\Omega|}\left(\int_{\Omega}\text{Tr}\widetilde{R}_{h}\bm{\tau}dx\right)\bm{I},

where $\widetilde{R}_{h}$ is applied to each row of $\bm{\tau}$ . We have the following super-approximation result of ${R}_{h}.$

Lemma 4.1

Assume (4.3) holds. For $\bm{\tau}\in\bm{\Sigma}\cap[W^{2}_{\infty}(\Omega)]^{2\times 2}$ , we have

\|\bm{\tau}-{R}_{h}{\Pi}_{h}\bm{\tau}\|\lesssim h^{1+\alpha}\|\bm{\tau}\|_{W^{2}_{\infty}(\Omega)}.

Proof

(4.3a) implies $R_{h}\bm{I}=\bm{O}$ and thus ${R}_{h}{\Pi}_{h}\bm{\tau}={R}_{h}\widetilde{\Pi}_{h}\bm{\tau}.$ Then

		$\displaystyle\bm{\tau}-{R}_{h}{\Pi}_{h}\bm{\tau}=\bm{\tau}-{R}_{h}\widetilde{\Pi}_{h}\bm{\tau}$		(4.4)
		$\displaystyle\quad=\bm{\tau}-\widetilde{R}_{h}\widetilde{\Pi}_{h}\bm{\tau}-\frac{1}{2\|\Omega\|}\int_{\Omega}\text{Tr}(\bm{\tau}-\widetilde{R}_{h}\widetilde{\Pi}_{h}\bm{\tau})dx,$		(4.4)

where we used $\int_{\Omega}\text{Tr}\bm{\tau}dx=0$ . Lemma 4.1 then follows from (4.4), the Cauchy–Schwarz inequality, and (4.3c). ∎

Combining Lemma 4.1 and Theorem 3.1, we obtain the following superconvergence estimate by postprocessing.

Theorem 4.2

Let $\bm{\sigma}_{h}$ be the solution to (1.12) based on the lowest order RT element and $\bm{\sigma}^{*}_{h}:=R_{h}\bm{\sigma}_{h}$ . Let the assumptions in Theorem 3.1 and Lemma 4.1 hold. We have

\|\bm{\sigma}-\bm{\sigma}^{*}_{h}\|\lesssim\max(h^{2}|\log h|^{\frac{1}{2}},h^{1+\alpha})\big{(}\|\bm{\sigma}\|_{W^{2}_{\infty}(\Omega)}+\|\bm{u}\|_{2}\big{)}.

Proof

Using the triangle inequality, Theorem 3.1, and Lemma 4.1, we have

	$\displaystyle\\|\bm{\sigma}-\bm{\sigma}^{*}_{h}\\|$	$\displaystyle\leq\\|\bm{\sigma}-R_{h}\Pi_{h}\bm{\sigma}\\|+\\|R_{h}(\Pi_{h}\bm{\sigma}-\bm{\sigma}_{h})\\|$
		$\displaystyle\lesssim\max(h^{2}\|\log h\|^{\frac{1}{2}},h^{1+\alpha})\big{(}\\|\bm{\sigma}\\|_{W^{2}_{\infty}(\Omega)}+\\|\bm{u}\\|_{2}\big{)},$

which completes the proof. ∎

Remark 4.1

One could reconstruct accurate numerical symmetric stress and pressure from the superconvergent pseudostress $\bm{\sigma}_{h}^{*}$ . In fact, let

\displaystyle\bm{\sigma}^{S*}_{h}:=\frac{1}{2}(\bm{\sigma}_{h}^{*}+\bm{\sigma}_{h}^{*\top}),\quad p_{h}^{*}:=-\frac{1}{d}{\rm Tr}\bm{\sigma}_{h}^{*}.

We conclude that $\bm{\sigma}^{*S}_{h}$ superconverges to the symmetric stress $\bm{\sigma}^{S}$ and $p_{h}^{*}$ superconverges to $p$ the pressure from Theorem 4.2 and the elementary inequalities

	$\displaystyle\\|\bm{\sigma}^{S}-\bm{\sigma}^{*S}_{h}\\|$	$\displaystyle\leq\frac{1}{2}\\|\bm{\sigma}-\bm{\sigma}^{}_{h}\\|+\frac{1}{2}\\|\bm{\sigma}^{\top}-\bm{\sigma}^{\top}_{h}\\|=\\|\bm{\sigma}-\bm{\sigma}^{*}_{h}\\|,$
	$\displaystyle\\|p-p^{*}_{h}\\|$	$\displaystyle\leq\frac{1}{d}\\|{\rm Tr}(\bm{\sigma}-\bm{\sigma}^{}_{h})\\|\leq\frac{1}{\sqrt{d}}\\|\bm{\sigma}-\bm{\sigma}^{}_{h}\\|.$

5 Experiments

In this section, the postprocessing procedure $\bm{\sigma}_{h}\mapsto\bm{\sigma}_{h}^{*}$ in Theorem 4.2 is based on the polynomial preserving recovery $\widetilde{R}_{h}$ for the lowest order RT element described in BankLi2019 . Roughly speaking, each row of $\bm{\sigma}_{h}^{*}$ is constructed as a continuous piecewise linear vector-valued polynomial with nodal values determined by least-squares fitted linear local polynomial on vertex patches. For linear Oseen equations, we set the viscosity $\nu$ to be $1$ and verify the a priori and superconvergence error estimates. The adaptivity performance of the error estimator base on $\bm{\sigma}_{h}^{*}$ , $\bm{u}_{h}^{*}$ is also under investigation. In the end, the proposed postprocessing is tested in the Kovasznay flow from the Navier–Stokes equation.

Table 1: Convergence history of the lowest order

RT

element

nt	$\\|\bm{u}-\bm{u}_{h}\\|$	$\\|\bm{e}_{h}\\|$	$\\|\bm{u}-\bm{u}^{*}_{h}\\|$	$\\|\bm{\sigma}-\bm{\sigma}_{h}\\|$	$\\|\bm{\xi}_{h}\\|$	$\\|\bm{\sigma}-\bm{\sigma}_{h}^{*}\\|$
19	3.210e-1	3.851e-2	9.334e-2	1.372	2.505e-1	2.113
76	1.620e-1	1.027e-2	2.438e-2	6.849e-1	7.172e-2	6.108e-1
304	8.121e-2	2.621e-3	6.176e-3	3.418e-1	1.957e-2	1.724e-1
1216	4.063e-2	6.598e-4	1.550e-3	1.707e-1	5.242e-3	4.353e-2
4864	2.032e-2	1.653e-4	3.880e-4	8.530e-2	1.390e-3	1.088e-2
19456	1.016e-2	4.136e-5	9.703e-5	4.270e-2	3.659e-4	2.719e-3
order	9.990e-1	1.990	1.994	1.001	1.904	1.961

Table 2: Convergence history of the lowest order

BDM

element

nt	$\\|\bm{u}-\bm{u}_{h}\\|$	$\\|\bm{e}_{h}\\|$	$\\|\bm{u}-\bm{u}^{*}_{h}\\|$	$\\|\bm{\sigma}-\bm{\sigma}_{h}\\|$	$\\|\bm{\xi}_{h}\\|$
19	3.192e-1	1.779e-2	6.316e-02	4.480e-1	4.755e-1
76	1.618e-1	6.026e-3	1.645e-02	1.249e-1	1.180e-1
304	8.118e-2	1.636e-3	4.155e-03	3.216e-2	2.928e-2
1216	4.062e-2	4.181e-4	1.042e-03	8.104e-3	7.288e-3
4864	2.032e-2	1.051e-4	2.606e-04	2.031e-3	1.818e-3
19456	1.016e-2	2.631e-5	6.515e-05	5.081e-4	4.541e-4
order	9.986e-1	1.964	1.996	1.987	2.005

5.1 A priori convergence

Consider the Oseen equation (1.2) on the unit square $\Omega=[0,1]^{2}$ with the smooth solutions

	$\displaystyle\bm{u}$	$\displaystyle=\begin{pmatrix}\sin(\pi(x_{1}+x_{2}))\\ -\sin(\pi(x_{1}+x_{2}))\end{pmatrix},$
	$\displaystyle p$	$\displaystyle=x_{1}+x_{2}-1.$

We set $c=0$ , $\bm{b}=(\cos(x_{2}),\sin(x_{1}))^{\top}$ , $\bm{g}=\bm{u}|_{\partial\Omega}=\bm{0}$ . $\bm{f}$ is computed from $\bm{u}$ and $\bm{b}$ . We start with the initial partition in Figure 1 (left). A sequence of piecewise uniform meshes is obtained by uniform quad-refinement, i.e., dividing each triangle into four similar subtriangles by connecting the midpoints of each edge. Numerical results are presented in Tables 1 and 2, where nt is short for “number of triangles”. The order of convergence is computed from the error quantities in those tables by least squares without using the data in the first rows.

The numerical rates of convergence coincide with a priori error estimates (2.26), (2.27) and the supercloseness estimates in (3.1) and Theorem 3.1. It is noted that for the lowest order BDM element, $\|\bm{\xi}_{h}\|\approx\mathcal{O}(h^{2})$ , which is predicted by a priori error estimates and thus not supersmall. Since the recovery procedure in BankLi2019 provides the super-approximation rate $\alpha=1$ in (4.3c) and Lemma 4.1, the recovery superconvergence estimate for the lowest order RT element predicted by Theorem 4.2 is $\|\bm{\sigma}-\bm{\sigma}_{h}^{*}\|=\mathcal{O}(h^{2}|\log h|^{\frac{1}{2}})$ , numerically confirmed by the last column in Table 1.

5.2 Adaptive mesh refinement for a non-smooth problem

Consider (1.2) on the L-shaped domain $\Omega=[-1,1]^{2}\backslash([0,1]\times[-1,0])$ with the smooth pressure $p=x_{1}+x_{2}$ and singular velocity

\displaystyle\bm{u}

\displaystyle=\begin{pmatrix}r^{\alpha}\sin(\alpha\theta)\\ r^{\alpha}\cos(\alpha\theta)\end{pmatrix},

where $\alpha=\frac{2}{3}$ , $(r,\theta)$ is the polar coordinate w.r.t. $(0,0)$ . Set $c=0$ , $\bm{b}=(1,2)^{\top}$ and $\bm{g}=\bm{u}|_{\partial\Omega}$ . Direct calculation shows that $\nabla\cdot\bm{u}=0$ and $\bm{f}=(1,1)^{\top}+\bm{b}\cdot\nabla\bm{u}.$ In this experiment, we use the classical adaptive feedback loop (cf. BabuskaRheinboldt1978 ; Dorfler1996 ; MNS2000 )

\texttt{SOLVE}\rightarrow\texttt{ESTIMATE}\rightarrow\texttt{MARK}\rightarrow\texttt{REFINE}

to obtain a sequence of adaptively refined grids $\{\mathcal{T}_{h_{\ell}}\}_{\ell\geq 0}$ and numerical solutions $(\{\bm{\sigma}_{h_{\ell}},\bm{u}_{h_{\ell}})\}_{\ell\geq 0}$ . In particular, the algorithm starts from the initial grid $\mathcal{T}_{h_{0}}$ presented in Figure 2(left). The module ESTIMATE computes the superconvergent recovery-based error indicator

\mathcal{E}_{\ell}(K):=\big{(}\|\bm{\sigma}_{h_{\ell}}^{*}-\bm{\sigma}_{h_{\ell}}\|^{2}_{L^{2}(K)}+\|\bm{u}_{h_{\ell}}^{*}-\bm{u}_{h_{\ell}}\|^{2}_{L^{2}(K)}\big{)}^{\frac{1}{2}}

for each triangle $K\in\mathcal{T}_{h_{\ell}}$ . The module MARK then selects a collection of triangles

\mathcal{M}_{\ell}:=\{K\in\mathcal{T}_{\ell}:\mathcal{E}_{\ell}(K)\geq 0.7\max_{K^{\prime}\in\mathcal{T}_{\ell}}\mathcal{E}_{\ell}(K^{\prime})\}

(5.1)

to be refined by local quad-refinement. To remove the newly created hanging nodes, minimal number of neighboring elements of $\mathcal{M}_{\ell}$ are bisected and the next level triangulation $\mathcal{T}_{h_{\ell+1}}$ is generated. See Figure 2(right) for an adaptively refined triangulation.

It can be observed from Figure 3 that $(\bm{\sigma}_{h_{\ell}},\bm{u}_{h_{\ell}})$ optimally converges to $(\bm{\sigma},\bm{u})$ . In addtion, the errors $\|\bm{\sigma}-\bm{\sigma}^{*}_{h_{\ell}}\|$ and $\|\bm{u}-\bm{u}^{*}_{h_{\ell}}\|$ are apparently superconvergent to 0. Let

	$\displaystyle\mathcal{E}_{\ell}$	$\displaystyle:=\big{(}\sum_{K\in\mathcal{T}_{h_{\ell}}}\mathcal{E}_{\ell}(K)^{2}\big{)}^{\frac{1}{2}}=\big{(}\\|\bm{\sigma}^{}_{h_{\ell}}-\bm{\sigma}_{h_{\ell}}\\|^{2}+\\|\bm{u}^{}_{h_{\ell}}-\bm{u}_{h_{\ell}}\\|^{2}\big{)}^{\frac{1}{2}},$
	$\displaystyle E_{\ell}$	$\displaystyle:=\big{(}\\|\bm{\bm{\sigma}}-\bm{\sigma}_{h_{\ell}}\\|^{2}+\\|\bm{u}-\bm{u}_{h_{\ell}}\\|^{2}\big{)}^{\frac{1}{2}}.$

Due to the observed superconvergence phenomena and the triangle inequality

\displaystyle|\mathcal{E}_{\ell}-E_{\ell}|\leq\big{(}\|\bm{\bm{\sigma}}-\bm{\sigma}^{*}_{h_{\ell}}\|^{2}+\|\bm{u}-\bm{u}^{*}_{h_{\ell}}\|^{2}\big{)}^{\frac{1}{2}},

the recovery-based error estimator $\mathcal{E}_{\ell}$ is asymptotically exact, i.e.,

\displaystyle\lim_{\ell\rightarrow\infty}\frac{\mathcal{E}_{\ell}}{E_{\ell}}=1.

5.3 Adaptive mesh refinement for dominant convection

Consider the Oseen equation (1.2) on the unit square $\Omega=[0,1]^{2}$ with $c=0$ , $\bm{b}=(500,1)^{\top}$ , $\bm{g}=\bm{0}$ , $\bm{f}=5000(x_{2},-x_{1})^{\top}$ . In the end, we use the same adaptive algorithm in Subsection 5.2 to solve this convection-dominated example. The initial grid is given in Figure 4 (left). The marked set in (5.1) is replaced by

\mathcal{M}_{\ell}=\{K\in\mathcal{T}_{\ell}:\mathcal{E}_{\ell}(K)\geq 0.3\max_{K^{\prime}\in\mathcal{T}_{\ell}}\mathcal{E}_{\ell}(K^{\prime})\}.

Due to the convection coefficient $\bm{b}=(500,1)^{\top}$ , the exact solution $\bm{u}$ near the axis $x_{1}=1$ is rapidly changing to preserve the homogeneous Dirichlet boundary condition. It can be observed from Figures 4 and 5 that the adaptive mixed method is able to capture the boundary layer by adaptively graded grids.

5.4 Incompressible Navier–Stokes equation

Similarly to the linear Oseen equation, (1.3) could be recast into the following pseudostress-velocity form

	$\displaystyle(\mathcal{A}\bm{\sigma},\bm{\tau})+\nu(\nabla\cdot\bm{\tau},\bm{u})$	$\displaystyle=\nu\langle\bm{g},\bm{\tau}\bm{n}\rangle,\quad\bm{\tau}\in\bm{\Sigma},$		(5.2)
	$\displaystyle-\nu(\nabla\cdot\bm{\sigma},\bm{v})+(\bm{u}\cdot\mathcal{A}\bm{\sigma},\bm{v})$	$\displaystyle=\nu(\bm{f},\bm{v}),\quad\bm{v}\in\bm{V}.$		(5.2)

Let $\Omega=[-0.5,1.5]\times[0,2]$ . The exact solution solutions are taken to be

\bm{u}=\begin{pmatrix}1-e^{\lambda x_{1}}\cos(2\pi x_{2})\\ \frac{\lambda}{2\pi}e^{\lambda x_{1}}\sin(2\pi x_{2})\end{pmatrix},\quad p=-\frac{1}{2}e^{2\lambda x_{1}}+\frac{1}{8\lambda}(e^{3\lambda}-e^{-\lambda}),

(5.3)

where $\lambda=\frac{1}{2\nu}-\sqrt{\frac{1}{4\nu^{2}}+4\pi^{2}}$ , $\nu=0.025$ . In fact (5.3) is a well-known benchmark problem known as the Kovasznay flow (cf. DiPietroErn2012 ; ChenLiCorinaCimbala2020 ). In this experiment, we use the lowest order RT element $\bm{\Sigma}_{h}\times\bm{V}_{h}$ to discretize (5.2):

	$\displaystyle(\mathcal{A}\bm{\sigma}_{h},\bm{\tau})+\nu(\nabla\cdot\bm{\tau},\bm{u}_{h})$	$\displaystyle=\nu\langle\bm{g},\bm{\tau}\bm{n}\rangle,\quad\bm{\tau}\in\bm{\Sigma}_{h},$
	$\displaystyle-\nu(\nabla\cdot\bm{\sigma}_{h},\bm{v})+(\bm{u}\cdot\mathcal{A}\bm{\sigma}_{h},\bm{v})$	$\displaystyle=\nu(\bm{f},\bm{v}),\quad\bm{v}\in\bm{V}_{h}.$

The initial triangulation is a uniform partition of the square $\Omega$ with 512 right triangles. A sequence of grids is then generated by subdividing each triangles into four congruent subtriangles.

Although our analysis is devoted to the linear Oseen equation, one could observe apparent superconvergence in both pseudostress and velocity of the Navier–Stokes equation from Table 3.

Table 3: Convergence history of the lowest order

RT

element for (5.2)

nt	$\\|\bm{u}-\bm{u}_{h}\\|$	$\\|\bm{e}_{h}\\|$	$\\|\bm{u}-\bm{u}^{*}_{h}\\|$	$\\|\bm{\sigma}-\bm{\sigma}_{h}\\|$	$\\|\bm{\xi}_{h}\\|$	$\\|\bm{\sigma}-\bm{\sigma}_{h}^{*}\\|$
512	2.767e-1	1.717e-1	1.984e-1	1.572e-1	3.317e-2	1.311e-1
2048	1.225e-1	5.578e-2	6.133e-2	6.094e-2	8.866e-3	3.765e-2
8192	5.902e-2	2.233e-2	2.326e-2	2.716e-2	2.410e-3	1.063e-2
32768	2.920e-2	1.030e-2	1.043e-2	1.303e-2	6.509e-4	2.960e-3
order	1.079	1.350	1.415	1.194	1.889	1.823

Appendix

Proof (Proof of Lemma 3.2)

Let $\{\Omega_{i}\}_{i=1}^{N}$ be a polygonal partition of $\Omega$ and $\mathcal{T}_{h}|_{\Omega_{i}}$ be a uniform grid for each $i$ . Let $\bm{f}_{1}$ , $\bm{f}_{2}$ be fixed unit normals to two edges of an arbitrary but fixed triangle in $\mathcal{T}_{h}|_{\Omega_{i}}$ . Let $\bm{S}=\begin{pmatrix}\bm{f}_{1}^{\top}\\ \bm{f}_{2}^{\top}\end{pmatrix}$ and recall $\bm{e}_{1}=(1,0)$ , $\bm{e}_{2}=(0,1).$ We have

		$\displaystyle\int_{\Omega_{i}}[\bm{\bm{D}}(\bm{v}-\widetilde{\Pi}_{h}\bm{v})]\cdot\nabla^{\perp}w_{h}dx=\int_{\Omega_{i}}\bm{S}^{-\top}[\bm{\bm{D}}(\bm{v}-\widetilde{\Pi}_{h}\bm{v})]\cdot\bm{S}\nabla^{\perp}w_{h}dx$		(5.4)
		$\displaystyle\qquad=\sum_{j=1}^{2}\int_{\Omega_{i}}[\bm{e}_{j}\bm{S}^{-\top}\bm{\bm{D}}(\bm{v}-\widetilde{\Pi}_{h}\bm{v})](\bm{f}_{j}^{\top}\nabla^{\perp}w_{h})dx:=\sum_{j=1}^{2}I_{j}.$		(5.4)

For $j=1,2$ , let $\mathcal{E}^{o}_{j}$ and $\mathcal{E}^{\partial}_{j}$ be the set of interior (inside $\Omega_{i}$ ) and boundary (on $\partial\Omega_{i}$ ) edges orthogonal to $\bm{f}_{j}$ , respectively. Let $N_{E}$ denote the region which is the union of triangles sharing $E$ as an edge. To estimate $I_{1}$ , let $\Omega_{i}$ be partitioned into parallelograms $N_{E}$ with $E\in\mathcal{E}_{1}^{o}$ and boundary triangles $N_{E}$ with $E\in\mathcal{E}^{\partial}_{1}$ , see Figure 6, where $E_{1},E_{2}\in\mathcal{E}_{1}^{\partial}$ , $E_{3},E_{4}\in\mathcal{E}_{1}^{o}$ . For any $E\in\mathcal{E}_{1}^{o}$ , note that $\bm{f}_{1}^{\top}\nabla^{\perp}w_{h}$ is single-valued across $E$ and thus constant on $N_{E}.$ For $E\in\mathcal{E}_{1}^{\partial}$ , let $P_{E,1}$ and $P_{E,2}$ be the two endpoints of $E$ and $h_{E}$ denote the length of $E.$ Then

$\displaystyle I_{1}$	$\displaystyle=\sum_{E\in\mathcal{E}^{o}_{1}\cup\mathcal{E}^{\partial}_{1}}\int_{N_{E}}[\bm{e}_{1}\bm{S}^{-\top}\bm{\bm{D}}(\bm{v}-\widetilde{\Pi}_{h}\bm{v})](\bm{f}_{1}^{\top}\nabla^{\perp}w_{h})dx$	(5.5)
	$\displaystyle=\sum_{E\in\mathcal{E}^{o}_{1}}(\bm{f}_{1}^{\top}\nabla^{\perp}w_{h})\bm{e}_{1}\bm{S}^{-\top}\bm{\bm{D}}\int_{N_{E}}(\bm{v}-\widetilde{\Pi}_{h}\bm{v})dx$
	$\displaystyle+\sum_{E\in\mathcal{E}^{\partial}_{1}}h_{E}^{-1}\bm{e}_{1}\bm{S}^{-\top}\bm{\bm{D}}\int_{N_{E}}(\bm{v}-\widetilde{\Pi}_{h}\bm{v})dx\big{(}w_{h}(P_{E,2})-w_{h}(P_{E,1})\big{)}$
	$\displaystyle:=I_{11}+I_{12}.$

When $N_{E}$ is an exact parallelogram, it is shown in Equation (3.15) of Brandts1994 that

\left|\int_{N_{E}}(\bm{v}-\widetilde{\Pi}_{h}\bm{v})dx\right|\lesssim h^{3}|\bm{v}|_{H^{2}(N_{E})}.

(5.6)

Using (5.6) and the Cauchy-Schwarz inequality, we have

	$\displaystyle I_{11}$	$\displaystyle\lesssim\sum_{E\in\mathcal{E}^{o}_{1}}h^{-1}\\|\nabla^{\perp}w_{h}\\|_{L^{2}(N_{E})}\left\|\int_{N_{E}}(\bm{v}-\widetilde{\Pi}_{h}\bm{v})dx\right\|$		(5.7)
		$\displaystyle\lesssim\sum_{E\in\mathcal{E}^{o}_{1}}h^{2}\\|\nabla^{\perp}w_{h}\\|_{L^{2}(N_{E})}\|\bm{v}\|_{H^{2}(N_{E})}\lesssim h^{2}\\|\nabla^{\perp}w_{h}\\|\|\bm{v}\|_{2}.$		(5.7)

Let $\mathcal{N}$ denote the collection of endpoints of edges in $\mathcal{E}_{1}^{\partial}$ ,

\mathcal{N}_{1}:=\{P\in\mathcal{N}:\text{ $P$ is not a corner of $\partial\Omega_{i}$}\},

and $\mathcal{N}_{2}:=\mathcal{N}\backslash\mathcal{N}_{1}$ . For instance, $\mathcal{N}_{2}=\{P_{1},P_{2}\}$ in Figure 6. For $P\in\mathcal{N}_{1},$ let $E_{1}$ , $E_{2}$ be the two boundary edges sharing $P$ and $\omega_{P}$ be the union of three triangles having $P$ as a vertex. For $P\in\mathcal{N}_{2},$ let $E$ denote the unique edge in $\mathcal{E}_{1}^{\partial}$ having $P$ as a vertex. Rearranging the summation in $I_{12}$ , we have

	$\displaystyle I_{12}$	$\displaystyle=\sum_{P\in\mathcal{N}_{1}}h_{E_{1}}^{-1}\bm{e}_{1}\bm{S}^{-\top}\bm{\bm{D}}\left(\int_{N_{E_{1}}}(\bm{v}-\widetilde{\Pi}_{h}\bm{v})dx-\int_{N_{E_{2}}}(\bm{v}-\widetilde{\Pi}_{h}\bm{v})dx\right)w_{h}(P)$		(5.8)
		$\displaystyle+\sum_{P\in\mathcal{N}_{2}}h_{E}^{-1}\bm{e}_{1}\bm{S}^{-\top}\bm{\bm{D}}\int_{N_{E}}(\bm{v}-\widetilde{\Pi}_{h}\bm{v})dxw_{h}(P).$		(5.8)

For $P\in\mathcal{N}_{1}$ , Equation (3.15) in HuMaMa2021 shows that

\left|\int_{N_{E_{1}}}(\bm{v}-\widetilde{\Pi}_{h}\bm{v})dx-\int_{N_{E_{2}}}(\bm{v}-\widetilde{\Pi}_{h}\bm{v})dx\right|\lesssim h^{3}|\bm{v}|_{H^{2}(\omega_{P})}.

(5.9)

Without loss of generality, let $\int_{\Omega}w_{h}dx=0.$ Then the discrete Sobolev and Poincaré inequalities yield

\|w_{h}\|_{L^{\infty}(\Omega)}\lesssim|\log h|^{\frac{1}{2}}\|w_{h}\|_{1}\lesssim|\log h|^{\frac{1}{2}}\|\nabla^{\perp}w_{h}\|.

(5.10)

Standard interpolation error estimate yields

\left|\int_{N_{E}}(\bm{v}-\widetilde{\Pi}_{h}\bm{v})dx\right|\lesssim h^{2}|\bm{v}|_{H^{1}(N_{E})}.

(5.11)

It then follows from (5.8)–(5.11), and the Cauchy–Schwarz inequality that

$\displaystyle I_{12}$	$\displaystyle\lesssim\big{(}\sum_{P\in\mathcal{N}_{1}}h^{2}\|\bm{v}\|_{H^{2}(\omega_{P})}+\sum_{P\in\mathcal{N}_{2}}h\|\bm{v}\|_{H^{1}(N_{E})}\big{)}\\|w_{h}\\|_{L^{\infty}(\Omega)}$	(5.12)
	$\displaystyle\lesssim\big{(}h^{2}\|\bm{v}\|_{W^{2}_{\infty}(\Omega)}\sum_{P\in\mathcal{N}_{1}}\|\omega_{P}\|^{\frac{1}{2}}+h\|\bm{v}\|_{W^{1}_{\infty}(\Omega)}\sum_{P\in\mathcal{N}_{2}}\|N_{E}\|^{\frac{1}{2}}\big{)}\|\log h\|^{\frac{1}{2}}\\|\nabla^{\perp}w_{h}\\|$
	$\displaystyle\lesssim h^{2}\\|\bm{v}\\|_{W^{2}_{\infty}(\Omega)}\|\log h\|^{\frac{1}{2}}\\|\nabla^{\perp}w_{h}\\|.$

We finally conclude the proof from (5.4), (5.7), (5.12), the same analysis for $I_{2}$ and all pieces in the partition $\{\Omega_{i}\}_{i=1}^{N}$ . ∎

Declarations

Funding The authors did not receive support from any organization for this work.

Data Availability Data sharing is not applicable to this article as no datasets were generated or analysed during the current study.

Conflicts of interest The authors have no relevant financial or non-financial interests to disclose.

Code availability The code used in this study is available from the authors upon request.

References

(1) Antonietti, P.F., Beirão da Veiga, L., Mora, D., Verani, M.: A stream virtual element formulation of the Stokes problem on polygonal meshes. SIAM J. Numer. Anal. 52(1), 386–404 (2014). DOI 10.1137/13091141X
(2) Arnold, D.N., Douglas Jr., J., Gupta, C.P.: A family of higher order mixed finite element methods for plane elasticity. Numer. Math. 45(1), 1–22 (1984). DOI 10.1007/BF01379659
(3) Arnold, D.N., Falk, R.S., Winther, R.: Finite element exterior calculus, homological techniques, and applications. Acta Numer. 15, 1–155 (2006). DOI 10.1017/S0962492906210018
(4) Arnold, D.N., Falk, R.S., Winther, R.: Finite element exterior calculus: from Hodge theory to numerical stability. Bull. Amer. Math. Soc. (N.S.) 47(2), 281–354 (2010). DOI 10.1090/S0273-0979-10-01278-4
(5) Arnold, D.N., Li, L.: Finite element exterior calculus with lower-order terms. Math. Comp. 86(307), 2193–2212 (2017). DOI 10.1090/mcom/3158
(6) Babuška, I., Rheinboldt, W.C.: Error estimates for adaptive finite element computations. SIAM J. Numer. Anal. 15(4), 736–754 (1978). DOI 10.1137/0715049
(7) Bank, R.E., Li, Y.: Superconvergent recovery of Raviart-Thomas mixed finite elements on triangular grids. J. Sci. Comput. 81(3), 1882–1905 (2019). DOI 10.1007/s10915-019-01068-0
(8) Bank, R.E., Xu, J.: Asymptotically exact a posteriori error estimators. I. Grids with superconvergence. SIAM J. Numer. Anal. 41(6), 2294–2312 (2003). DOI 10.1137/S003614290139874X
(9) Barrios, T.P., Cascón, J.M., González, M.: Augmented mixed finite element method for the Oseen problem: a priori and a posteriori error analyses. Comput. Methods Appl. Mech. Engrg. 313, 216–238 (2017). DOI 10.1016/j.cma.2016.09.012
(10) Barrios, T.P., Cascón, J.M., González, M.: On an adaptive stabilized mixed finite element method for the Oseen problem with mixed boundary conditions. Comput. Methods Appl. Mech. Engrg. 365, 113,007, 21 (2020). DOI 10.1016/j.cma.2020.113007
(11) Boffi, D., Brezzi, F., Fortin, M.: Mixed finite element methods and applications, Springer Series in Computational Mathematics, vol. 44. Springer, Heidelberg (2013). DOI 10.1007/978-3-642-36519-5
(12) Bramble, J.H., Xu, J.: A local post-processing technique for improving the accuracy in mixed finite-element approximations. SIAM J. Numer. Anal. 26(6), 1267–1275 (1989). DOI 10.1137/0726073
(13) Brandts, J.H.: Superconvergence and a posteriori error estimation for triangular mixed finite elements. Numer. Math. 68(3), 311–324 (1994). DOI 10.1007/s002110050064
(14) Brandts, J.H.: Superconvergence for triangular order $k=1$ Raviart-Thomas mixed finite elements and for triangular standard quadratic finite element methods. Appl. Numer. Math. 34(1), 39–58 (2000). DOI 10.1016/S0168-9274(99)00034-3
(15) Brezzi, F., Douglas Jr., J., Durán, R., Fortin, M.: Mixed finite elements for second order elliptic problems in three variables. Numer. Math. 51(2), 237–250 (1987). DOI 10.1007/BF01396752
(16) Brezzi, F., Douglas Jr., J., Fortin, M., Marini, L.D.: Efficient rectangular mixed finite elements in two and three space variables. RAIRO Modél. Math. Anal. Numér. 21(4), 581–604 (1987). DOI 10.1051/m2an/1987210405811
(17) Brezzi, F., Douglas Jr., J., Marini, L.D.: Two families of mixed finite elements for second order elliptic problems. Numer. Math. 2(47), 217–235 (1985)
(18) Cáceres, E., Gatica, G.N.: A mixed virtual element method for the pseudostress-velocity formulation of the Stokes problem. IMA J. Numer. Anal. 37(1), 296–331 (2017). DOI 10.1093/imanum/drw002
(19) Cáceres, E., Gatica, G.N., Sequeira, F.A.: A mixed virtual element method for the Brinkman problem. Math. Models Methods Appl. Sci. 27(4), 707–743 (2017). DOI 10.1142/S0218202517500142
(20) Cai, Z., Tong, C., Vassilevski, P.S., Wang, C.: Mixed finite element methods for incompressible flow: stationary Stokes equations. Numer. Methods Partial Differential Equations 26(4), 957–978 (2010). DOI 10.1002/num.20467
(21) Cai, Z., Wang, C., Zhang, S.: Mixed finite element methods for incompressible flow: stationary Navier-Stokes equations. SIAM J. Numer. Anal. 48(1), 79–94 (2010). DOI 10.1137/080718413
(22) Cai, Z., Wang, Y.: A multigrid method for the pseudostress formulation of Stokes problems. SIAM J. Sci. Comput. 29(5), 2078–2095 (2007). DOI 10.1137/060661429
(23) Carrero, J., Cockburn, B., Schötzau, D.: Hybridized globally divergence-free LDG methods. I. The Stokes problem. Math. Comp. 75(254), 533–563 (2006). DOI 10.1090/S0025-5718-05-01804-1
(24) Carstensen, C., Gallistl, D., Schedensack, M.: Quasi-optimal adaptive pseudostress approximation of the Stokes equations. SIAM J. Numer. Anal. 51(3), 1715–1734 (2013)
(25) Carstensen, C., Kim, D., Park, E.J.: A priori and a posteriori pseudostress-velocity mixed finite element error analysis for the Stokes problem. SIAM J. Numer. Anal. 49(6), 2501–2523 (2011)
(26) Cesmelioglu, A., Cockburn, B., Nguyen, N.C., Peraire, J.: Analysis of HDG methods for Oseen equations. J. Sci. Comput. 55(2), 392–431 (2013). DOI 10.1007/s10915-012-9639-y
(27) Cesmelioglu, A., Cockburn, B., Qiu, W.: Analysis of a hybridizable discontinuous Galerkin method for the steady-state incompressible Navier-Stokes equations. Math. Comp. 86(306), 1643–1670 (2017). DOI 10.1090/mcom/3195
(28) Chen, X., Li, Y., Drapaca, C., Cimbala, J.: A unified framework of continuous and discontinuous galerkin methods for solving the incompressible navier-stokes equation. J. Comp. Phys. 422, 109,799 (2020). DOI 10.1016/j.jcp.2020.109799
(29) Christiansen, S.H., Winther, R.: Smoothed projections in finite element exterior calculus. Math. Comp. 77(262), 813–829 (2008)
(30) Cockburn, B., Sayas, F.J.: Divergence-conforming HDG methods for Stokes flows. Math. Comp. 83(288), 1571–1598 (2014). DOI 10.1090/S0025-5718-2014-02802-0
(31) Cui, M., Ye, X.: Superconvergence of finite volume methods for the Stokes equations. Numer. Methods Partial Differential Equations 25(5), 1212–1230 (2009). DOI 10.1002/num.20399
(32) Demlow, A.: Suboptimal and optimal convergence in mixed finite element methods. SIAM J. Numer. Anal. 39(6), 1938–1953 (2002). DOI 10.1137/S0036142900376900
(33) Di Pietro, D.A., Ern, A.: Mathematical aspects of discontinuous Galerkin methods, Mathématiques & Applications (Berlin) [Mathematics & Applications], vol. 69. Springer, Heidelberg (2012). DOI 10.1007/978-3-642-22980-0
(34) Dörfler, W.: A convergent adaptive algorithm for Poisson’s equation. SIAM J. Numer. Anal. 33(3), 1106–1124 (1996). DOI 10.1137/0733054
(35) Douglas Jr., J., Roberts, J.E.: Global estimates for mixed methods for second order elliptic equations. Math. Comp. 44(169), 39–52 (1985). DOI 10.2307/2007791
(36) Durán, R.: Superconvergence for rectangular mixed finite elements. Numer. Math. 58(3), 287–298 (1990). DOI 10.1007/BF01385626
(37) Eichel, H., Tobiska, L., Xie, H.: Supercloseness and superconvergence of stabilized low-order finite element discretizations of the Stokes problem. Math. Comp. 80(274), 697–722 (2011). DOI 10.1090/S0025-5718-2010-02404-4
(38) Ewing, R.E., Liu, M.M., Wang, J.: Superconvergence of mixed finite element approximations over quadrilaterals. SIAM J. Numer. Anal. 36(3), 772–787 (1999)
(39) Falk, R.S., Neilan, M.: Stokes complexes and the construction of stable finite elements with pointwise mass conservation. SIAM J. Numer. Anal. 51(2), 1308–1326 (2013). DOI 10.1137/120888132
(40) Fu, G., Jin, Y., Qiu, W.: Parameter-free superconvergent $H({\rm div})$ -conforming HDG methods for the Brinkman equations. IMA J. Numer. Anal. 39(2), 957–982 (2019). DOI 10.1093/imanum/dry001
(41) Gatica, G.N., Gatica, L.F., Márquez, A.: Analysis of a pseudostress-based mixed finite element method for the Brinkman model of porous media flow. Numer. Math. 126(4), 635–677 (2014). DOI 10.1007/s00211-013-0577-x
(42) Girault, V., Raviart, P.A.: Finite element methods for Navier-Stokes equations, Springer Series in Computational Mathematics, vol. 5. Springer-Verlag, Berlin (1986). DOI 10.1007/978-3-642-61623-5. Theory and algorithms
(43) Guzmán, J., Neilan, M.: Conforming and divergence-free Stokes elements on general triangular meshes. Math. Comp. 83(285), 15–36 (2014). DOI 10.1090/S0025-5718-2013-02753-6
(44) Hu, J., Ma, L., Ma, R.: Optimal superconvergence analysis for the Crouzeix-Raviart and the Morley elements. Advances in Computational Mathematics 47 (2021). DOI 10.1007/s10444-021-09874-7
(45) Hu, J., Yu, G.: A unified analysis of quasi-optimal convergence for adaptive mixed finite element methods. SIAM J. Numer. Anal. 56(1), 296–316 (2018)
(46) Huang, Y., Li, J., Wu, C.: Averaging for superconvergence: verification and application of 2d edge elements to maxwell’s equations in metamaterials. Comput. Methods Appl. Mech. Engrg. 255, 121–132 (2013)
(47) Kundu, P.K., Dowling, D.R., Tryggvason, G., Cohen, I.M.: Fluid mechanics, 6th edn. Academic Press (2015)
(48) Li, Y.: Quasi-optimal adaptive mixed finite element methods for controlling natural norm errors. Math. Comp. 90, 565–593 (2021). DOI 10.1090/mcom/3590
(49) Li, Y.: Superconvergent flux recovery of the Rannacher-Turek nonconforming element. J. Sci. Comput. 87(1), Paper No. 32, 19 (2021). DOI 10.1007/s10915-021-01445-8
(50) Li, Y.W.: Global superconvergence of the lowest-order mixed finite element on mildly structured meshes. SIAM J. Numer. Anal. 56(2), 792–815 (2018). DOI 10.1137/17M112587X
(51) Liu, H., Yan, N.: Superconvergence analysis of the nonconforming quadrilateral linear-constant scheme for Stokes equations. Adv. Comput. Math. 29(4), 375–392 (2008). DOI 10.1007/s10444-007-9054-3
(52) Liu, X., Li, J., Chen, Z.: A weak Galerkin finite element method for the Oseen equations. Adv. Comput. Math. 42(6), 1473–1490 (2016). DOI 10.1007/s10444-016-9471-2
(53) Lovadina, C., Stenberg, R.: Energy norm a posteriori error estimates for mixed finite element methods. Math. Comp. 75(256), 1659–1674 (2006). DOI 10.1090/S0025-5718-06-01872-2
(54) Morin, P., Nochetto, R.H., Siebert, K.G.: Data oscillation and convergence of adaptive FEM. SIAM J. Numer. Anal. 38(2), 466–488 (2000)
(55) Mu, L., Wang, J., Ye, X.: A stable numerical algorithm for the Brinkman equations by weak Galerkin finite element methods. J. Comput. Phys. 273, 327–342 (2014). DOI 10.1016/j.jcp.2014.04.017
(56) Pan, J.: Global superconvergence for the bilinear-constant scheme for the Stokes problem. SIAM J. Numer. Anal. 34(6), 2424–2430 (1997). DOI 10.1137/S0036142995286167
(57) Park, E.J., Seo, B.: An upstream pseudostress-velocity mixed formulation for the Oseen equations. Bull. Korean Math. Soc. 51(1), 267–285 (2014). DOI 10.4134/BKMS.2014.51.1.267
(58) Qian, Y., Wu, S., Wang, F.: A mixed discontinuous Galerkin method with symmetric stress for Brinkman problem based on the velocity-pseudostress formulation. arXiv e-prints arXiv:1907.01246 (2019)
(59) Raviart, P.A., Thomas, J.M.: A mixed finite element method for 2nd order elliptic problems. In: Mathematical aspects of finite element methods, pp. 292–315. Lecture Notes in Math., Vol. 606. (Proc. Conf., Consiglio Naz. delle Ricerche (C.N.R.), Rome (1977)
(60) Stenberg, R.: Postprocessing schemes for some mixed finite elements. RAIRO Modél. Math. Anal. Numér. 25(1), 151–167 (1991). DOI 10.1051/m2an/1991250101511
(61) Wang, J., Ye, X.: Superconvergence of finite element approximations for the Stokes problem by projection methods. SIAM J. Numer. Anal. 39(3), 1001–1013 (2001). DOI 10.1137/S003614290037589X
(62) Wang, J., Ye, X.: A weak Galerkin finite element method for the stokes equations. Adv. Comput. Math. 42(1), 155–174 (2016). DOI 10.1007/s10444-015-9415-2
(63) Xu, J., Zhang, Z.: Analysis of recovery type a posteriori error estimators for mildly structured grids. Math. Comp. 73(247), 1139–1152 (2004). DOI 10.1090/S0025-5718-03-01600-4
(64) Ye, X.: Superconvergence of nonconforming finite element method for the Stokes equations. Numer. Methods Partial Differential Equations 18(2), 143–154 (2002). DOI 10.1002/num.1036.abs
(65) Zhang, Z., Naga, A.: A new finite element gradient recovery method: superconvergence property. SIAM J. Sci. Comput. 26(4), 1192–1213 (2005). DOI 10.1137/S1064827503402837
(66) Zienkiewicz, O.C., Zhu, J.Z.: A simple error estimator and adaptive procedure for practical engineering analysis. Internat. J. Numer. Methods Engrg. 24(2), 337–357 (1987). DOI 10.1002/nme.1620240206
(67) Zienkiewicz, O.C., Zhu, J.Z.: The superconvergent patch recovery and a posteriori error estimates. I. The recovery technique. Internat. J. Numer. Methods Engrg. 33(7), 1331–1364 (1992). DOI 10.1002/nme.1620330702


$\displaystyle\\|\bm{v}-P_{h}\bm{v}\\|$	$\displaystyle\lesssim h^{s}\|\bm{v}\|_{s},$	(2.5a)
$\displaystyle\\|\bm{\tau}-\Pi_{h}\bm{\tau}\\|$	$\displaystyle\lesssim h^{s}\|\bm{\tau}\|_{s},$	(2.5b)
$\displaystyle\\|\nabla\cdot(\bm{\tau}-{\Pi}_{h}\bm{\tau})\\|$	$\displaystyle\lesssim h^{s}\|\operatorname{div}\bm{\tau}\|_{s},$	(2.5c)


$\displaystyle\\|\nabla\cdot\bm{\xi}_{h}\\|$	$\displaystyle\lesssim\\|\bm{\zeta}_{h}\\|+\\|\bm{\xi}_{h}\\|+\\|\bm{e}_{h}\\|+h\\|\bm{\rho}_{h}\\|,$	(2.17a)
$\displaystyle\\|\nabla\cdot\bm{\xi}_{h}\\|_{-1}$	$\displaystyle\lesssim h\big{(}{\left\|\kern-1.07639pt\left\|\kern-1.07639pt\left\|\bm{\zeta}_{h}\right\|\kern-1.07639pt\right\|\kern-1.07639pt\right\|}+{\left\|\kern-1.07639pt\left\|\kern-1.07639pt\left\|\bm{\xi}_{h}\right\|\kern-1.07639pt\right\|\kern-1.07639pt\right\|}+\\|\bm{\rho}_{h}\\|\big{)}+\\|\bm{e}_{h}\\|.$	(2.17b)


$\displaystyle\\|\bm{u}-\bm{u}_{h}\\|$	$\displaystyle\lesssim\\|\bm{\rho}_{h}\\|+h{\left\|\kern-1.07639pt\left\|\kern-1.07639pt\left\|\bm{\zeta}_{h}\right\|\kern-1.07639pt\right\|\kern-1.07639pt\right\|},$	(2.18a)
$\displaystyle{\left\|\kern-1.07639pt\left\|\kern-1.07639pt\left\|\bm{\sigma}-{\bm{\sigma}_{h}}\right\|\kern-1.07639pt\right\|\kern-1.07639pt\right\|}$	$\displaystyle\lesssim{\left\|\kern-1.07639pt\left\|\kern-1.07639pt\left\|\bm{\zeta}_{h}\right\|\kern-1.07639pt\right\|\kern-1.07639pt\right\|}+h\\|\bm{\rho}_{h}\\|,$	(2.18b)
$\displaystyle\\|\bm{\sigma}-{\bm{\sigma}_{h}}\\|$	$\displaystyle\lesssim\\|\bm{\zeta}_{h}\\|+h\\|\nabla\cdot\bm{\zeta}_{h}\\|+h\\|\bm{\rho}_{h}\\|.$	(2.18c)

	$\displaystyle\\|\bm{u}-\bm{u}_{h}\\|$	$\displaystyle\leq\\|\bm{\rho}_{h}\\|+\\|\bm{e}_{h}\\|,$
	$\displaystyle{\left\|\kern-1.07639pt\left\|\kern-1.07639pt\left\|\bm{\sigma}-\bm{\sigma}_{h}\right\|\kern-1.07639pt\right\|\kern-1.07639pt\right\|}$	$\displaystyle\leq{\left\|\kern-1.07639pt\left\|\kern-1.07639pt\left\|\bm{\zeta}_{h}\right\|\kern-1.07639pt\right\|\kern-1.07639pt\right\|}+{\left\|\kern-1.07639pt\left\|\kern-1.07639pt\left\|\bm{\xi}_{h}\right\|\kern-1.07639pt\right\|\kern-1.07639pt\right\|},$

	$\displaystyle\\|\bm{\xi}_{h}\\|^{2}$	$\displaystyle\lesssim(\mathcal{A}\bm{\xi}_{h},\bm{\xi}_{h})+\\|\nabla\cdot\bm{\xi}_{h}\\|_{-1}^{2}$
		$\displaystyle\lesssim\\|\bm{\zeta}_{h}\\|\\|\bm{\xi}_{h}\\|+\\|\nabla\cdot\bm{\xi}_{h}\\|\\|\bm{e}_{h}\\|+\\|\nabla\cdot\bm{\xi}_{h}\\|_{-1}^{2}$
		$\displaystyle\leq\frac{\varepsilon^{2}}{2}\\|\bm{\xi}_{h}\\|^{2}+\frac{\varepsilon^{-2}}{2}\\|\bm{\zeta}_{h}\\|^{2}+\frac{\varepsilon^{2}}{2}\\|\nabla\cdot\bm{\xi}_{h}\\|^{2}+\frac{\varepsilon^{-2}}{2}\\|\bm{e}_{h}\\|^{2}+\\|\nabla\cdot\bm{\xi}_{h}\\|_{-1}^{2},$