Analysis of recovery type a posteriori error estimators for mildly structured grids

Abstract

Some recovery type error estimators for linear finite elements are analyzed under $O(h^{1+\alpha })$ $(\alpha > 0)$ regular grids. Superconvergence of order $O(h^{1+\rho })$ $(0 < \rho \le \alpha )$ is established for recovered gradients by three different methods. As a consequence, a posteriori error estimators based on those recovery methods are asymptotically exact.

1. Introduction

A posteriori error estimates have become standard in modern engineering and scientific computation. There are two types of popular error estimators: the residual type (see, e.g., 2, 4) and the recovery type (see, e.g., 21). The most representative recovery type error estimator is the Zienkiewicz-Zhu error estimator, especially the estimator based on gradient patch recovery by local discrete least-squares fitting 22, 23. The method is now widely used in engineering practice for its robustness in a posteriori error estimates and its efficiency in computer implementation. It is a common belief that the robustness of the ZZ estimator is rooted in the superconvergence property of the associated gradient recovery under structured meshes. Superconvergence properties of the ZZ recovery based on local least-squares fitting are proven by Zhang 17 for all popular elements under rectangular meshes, by Li-Zhang 11 for linear elements under strongly regular triangular meshes, and by Zhang-Victory 18 for tensor product elements under strongly regular quadrilateral meshes.

While there is a sizable literature on theoretical investments for residual type error estimators (see, e.g., 1, 3, 10, 14 and references therein), there have not been many theoretical results on recovery type error estimators. Nevertheless, the recovery type error estimators perform astonishingly well even for unstructured grids. The current paper intends to explain this phenomenon. We observe that for an unstructured mesh, when adaptive procedure is used, a mesh refinement will usually bring in some kind of local structure. It is then reasonable to assume that for most of the domain, every two adjacent triangles form an $O(h^{1+\alpha })$ approximate parallelogram. Under this assumption, we are able to establish superconvergence of the gradient recovery operator for three popular methods: weighted averaging, local $L^2$-projection, and the ZZ patch recovery. Furthermore, by utilizing an integral identity for linear elements on one triangular element developed by Bank and Xu 5, we are able to generalize their superconvergence result between the finite element solution and the linear interpolation from an $O(h^2)$ regular grid to an $O(h^{1+\alpha })$ regular grid. Finally, we are able to prove asymptotic exactness of the three recovery error estimators.

The topic of a posteriori error estimates has recently attracted more and more attention in the scientific community (see, e.g., 5, 6, 7, 9, 16, 20; also see recent books 1, 3 for some general discussions). The literature regarding finite element superconvergence theory can be found in the books 8, 10, 12, 15, 19.

2. Geometry identities of a triangle

In this section, we shall generalize the result in 5 for $\alpha = 1$ to all $\alpha > 0$. Following the argument in 5, we consider in Figure 1, a triangle $\tau$ with vertices $\mathbf{p}_k^{t}=(x_k,y_k)$, $1\leq k \leq 3$, oriented counterclockwise, and corresponding nodal basis functions (barycentric coordinates) $\{ \phi _k \}_{k=1}^3$. Let $\{ e_k \}_{k=1}^3$ denote the edges of element $\tau$, $\{ \theta _k \}_{k=1}^3$ the angles, $\{ \mathbf{n}_k \}_{k=1}^3$ the unit outward normal vectors, $\{ \mathbf{t}_k \}_{k=1}^3$ the unit tangent vectors with counterclockwise orientation, $\{ \ell _k \}_{k=1}^3$ the edge lengths, and $\{ d_k \}_{k=1}^3$ the perpendicular heights. Let $\tilde{\mathbf{p}}$ be the point of intersection for the perpendicular bisectors of the three sides of $\tau$. Let $|s_k|$ denote the distance between $\tilde{\mathbf{p}}$ and side $k$. If $\tau$ has no obtuse angles, then the $s_k$ will be nonnegative. Otherwise, the distance to the side opposite the obtuse angle will be negative.

Let ${\mathcal{D}}_{\tau }$ be a symmetric $2\times 2$ matrix with constant entries. We define

$$\begin{equation*} \xi _k=-\mathbf{n}_{k+1}\cdot {\mathcal{D}}_{\tau }\mathbf{n}_{k-1}. \end{equation*}$$

The important special case ${\mathcal{D}}_{\tau }=I$ corresponds to $-\Delta$, and in this case $\xi _k=\cos \theta _k$. Let $q_k=\phi _{k+1}\phi _{k-1}$ denote the quadratic bump function associated with edge $e_k$ and let $\psi _k=\phi _k(1-\phi _k)$.

The following fundamental identity is proved in 5 for $v_h \in P_1(\tau )$:

$$\begin{equation} \begin{split} \int _{\tau }&\nabla (u-u_I)\cdot {\mathcal{D}}_{\tau } \nabla v_h = \sum _{k=1}^3 \int _{e_k} \frac{\xi _k q_k}{2\sin \theta _k} \left\{ (\ell _{k+1}^2-\ell _{k-1}^2) \frac{\partial ^2u}{\partial \mathbf{t}_k^2}+4|\tau | \frac{\partial ^2u}{\partial \mathbf{t}_k\partial \mathbf{n}_k} \right\}\frac{\partial v_h}{\partial \mathbf{t}_k}\\ &-\int _{\tau } \sum _{k=1}^3 \frac{\ell _k\xi _k}{2\sin ^2\theta _k} \left\{ \ell _{k+1}\psi _{k-1} \frac{\partial ^3 u}{\partial ^2\mathbf{t}_{k+1}\partial \mathbf{t}_{k-1}} + \ell _{k-1}\psi _{k+1} \frac{\partial ^3 u}{\partial ^2\mathbf{t}_{k-1}\partial \mathbf{t}_{k+1}} \right\}\frac{\partial v_h}{\partial \mathbf{t}_k}, \end{split} {\tag{2.1}} \end{equation}$$

where $u_I\in P_1(\tau )$ is the linear interpolation of $u$ on $\tau$.

We say that two adjacent triangles (sharing a common edge) form an $O(h^{1+\alpha })$ ($\alpha > 0$) approximate parallelogram if the lengths of any two opposite edges differ only by $O(h^{1+\alpha })$.

Definition.

The triangulation ${\mathcal{T}}_h={\mathcal{T}}_{1,h}\cup {\mathcal{T}}_{2,h}$ is said to satisfy Condition $(\alpha ,\sigma )$ if there exist positive constants $\alpha$ and $\sigma$ such that every two adjacent triangles inside ${\mathcal{T}}_{1,h}$ form an $O(h^{1+\alpha })$ parallelogram and

$$\begin{equation*} \bar{\Omega }_{1,h} \cup \bar{\Omega }_{2,h} = \bar{\Omega }, \quad |\Omega _{2,h}|=O(h^\sigma ),\quad \bar{\Omega }_{i,h}\equiv \bigcup _{\tau \in {\mathcal{T}}_{i,h}}\bar{\tau }, \quad i=1,2. \end{equation*}$$

Remark.

There are two important ingredients in an automatic mesh generation code. One, called swap diagonal, changes the direction of some diagonal edges in order to obtain near parallel directions for adjacent element edges and to make as many nodes as possible have six triangles attached. Another, known as Lagrange smoothing, iteratively relocates nodes to place each node near a mesh symmetry center (see condition (3.1) in Section 3).

Clearly, both swap diagonal and Lagrange smoothing are intended to make every two adjacent triangles form an $O(h^{1+\alpha })$ parallelogram. Eventually, only a small portion of elements (including boundary elements) do not satisfy this condition. These elements then belong to $\Omega _{2,h}$, which has a small measure. Therefore, Condition $(\alpha ,\sigma )$ is a reasonable condition in practice and can be satisfied by most meshes produced by automatic mesh generation codes.

Denote ${\mathcal{V}}_h \subset H^1(\Omega )$, the $C^0$ linear finite element space associated with ${\mathcal{T}}_h$.

Lemma 2.1.

Assume that ${\mathcal{T}}_h$ satisfy Condition $(\alpha ,\sigma )$. Let ${\mathcal{D}}_{\tau }$ be a piecewise constant matrix function defined on ${\mathcal{T}}_h$, whose elements ${\mathcal{D}}_{\tau ij}$ satisfy

$$\begin{gather*} |{\mathcal{D}}_{\tau ij}|\lesssim 1, \quad |{\mathcal{D}}_{\tau ij}-{\mathcal{D}}_{\tau ' ij} |\lesssim h^\alpha , \quad i=1,2; \; j=1,2. \end{gather*}$$

Here $\tau$ and $\tau '$ are a pair of triangles sharing a common edge. Then for any $v_h \in {\mathcal{V}}_h$

$$\begin{equation} \left| \sum _{\tau \in {\mathcal{T}}_h}\int _\tau \nabla (u-u_I)\cdot {\mathcal{D}}_{\tau }\nabla v_h \right| \lesssim h^{1+\rho } ( \|u\|_{3,\Omega } + | u |_{2,\infty ,\Omega } ) | v |_{1,\Omega }, \quad \rho = \min (\alpha ,\frac{\sigma }{2},\frac{1}{2}), {\tag{2.2}} \end{equation}$$

where $u_I \in {\mathcal{V}}_h$ is the interpolation of $u$.

Proof.

Applying (2.1),

$$\begin{equation} \sum _{\tau \in {\mathcal{T}}_h}\int _\tau \nabla (u-u_I)\cdot {\mathcal{D}}_{\tau }\nabla v_h =I_1+I_2 {\tag{2.3}} \end{equation}$$

where

$$\begin{eqnarray*} I_1&=&\sum _{\tau \in {\mathcal{T}}_h}\sum _{k=1}^3 \int _{e_k}\frac{\xi _k q_k}{2\sin \theta _k} \left\{ (\ell _{k+1}^2-\ell _{k-1}^2) \frac{\partial ^2u}{\partial \mathbf{t}_k^2}+4|\tau | \frac{\partial ^2u}{\partial \mathbf{t}_k\partial \mathbf{n}_k} \right\}\frac{\partial v_h}{\partial \mathbf{t}_k}, \\ I_2&=&-\sum _{\tau \in {\mathcal{T}}_h} \int _{\tau } \sum _{k=1}^3 \frac{\ell _k\xi _k}{2\sin ^2\theta _k}\\ &&\qquad \times \left\{ \ell _{k+1}\psi _{k-1} \frac{\partial ^3 u}{\partial ^2\mathbf{t}_{k+1}\partial \mathbf{t}_{k-1}} + \ell _{k-1}\psi _{k+1} \frac{\partial ^3 u}{\partial ^2\mathbf{t}_{k-1}\partial \mathbf{t}_{k+1}} \right\}\frac{\partial v_h}{\partial \mathbf{t}_k}. \\ \end{eqnarray*}$$

$I_2$ is easily estimated by

$$\begin{equation} |I_2|\lesssim h^2|\!| u |\!|_{3,\Omega }| v_h |_{1,\Omega }. {\tag{2.4}} \end{equation}$$

To estimate $I_1$, we separate all interior edges into two different groups. ${\mathcal{E}}_1$ is the set of edges $e$ such that the two adjacent triangles sharing $e$ form an $O(h^{1+\alpha })$ approximate parallelogram and ${\mathcal{E}}_2$ is the set of the remaining interior edges. The set of all interior edges is given by ${\mathcal{E}}={\mathcal{E}}_1+{\mathcal{E}}_2$.

For each $e\in {\mathcal{E}}$, there are two triangles, say $\tau$ and $\tau '$, that share $e$ as a common edge. Denote, with respect to $\tau$,

$$\begin{equation*} \alpha _e=\frac{\xi _k}{2\sin \theta _k} (\ell _{k+1}^2-\ell _{k-1}^2),\hspace{0.3cm} \beta _e=\frac{\xi _k}{2\sin \theta _k}4|\tau |, \end{equation*}$$

and with respect to $\tau '$,

$$\begin{equation*} \alpha '_e=\frac{\xi _{k'}}{2\sin \theta _{k'}} (\ell _{k'+1}^2-\ell _{k'-1}^2),\hspace{0.3cm} \beta '_e=\frac{\xi _{k'}}{2\sin \theta _{k'}}4|\tau '|. \end{equation*}$$

Taking $\mathbf{n}$ and $\mathbf{t}$ to correspond to $\tau$, we can write

$$\begin{equation*} I_1=I_{11}+I_{12}+I_{13}, \end{equation*}$$

where

$$\begin{equation*} I_{1j}=\sum _{e\in {\mathcal{E}}_j} \int _e q_e \left\{ (\alpha _e-\alpha '_e)\frac{\partial ^2u}{\partial \mathbf{t}^2} +(\beta _e-\beta '_e)\frac{\partial ^2u}{\partial \mathbf{t}\partial \mathbf{n}} \right\}\frac{\partial v_h}{\partial \mathbf{t}} \end{equation*}$$

for $j=1,2$, and

$$\begin{equation*} I_{13}=\sum _{e\subset \partial \Omega }\int _e q_e \left\{ \alpha _e\frac{\partial ^2u}{\partial \mathbf{t}^2} +\beta _e\frac{\partial ^2u}{\partial \mathbf{t}\partial \mathbf{n}} \right\}\frac{\partial v_h}{\partial \mathbf{t}}. \end{equation*}$$

It is easy to see that, if $v_h=0$ on $\partial \Omega$, then $I_{13}=0$. Otherwise, we have the following estimate:

$$\begin{equation} |I_{13}|\lesssim h^{3/2}| u |_{2,\infty ,\partial \Omega } | v_h |_{1,\Omega }. {\tag{2.5}} \end{equation}$$

Setting $\mathbf{z}=\mathbf{t}$ and $\mathbf{z}=\mathbf{n}$, we estimate

$$\begin{equation} \left|\int _e q_e \frac{\partial ^2u}{\partial \mathbf{t}\partial \mathbf{z}} \frac{\partial v_h}{\partial \mathbf{t}}\right| \lesssim h^{-1}| u |_{2,\infty ,\Omega }\int _\tau |\nabla v_h|. {\tag{2.6}} \end{equation}$$

By definition, for $e\in {\mathcal{E}}_1$, $\alpha '_e=\alpha _e(1+O(h^\alpha ))$ and $\beta '_e=\beta _e(1+O(h^\alpha ))$. Therefore

$$\begin{gather*} |\alpha _e-\alpha '_e|\lesssim h^{2+\alpha }, \quad |\beta _e-\beta '_e|\lesssim h^{2+\alpha }. \end{gather*}$$

Combining this with 2.6, we have

$$\begin{equation} |I_{11}| \lesssim h^{1+\alpha } |u|_{2,\infty ,\Omega } \int _\Omega |\nabla v_h| \lesssim h^{1+\alpha } | u |_{2,\infty ,\Omega }| v_h |_{1,\Omega }. {\tag{2.7}} \end{equation}$$

Now we turn to the estimate for $I_{12}$. Since adjacent elements in $\Omega _{2,h}$ do not form an $O(h^{1+\alpha })$ approximate parallelogram, we simply estimate

$$\begin{gather*} |\alpha _e-\alpha '_e|\leq |\alpha _e|+|\alpha '_e|\lesssim h^2, \quad |\beta _e-\beta '_e|\leq |\beta _e|+|\beta '_e|\lesssim h^2. \end{gather*}$$

Similarly to 2.7, this leads to

$$\begin{equation*} |I_{12}|\lesssim h |u|_{2,\infty ,\Omega } \sum _{\tau \in {\mathcal{T}}_{2,h}} \int _\tau |\nabla v_h| \lesssim h |u|_{2,\infty ,\Omega } \|\nabla v_h\|_{0,\Omega _{2,h}} h^{\sigma /2}. \end{equation*}$$

Combining this with 2.5 and 2.7 leads to

$$\begin{equation} |I_1|\lesssim h^{1+\rho } | u |_{2,\infty ,\Omega } | v_h |_{1,\Omega }. {\tag{2.8}} \end{equation}$$

Finally, applying 2.4 and 2.8 to 2.3, we obtain 2.2.

■

3. Gradient recovery operators

We define ${\mathcal{N}}_h$ as the nodal set of a quasi-uniform triangulation ${\mathcal{T}}_h$. Given $z\in {\mathcal{N}}_h$, we consider an element patch $\omega$ around $z$, which we choose as the origin of a local coordinates. Let $(x_j,y_j)$ be the barycenter of a triangle $\tau _j\subset \omega$, $j=1,2,\ldots ,m$. We require that one of the following two geometric conditions be satisfied for $\alpha \ge 0$:

$$\begin{equation} \frac{1}{m} \sum _{j=1}^m (x_j,y_j) = O(h^{1+\alpha }) (1,1). {\tag{3.1}} \end{equation}$$

$$\begin{equation} \sum _{j=1}^m \frac{|\tau _j|}{|\omega |} (x_j,y_j) = O(h^{1+\alpha }) (1,1). {\tag{3.2}} \end{equation}$$ Here we use $(x_j,y_j)$ to represent a vector in conditions (3.1) and (3.2).

Remark.

Condition $(\alpha ,\sigma )$ implies both conditions (3.1) and (3.2) for $z \in {\mathcal{N}}_h \cap \Omega _{1,h}$. Indeed, conditions (3.1) and (3.2) are trivially (with $\alpha = \infty$) satisfied by uniform meshes of the regular pattern, the Union Jack pattern, and the criss-cross pattern, and allow an $O(h^{1+\alpha })$ deviation from those meshes. For example, a strongly regular mesh is an $O(h^2)$ deviation from a uniform mesh of the regular pattern. Note that the condition (3.1) depends only on relative positions of the barycenters of the triangles and is independent of the shapes, sizes, and numbers of those triangles.

A boundary node $z$ usually leads to $\alpha = 0$. However, if $z$ is an interior node with $\alpha = 0$, then there are no restrictions and we have a completely unstructured mesh around $z$.

Let $u_I\in {\mathcal{V}}_h$ be the linear interpolation of a given function $u$. We shall discuss a gradient recovery operator $G_h$ and prove the superconvergent property between $\nabla u$ and $G_hu_I$.

The value of $G_hu_I$ is first determined at a vertex and then linearly interpolated over the whole domain. There are three popular ways to generate $G_hu_I$ at a vertex $z$.

(a) Weighted averaging.

$$\begin{equation} G_hu_I(z) = \sum _{j=1}^m \frac{|\tau _j|}{|\omega |} \nabla u_I(x_j,y_j). {\tag{3.3}} \end{equation}$$

(b) Local $L^2$-projection. We seek linear functions $p_l\in P_1(\omega )$ ($l=1,2$), such that

$$\begin{equation} \int _\omega [ p_l(x,y) - \partial _lu_I(x,y) ]q(x,y) dxdy = 0, \quad \forall q \in P_1(\omega ), \quad l=1,2. {\tag{3.4}} \end{equation}$$

Then we define $G_hu_I(z) = (p_1(0,0),p_2(0,0))$.

(c) Local discrete least-squares fitting proposed by Zienkiewicz-Zhu 22. We seek linear functions $p_l\in P_1(\omega )$ ($l=1,2$), such that

$$\begin{equation} \sum _{j=1}^m [ p_l(x_j,y_j) - \partial _lu_I(x_j,y_j) ]q(x_j,y_j) = 0, \quad \forall q \in P_1(\omega ), \quad l=1,2. {\tag{3.5}} \end{equation}$$

Then we define $G_hu_I(z) = (p_1(0,0),p_2(0,0))$.

Note that (c) is a discrete version of (b). The existence and uniqueness of the minimizers in (b) and (c) can be found in 11, Lemma 1. The following theorem generalizes the result in 11 from $\alpha = 1$ to $\alpha > 0$.

Theorem 3.1.

Let $\omega$ be an element patch around a node $z\in {\mathcal{N}}_h$, let $u \in W^3_\infty (\omega )$, and let $G_hu_I(z)$ be produced by either the local $L^2$-projection or the weighted averaging under condition 3.2, or by the local discrete least-squares fitting under condition 3.1. Then

$$\begin{equation*} | G_hu_I(z) - \nabla u(z) | \lesssim h^{1+\alpha } \|u\|_{3,\infty ,\omega }. \end{equation*}$$

Proof.

(a) For the weighted averaging, we have

$$\begin{eqnarray*} && \sum _{j=1}^m \frac{|\tau _j|}{|\omega |} \partial _lu_I(x_j,y_j) - \partial _lu(0,0) \\ &=& \sum _{j=1}^m \frac{|\tau _j|}{|\omega |} \partial _l (u_I-u)(x_j,y_j) + \sum _{j=1}^m \frac{|\tau _j|}{|\omega |} [\partial _lu(x_j,y_j) - \partial _lu(0,0)] \\ &=& \sum _{j=1}^m \frac{|\tau _j|}{|\omega |} \partial _l (u_I-u)(x_j,y_j) + \nabla \partial _lu(0,0) \cdot \sum _{j=1}^m \frac{|\tau _j|}{|\omega |} (x_j,y_j) + R_1(u), \end{eqnarray*}$$

where, by the Taylor expansion,

$$\begin{equation*} |R_1(u)| \lesssim h^2|u|_{3,\infty ,\omega }. \end{equation*}$$

Since the barycenter is the derivative superconvergent point for the linear interpolation, then

$$\begin{equation*} | \partial _l (u_I-u)(x_j,y_j)| \lesssim h^2|u|_{3,\infty ,\omega }, \quad j=1,2,\ldots ,m. \end{equation*}$$

Recall the condition (3.2), and we derive

$$\begin{equation*} | \nabla \partial _lu(0,0) \cdot \sum _{j=1}^m \frac{|\tau _j|}{|\omega |}(x_j,y_j) | \lesssim h^{1+\alpha }|u|_{2,\infty ,\omega }. \end{equation*}$$

Therefore,

$$\begin{equation} | \sum _{j=1}^m \frac{|\tau _j|}{|\omega |} \partial _lu_I(x_j,y_j) - \partial _lu(0,0) | \lesssim h^{1+\alpha } \|u\|_{3,\infty ,\omega }. {\tag{3.6}} \end{equation}$$

(b) For the local $L^2$-projection, we set $q=1$ in (3.4) to obtain

$$\begin{equation*} \sum _{j=1}^m |\tau _j| p_l(x_j,y_j) = \sum _{j=1}^m |\tau _j| \partial _lu_I(x_j,y_j). \end{equation*}$$

Therefore,

$$\begin{align*} p_l(0,0) - \sum _{j=1}^m \frac{|\tau _j|}{|\omega |} \partial _lu_I(x_j,y_j) &= p_l(0,0) - \sum _{j=1}^m \frac{|\tau _j|}{|\omega |} p_l(x_j,y_j)\\ &= - \nabla p_l(0,0) \cdot \sum _{j=1}^m \frac{|\tau _j|}{|\omega |} (x_j,y_j). \end{align*}$$

Using (see 11, Lemma 2)

$$\begin{equation} |\nabla p_l(0,0)| \lesssim \|u\|_{3,\infty ,\omega } {\tag{3.7}} \end{equation}$$

and condition (3.2), we obtain

$$\begin{equation} | p_l(0,0) - \sum _{j=1}^m \frac{|\tau _j|}{|\omega |} \partial _lu_I(x_j,y_j) | \lesssim h^{1+\alpha } \|u\|_{3,\infty ,\omega }. {\tag{3.8}} \end{equation}$$

Combining (3.6) and (3.8), we have proved

$$\begin{equation} |p_l(0,0) - \partial _lu(0,0)| \lesssim h^{1+\alpha } \|u\|_{3,\infty ,\omega }. {\tag{3.9}} \end{equation}$$

(c) For the local discrete least-squares fitting, we set $q=1$ in (3.5) to obtain

$$\begin{equation*} \sum _{j=1}^m p_l(x_j,y_j) = \sum _{j=1}^m \partial _lu_I(x_j,y_j). \end{equation*}$$

Therefore,

$$\begin{align*} p_l(0,0) - \frac{1}{m} \sum _{j=1}^m \partial _lu_I(x_j,y_j) &= p_l(0,0) - \frac{1}{m} \sum _{j=1}^m p_l(x_j,y_j)\\ &= - \frac{1}{m} \nabla p_l(0,0) \cdot \sum _{j=1}^m (x_j,y_j). \end{align*}$$

Using (3.7) and condition (3.1), we obtain

$$\begin{equation} | p_l(0,0) - \frac{1}{m} \sum _{j=1}^m \partial _lu_I(x_j,y_j) | \lesssim h^{1+\alpha } \|u\|_{3,\infty ,\omega }. {\tag{3.10}} \end{equation}$$

Next,

$$\begin{eqnarray*} && \frac{1}{m} \sum _{j=1}^m \partial _lu_I(x_j,y_j) - \partial _lu(0,0) \\ &=& \frac{1}{m} \sum _{j=1}^m \partial _l (u_I-u)(x_j,y_j) + \frac{1}{m} \sum _{j=1}^m [\partial _lu(x_j,y_j) - \partial _lu(0,0)] \\ &=& \frac{1}{m} \sum _{j=1}^m \partial _l (u_I-u)(x_j,y_j) + \frac{1}{m} \nabla \partial _lu(0,0) \cdot \sum _{j=1}^m (x_j,y_j) + R_2(u), \end{eqnarray*}$$

with $|R_2(u)| \lesssim h^2|u|_{3,\infty ,\omega }$. Therefore,

$$\begin{equation} | \frac{1}{m} \sum _{j=1}^m \partial _lu_I(x_j,y_j) - \partial _l u(0,0) | \lesssim h^{1+\alpha }\|u\|_{3,\infty ,\omega }. {\tag{3.11}} \end{equation}$$

Combining (3.10) and (3.11), we obtain (3.9) for the current case.

■

Theorem 3.2.

The recovery operator $G_h$ satisfies

$$\begin{equation*} G_h v(z) = \sum _{j=1}^m c_j \nabla v(x_j,y_j), \quad \sum _{j=1}^m c_j = 1, \end{equation*}$$

in all three cases unconditionally. Furthermore, $c_j > 0$ for

(a) the weighted averaging unconditionally;

(b) the local $L^2$-projection under the condition 3.2;

(c) the local discrete least-squares fitting under the condition 3.1.

Proof.

The assertion is obvious for the weighted averaging case.

Choose $v = x+y$. Then the minimizer $p_1 = 1$ and $p_2 = 1$ in both cases (b) and (c). Therefore,

$$\begin{equation*} G_hv(z) = (1,1) = \sum _{j=1}^m c_j \nabla (x+y) = \sum _{j=1}^m c_j(1,1). \end{equation*}$$

Now we let $p_l(x,y) = a_0 + a_1x + a_2y$. Then for the local discrete least-squares fitting, $a_i$’s are given by

$$\begin{equation} \begin{pmatrix} m & \sum _j x_j & \sum _j y_j \\ \sum _j x_j & \sum _j x_j^2 & \sum _j x_jy_j \\ \sum _j y_j & \sum _j x_jy_j & \sum _j y_j^2 \end{pmatrix} \begin{pmatrix} a_0 \\ a_1 \\ a_2 \end{pmatrix} = \begin{pmatrix} \sum _j \partial _l u_h(x_j,y_j) \\ \sum _j x_j\partial _l u_h(x_j,y_j) \\ \sum _j y_j\partial _l u_h(x_j,y_j) \end{pmatrix}. {\tag{3.12}} \end{equation}$$

Note that

$$\begin{equation*} \sum _j x_j^2 = O(h^2), \quad \sum _j x_jy_j = O(h^2), \quad \sum _j y_j^2 = O(h^2); \end{equation*}$$

and under condition (3.1),

$$\begin{equation*} \sum _j x_j = O(h^{1+\alpha }), \quad \sum _j y_j = O(h^{1+\alpha }). \end{equation*}$$

By scaling argument we see that

$$\begin{equation*} a_1 = O(h^{\alpha -1}), \quad a_2 = O(h^{\alpha -1}). \end{equation*}$$

Therefore,

$$\begin{eqnarray*} a_0 &=& \frac{1}{m} \sum _j \partial _l u_h(x_j,y_j) - \frac{a_1}{m} \sum _j x_j - \frac{a_2}{m} \sum _j y_j \\ &=& \sum _j c_j \partial _l u_h(x_j,y_j) \end{eqnarray*}$$

with

$$\begin{equation*} c_j = \frac{1}{m} + O(h^{2\alpha }) > 0. \end{equation*}$$

A similar argument shows that

$$\begin{equation*} c_j = \frac{|\tau _j|}{|\omega |} + O(h^{2\alpha }) > 0 \end{equation*}$$

for the local $L^2$-projection when condition (3.2) is satisfied.

■

Under the given condition, the recovered gradient at a vertex $z$ is a convex combination of gradient values on the element patch surrounding $z$.

4. Superconvergence of the recovery operators

We consider the non-self-adjoint problem: find $u\in H^1(\Omega )$ such that

$$\begin{equation} B(u,v) = \int _\Omega [ ({\mathcal{D}} \nabla u + \pmb{b} u)\cdot \nabla v + cuv ] = f(v), \quad \forall v \in H^1(\Omega ). {\tag{4.1}} \end{equation}$$

Here $\mathcal{D}$ is a $2\times 2$ symmetric, positive definite matrix, and $f(\cdot )$ is a linear functional. We assume that all the coefficient functions are smooth, and the bilinear form $B(\cdot ,\cdot )$ is continuous and satisfies the inf-sup condition on $H^1(\Omega )$. These conditions insure that (4.1) has a unique solution.

The finite element solution $u_h \in {\mathcal{V}}_h$ satisfies

$$\begin{equation} B(u_h,v_h) = f(v_h) \quad \forall v_h \in {\mathcal{V}}_h. {\tag{4.2}} \end{equation}$$

To insure a unique solution for (4.2), we further assume the inf-sup condition of $B$ to be satisfied on ${\mathcal{V}}_h$.

We define the piecewise constant matrix function ${\mathcal{D}}_{\tau }$ in terms of the diffusion matrix ${\mathcal{D}}$ as follows:

$$\begin{equation*} {\mathcal{D}}_{\tau ij}=\frac{1}{|\tau |}\int _{\tau } {\mathcal{D}}_{ij}\,dx. \end{equation*}$$

Note that ${\mathcal{D}}_{\tau }$ is symmetric and positive definite.

Theorem 4.1.

Let the solution of 4.1 satisfy $u \in H^3(\Omega ) \cap W^2_\infty (\Omega )$, let $u_h$ be the solution of 4.2 and let $u_I \in {\mathcal{V}}_h$ be the linear interpolation of $u$. Assume that the triangulation ${\mathcal{T}}_h$ satisfies Condition $(\alpha ,\sigma )$. Then

$$\begin{equation*} \|u_h-u_I\|_{1,\Omega } \lesssim h^{1+\rho } (\|u\|_{3,\Omega } + |u|_{2,\infty ,\Omega }), \quad \rho = \min (\alpha , \frac{1}{2}, \frac{\sigma }{2}). \end{equation*}$$

Proof.

We begin with the identity

$$\begin{multline*} B(u-u_I,v_h)=\sum _{\tau \in {\mathcal{T}}_h} \int _{\tau } \nabla (u-u_I)\cdot {\mathcal{D}}_{\tau }\nabla v_h\,dx +\sum _{\tau \in {\mathcal{T}}_h} \int _{\tau } \nabla (u-u_I)\cdot ({\mathcal{D}}-{\mathcal{D}}_{\tau })\nabla v_h\,dx \\ +\int _{\Omega } (u-u_I)( \mathbf{b}\cdot \nabla v_h+c v)\,dx =I_1+I_2+I_3. \end{multline*}$$

The first term $I_1$ is estimated using Lemma 2.1 and $I_2$ and $I_3$ can be easily estimated by

$$\begin{equation*} |I_2|+|I_3|\lesssim h^2|\!| u |\!|_{2,\Omega }|\!| v |\!|_{1,\Omega }. \end{equation*}$$

Thus

$$\begin{equation*} |B(u-u_I,v_h)| \lesssim h^{1+\rho }\left( |\!| u |\!|_{3,\Omega }+| u |_{2,\infty ,\Omega }\right) |\!| v_h |\!|_{1,\Omega }. \end{equation*}$$

We complete the proof using the inf-sup condition in

$$\begin{align*} \hspace{4pc} |\!| u_h-u_I |\!|_{1,\Omega } &\lesssim \sup _{v_h\in {\mathcal{V}}_h} \frac{B(u_h-u_I,v_h)}{|\!| v_h |\!|_{1,\Omega }} =\sup _{v_h\in {\mathcal{V}}_h} \frac{B(u-u_I,v_h)}{|\!| v_h |\!|_{1,\Omega }} \\ &\lesssim h^{1+\rho }\left( |\!| u |\!|_{3,\Omega }+| u |_{2,\infty ,\Omega }\right). \hspace{10pc} \end{align*}$$■

Theorem 4.2.

Let the solution of 4.1 satisfy $u \in W^3_\infty (\Omega )$, let $u_h$ be the solution of 4.2, and let $G_h$ be a recovery operator defined by one of the three: (a) the weighted averaging, (b) the local $L^2$-projection, and (c) the local discrete least-squares fitting. Assume that the triangulation ${\mathcal{T}}_h$ satisfies Condition $(\alpha ,\sigma )$. Then

$$\begin{equation*} \|\nabla u - G_hu_h\|_{0,\Omega } \lesssim h^{1+\rho } \|u\|_{3,\infty ,\Omega }. \end{equation*}$$

Proof.

We decompose

$$\begin{equation} \nabla u - G_hu_h = (\nabla u - (\nabla u)_I) + ((\nabla u)_I - G_hu_I) + G_h(u_I-u_h), {\tag{4.3}} \end{equation}$$

where $(\nabla u)_I \in {\mathcal{V}}_h \times {\mathcal{V}}_h$ is the linear interpolation of $\nabla u$. By the standard approximation theory,

$$\begin{equation} \|\nabla u - (\nabla u)_I\|_{0,\Omega } \lesssim h^2|u|_{3,\Omega }. {\tag{4.4}} \end{equation}$$

We observe that when we pick an element patch on ${\mathcal{T}}_{1,h}$, both conditions (3.1) and (3.2) are satisfied. Therefore, using Theorem 3.1, we have

$$\begin{eqnarray} \|(\nabla u)_I - G_hu_I\|_{0,\Omega _{1,h}} &\le & \left( \sum _{\tau \in \Omega _{1,h}} |\tau | \sum _{z\in {\mathcal{N}}_h\cap \bar{\tau }} |G_hu_I(z)-\nabla u(z)|^2 \right)^{1/2} \\ &\lesssim & h^{1+\alpha } \|u\|_{3,\infty ,\Omega } |\Omega _{1,h}|^{1/2} \lesssim h^{1+\alpha } \|u\|_{3,\infty ,\Omega }. {\tag{4.5}} \end{eqnarray}$$

On the other hand,

$$\begin{equation} \|(\nabla u)_I - G_hu_I\|_{0,\Omega _{2,h}} \lesssim h\|u\|_{3,\infty ,\Omega } |\Omega _{2,h}|^{1/2} \lesssim h^{1+\sigma /2}\|u\|_{3,\infty ,\Omega } {\tag{4.6}} \end{equation}$$

by Condition $(\alpha ,\sigma )$. Combining (4.5) with (4.6), we have

$$\begin{equation} \|(\nabla u)_I - G_hu_I\|_{0,\Omega } \lesssim h^{1+\min (\alpha ,\sigma /2)}\|u\|_{3,\infty ,\Omega }. {\tag{4.7}} \end{equation}$$

Similarly as in (4.5), we have, by using the fact proved in Theorem 3.2, that $G_h v(z)$ is a convex combination of $\nabla v|_{\tau _z}$’s;

$$\begin{eqnarray} \|G_h(u_I-u_h)\|_{0,\Omega _{1,h}} &\le & \left( \sum _{\tau \in {\mathcal{T}}_{1,h}} |\tau | \sum _{z\in {\mathcal{N}}_h\cap \bar{\tau }} |G_h(u_I-u_h)(z)|^2 \right)^{1/2} \\ &\lesssim & \left( \sum _{\tau \in {\mathcal{T}}_{1,h}} |\tau | |\nabla (u_I-u_h)|_\tau |^2 \right)^{1/2} = \|\nabla (u_I-u_h)\|_{0,\Omega _{1,h}} \\ &\lesssim & h^{1+\rho } \|u\|_{3,\infty ,\Omega }, {\tag{4.8}} \end{eqnarray}$$

by Theorem 4.1. In addition,

$$\begin{eqnarray} \|G_h(u_I-u_h)\|_{0,\Omega _{2,h}} &\le & \left( \sum _{\tau \in {\mathcal{T}}_{2,h}} |\tau | \sum _{z\in {\mathcal{N}}_h\cap \bar{\tau }} |G_h(u_I-u_h)(z)|^2 \right)^{1/2} \\ &\lesssim & h\|u\|_{3,\infty ,\Omega } \left( \sum _{\tau \in {\mathcal{T}}_{2,h}} |\tau | \right)^{1/2} \\ &\lesssim & h^{1+\sigma /2} \|u\|_{3,\infty ,\Omega }. {\tag{4.9}} \end{eqnarray}$$

Combining (4.8) and (4.9) yields

$$\begin{equation} \|G_h(u_I-u_h)\|_{0,\Omega } \lesssim h^{1+\rho }\|u\|_{3,\infty ,\Omega }. {\tag{4.10}} \end{equation}$$

The conclusion follows by applying (4.4), (4.7), and (4.10) to the right-hand side of (4.3).

■

Theorem 4.2 requires the global regularity $u\in W^3_\infty (\Omega )$ which is too restrictive in practice. The next theorem turns to interior maximum norm estimates and relaxes the global regularity assumption on the solution.

Theorem 4.3.

Consider an interior patch $\omega _z \subset \subset \Omega _d\subset \Omega _{1,h}$ with $d=\text{dist}(\omega _z,\partial \Omega _d) \ge Kh$ for some constant $K > 0$. Let $u \in W^2_\infty (\Omega )\cap W^3_\infty (\Omega _d)$ be the solution of 4.1, let $u_h$ be the solution of 4.2, and let $G_h$ be a recovery operator defined by one of the three: (a) the weighted averaging, (b) the local $L^2$-projection, and (c) the local discrete least-squares fitting. Then we have

$$\begin{equation*} |(\nabla u - G_hu_h)(z)| \lesssim h^{1+\min (1,\alpha )} \|u\|_{3,\infty ,\omega _z} + d^{-1}h^2\ln \frac{1}{h} \|u\|_{2,\infty ,\Omega } + h^{1+\alpha } \ln \frac{d}{h} |u|_{2,\infty ,\Omega _d}. \end{equation*}$$

Proof.

We denote ${\mathcal{V}}_h^0(\Omega _d)$ as the finite element subspace that has a compact support on $\Omega _d$ and start from

$$\begin{equation*} B(u_h-u_I,\chi ) = B(u-u_I,\chi ) = F(\chi ), \quad \forall \chi \in {\mathcal{V}}_h^0(\Omega _d), \end{equation*}$$

with

$$\begin{equation*} F(\chi ) = \sum _{e\in {\mathcal{E}}_d} \int _e q_e \left\{ (\alpha _e-\alpha '_e)\frac{\partial ^2u}{\partial \pmb{t}^2} +(\beta _e-\beta '_e)\frac{\partial ^2u}{\partial \pmb{t}\partial \pmb{n}} \right\}\frac{\partial \chi }{\partial \pmb{t}}, \end{equation*}$$

where ${\mathcal{E}}_d$ is the edge set of $\Omega _d$. By the same argument as in (2.7), we have

$$\begin{equation*} |F(\chi )| \lesssim h^{1+\alpha } |u|_{2,\infty ,\Omega _d} \int _{\Omega _d} |\nabla \chi |. \end{equation*}$$

Therefore,

$$\begin{equation} |||F|||_{-1,\infty ,\Omega _d} = \sup _{ \chi \in W^1_1(\Omega _d), |\chi |_{W^1_1(\Omega _d)}=1 } F(\chi ) \lesssim h^{1+\alpha } |u|_{2,\infty ,\Omega _d}. {\tag{4.11}} \end{equation}$$

Recall Theorem 1.2 of Schatz-Wahlbin 13 (it is straightforward to verify that all conditions of that theorem are satisfied under the current situation):

$$\begin{eqnarray*} |e|_{W^1_\infty (\Omega _0)} + d^{-1}\|e\|_{L_\infty (\Omega _0)} &\le & C \min _{\chi \in S^h} ( |w-\chi |_{W^1_\infty (\Omega _d)} + d^{-1}\|w-\chi \|_{L_\infty (\Omega _d)} ) \\ &+& Cd^{-1-s-N/q}\|e\|_{W^{-s}_q(\Omega _d)} + C\ln \frac{d}{h} |||F|||_{-1,\infty ,\Omega _d}, \end{eqnarray*}$$

where $e=w-w_h$ satisfies $B(e,\chi ) = F(\chi )$. Now, setting $q=\infty , s=0$, $w=0$, $\Omega _0 = \omega _z$, and $w_h = u_I-u_h$, we obtain

$$\begin{equation*} |u_h-u_I|_{1,\infty ,\omega _z} \lesssim d^{-1} \|u_h-u_I\|_{L_\infty (\Omega _d)} + \ln \frac{d}{h} |||F|||_{-1,\infty ,\Omega _d}. \end{equation*}$$

Applying (4.11) and $\|u_h-u_I\|_{L_\infty (\Omega _d)} \lesssim h^2\ln \frac{1}{h}|u|_{2,\infty ,\Omega }$ results in

$$\begin{equation} |u_h-u_I|_{1,\infty ,\omega _z} \lesssim d^{-1}h^2 \ln \frac{1}{h}\|u\|_{2,\infty ,\Omega } + h^{1+\alpha } \ln \frac{d}{h} |u|_{2,\infty ,\Omega _d}. {\tag{4.12}} \end{equation}$$

Now we decompose

$$\begin{equation*} (\nabla u - G_hu_h)(z) = (\nabla u - G_hu_I)(z) + G_h(u_I-u_h)(z). \end{equation*}$$

By Theorem 3.2, $G_hv(z)$ is a convex combination of values of $\nabla v$ on $\tau \in \omega _z$. Consequently, $G_h$ is a bounded operator in the sense

$$\begin{equation*} |G_hv_h(z)| \lesssim |v_h|_{1,\infty ,\omega _z}, \quad \forall v_h \in {\mathcal{V}}_h. \end{equation*}$$

Therefore,

$$\begin{equation} |(\nabla u - G_hu_h)(z)| \lesssim |(\nabla u - G_hu_I)(z)| + |u_I-u_h|_{1,\infty ,\omega _z}. {\tag{4.13}} \end{equation}$$

The conclusion follows by applying Theorem 3.1 and (4.12) to the right-hand side of (4.13).

■

Remark.

When $\alpha < 1$, we choose $d = h^{1-\alpha }$ and obtain

$$\begin{equation*} |(\nabla u - G_hu_h)(z)| \lesssim h^{1+\alpha }\ln \frac{1}{h}. \end{equation*}$$

When $\alpha \ge 1$, we choose $d = h^{1-\beta }$ with $\beta \in (0,1]$ and obtain

$$\begin{equation*} |(\nabla u - G_hu_h)(z)| \lesssim h^{1+\beta }\ln \frac{1}{h}. \end{equation*}$$

We see that when $\alpha \ge 1$, the recovery is more accurate as $z$ leaves the boundary.

5. Asymptotic exactness of the recovery type error estimators

With preparation in the previous sections, it is now straightforward to prove the asymptotic exactness of error estimators based on the recovery operator $G_h$. The global error estimator is naturally defined by

$$\begin{equation} \eta _h = \|G_hu_h-\nabla u_h\|_{0,\Omega }. {\tag{5.1}} \end{equation}$$

Theorem 5.1.

Assume the hypotheses of Theorem $4.2$. Furthermore, assume that there exists a constant $c(u) > 0$ such that

$$\begin{equation} \|\nabla (u-u_h)\| \ge c(u)h. {\tag{5.2}} \end{equation}$$

Then

$$\begin{equation*} \left| \frac{\eta _h}{\|\nabla (u-u_h)\|_{0,\Omega }} - 1\right| \lesssim h^\rho , \quad \quad \rho = \min (\frac{1}{2},\frac{\sigma }{2},\alpha ). \end{equation*}$$

Proof.

By Theorem 4.2 and hypothesis (5.2), we have

$$\begin{equation*} \left| \frac{\eta _h}{\|\nabla (u-u_h)\|_{0,\Omega }} - 1\right| \le \frac{\|G_hu_h-\nabla u\|_{0,\Omega }}{\|\nabla (u-u_h)\|_{0,\Omega }} \lesssim \frac{h^{1+\rho }\|u\|_{3,\infty ,\Omega }}{c(u)h} \lesssim h^\rho . \end{equation*}$$ ■

The pointwise error estimator at a vertex $z \in \bar{\tau } \subset \Omega _{1,h}$ is naturally defined by

$$\begin{equation} \eta _h^z = |G_hu_h(z)-\nabla u_h(\tau )|. {\tag{5.3}} \end{equation}$$

The next theorem shows that the pointwise error estimator is asymptotically exact.

Theorem 5.2.

Assume the hypotheses of Theorem $4.3$. Let $z$ be a vertex of elements $\tau \subset \Omega _{1,h}$ and assume that there exists a constant $c(u)> 0$ such that

$$\begin{equation} |\nabla u(z) - \nabla u_h(\tau )| \ge c(u)h. {\tag{5.4}} \end{equation}$$

Then we have (a) when $\alpha \in (0,1)$,

$$\begin{equation*} \left| \frac{\eta _h^z}{|\nabla u(z)-\nabla u_h(\tau )|} - 1 \right| \lesssim h^\alpha , \end{equation*}$$

with $\text{dist}(z,\partial \Omega _{1,h}) \ge Kh^{1-\alpha }$; and (b) when $\alpha \ge 1$,

$$\begin{equation*} \left| \frac{\eta _h^z}{|\nabla u(z)-\nabla u_h(\tau )|} - 1 \right| \lesssim h^\beta , \quad \forall \beta \in (0,1], \end{equation*}$$

with $\text{dist}(z,\partial \Omega _{1,h}) \ge Kh^{1-\beta }$.

Proof.

We only prove the case when $\alpha \in (0,1)$. By Theorem 4.3 and hypothesis (5.4), we have

$$\begin{equation*} \left| \frac{\eta _h^z}{|\nabla u(z)-\nabla u_h(\tau )|} - 1 \right| \le \frac{|G_hu_h(z)-\nabla u(z)|}{|\nabla u(z)-\nabla u_h(\tau )|} \lesssim \frac{h^{1+\alpha }}{h} = h^\alpha . \end{equation*}$$ ■

We see that the error estimators (5.1) and (5.3) based on the gradient recovery operator are asymptotically exact under Condition $(\alpha ,\sigma )$. As we mentioned above, this condition is not a very restrictive condition in practice. An automatic mesh generator usually produces some grids which are mildly structured. In practice, a completely unstructured mesh is seldom seen. Our analysis explains in part the good performance of the ZZ error estimator based on the local discrete least-squares fitting for general grids.

Figure 1.

Parameters associated with the triangle $\tau$