Each averaging technique yields reliable a posteriori error control in FEM on unstructured grids. Part I: Low order conforming, nonconforming, and mixed FEM

Abstract

Averaging techniques are popular tools in adaptive finite element methods for the numerical treatment of second order partial differential equations since they provide efficient a posteriori error estimates by a simple postprocessing. In this paper, their reliablility is shown for conforming, nonconforming, and mixed low order finite element methods in a model situation: the Laplace equation with mixed boundary conditions. Emphasis is on possibly unstructured grids, nonsmoothness of exact solutions, and a wide class of averaging techniques. Theoretical and numerical evidence supports that the reliability is up to the smoothness of given right-hand sides.

1. Introduction

Error control and efficient mesh-design in finite element simulations of computational engineering and scientific computing finite element simulations is frequently based on a posteriori error estimates. One of the more popular techniques is local or global averaging, e.g., in form of the ZZ-error indicator ZZ. Efficiency and reliability of this estimator were known only for very structured grids and for solutions of higher regularity and then we have even asymptotic exactness V. Numerical experiments in Baetal showed that averaging techniques were quite more reliable on irregular meshes than expected. For homogeneous Dirichlet conditions and conforming finite element methods, the reliability and efficiency of the ZZ-estimator is proven on unstructured, merely shape-regular grids R2.

This work is devoted to give theoretical and numerical support for the robust reliability of all averaging techniques, robust with respect to violated (local) symmetry of meshes and superconvergence and robust with respect to other boundary conditions or other finite element methods.

For a more precise description of averaging techniques, let us discuss a discretisation of a conservation equation

$$\begin{equation} f+\operatorname {div}p=0 {\tag{1.1}} \end{equation}$$

with a given right-hand side $f\in L^2(\Omega )$ and a known approximation $p_h\in L^2(\Omega )^d$ to the unknown exact flux $p\in H(\operatorname {div};\Omega )$ in the bounded Lipschitz domain $\Omega \subset \mathbb{R}^d$ with piecewise flat boundary. The test function finite element space should include the continuous piecewise affines $\mathcal{S}^1_D(\mathcal{T})$ (with homogeneous Dirichlet boundary conditions) based on a regular (in the sense of Ci, cf. Section 2) triangulation $\mathcal{T}$ of $\Omega$. Suppose a Galerkin property for $p-p_h$ with $\mathcal{S}^1_D(\mathcal{T})$, i.e.,

$$\begin{equation} \int _\Omega p_h\cdot \nabla v_h\,dx=\int _\Omega fv_h\,dx\qquad \text{for all } v_h\in \mathcal{S}^1_D(\mathcal{T}). {\tag{1.2}} \end{equation}$$

What can be said about the error $\Vert p-p_h\Vert _{L^2(\Omega )}$ when we regard $p$ as an unknown and $p_h$ as a known variable?

In averaging techniques, the error estimator is based on a smoother approximation, e.g., in $\mathcal{S}^1(\mathcal{T})^d$, the continuous $\mathcal{T}$-piecewise linears, to the (components of the) discrete solution $p_h$. For instance,

$$\begin{equation} \eta _Z:=\min _{q_h\in \mathcal{S}^1(\mathcal{T})^d}\Vert p_h-q_h\Vert _{L^2(\Omega )} {\tag{1.3}} \end{equation}$$

may serve as a computable error estimator and the elementwise contributions as local error indicators in an adaptive mesh-refining algorithm.

The triangle inequality shows that $\eta _Z$ is efficient with constant $1$ up to higher order terms of the exact solution $p$, indeed,

$$\begin{equation} \eta _Z\le \Vert p-p_h\Vert _{L^2(\Omega )}+ \min _{q_h\in \mathcal{S}^1(\mathcal{T})^d}\Vert p-q_h\Vert _{L^2(\Omega )}. {\tag{1.4}} \end{equation}$$

The last term converges as $O(h^{2})$ (provided $p$ is smooth enough and $h$ denotes the maximal mesh–size in $\mathcal{T}$) and so, generically, is of higher order than the error $\Vert p-p_h\Vert _{L^2(\Omega )}=O(h)$ in the lowest order finite element method. If the second term in the right-hand side of 1.4 fails to be of higher order, one can still prove efficiency of $\eta _Z$ using equivalence of global and local averaging (cf. Theorem 3.2) and that local averaging is equivalent to weighted jumps across interelement boundaries. An efficiency estimate with higher order terms that depend on local smoothness of right-hand sides but with unknown constants then follows as in V.

In practise, we may apply an averaging operator $\mathcal{A}:L^2(\Omega )^d\to \mathcal{S}^{1}(\mathcal{T})^d$ to $p_h$ and compute the upper bound $\Vert p_h-\mathcal{A}p_h\Vert _{L^2(\Omega )}$ of $\eta _Z$. Then, efficiency depends strongly on the approximation properties of $\mathcal{A}$ and deserves further investigation.

In this paper, the focus is on the reliability of $\eta _Z$, i.e., we investigate under which conditions an estimate

$$\begin{equation} \Vert p-p_h\Vert _{L^2(\Omega )}\le c_1\,\eta _Z+\text{ h.o.t.} \le c_1\,\Vert p_h-\mathcal{A}p_h\Vert _{L^2(\Omega )}+\text{ h.o.t.} {\tag{1.5}} \end{equation}$$

holds, we study what the constant $c_1>0$ depends on, what affects the higher-order contributions “$\text{h.o.t.}$”, how to modify the definition of $\eta _Z$ in the presence of mixed boundary conditions, and how to modify the general setting presented for nonconforming and mixed lowest order finite element methods.

Recall from 1.5 that any averaging technique, described by $\mathcal{A}$, then is reliable up to higher order terms. We also prove equivalence to local modifications of $\eta _Z$ where the minimisation is over smaller domains, e.g., patches of nodes or edges.

The outline of the paper is as follows. Preliminaries and notation are introduced in Section 2 where we state and prove stability and first order estimates for a certain approximation operator $\mathcal{J}:H^1_D(\Omega )\to \mathcal{S}^1_D(\mathcal{T})$ essentially designed to yield further local orthogonality properties as in CVC2. Basic estimates are provided in Section 3 for a local and global averaging technique and their equivalence. The subsequent Sections 4, 5, and 6 display the consequences to averaging techniques in a posteriori error control for first order conforming, nonconforming and mixed finite element schemes. Numerical evidence, reported in Section 7, supports the theoretical results for adaptively refined and evenly perturbed meshes. Although asymptotic exactness is not claimed in this paper, our numerical experiments illustrate that $\eta _Z$ is a very good approximate to $\|p-p_h\|_{L^2(\Omega )}$ even on perturbed grids.

The proofs are given for a simple elliptic model example with mixed boundary conditions for conforming, nonconforming, and mixed finite elements in two dimensions for notational simplicity. More interesting examples such as higher-order schemes, the application to the Stokes problem or the Navier-Lamé equations without incompressibility-locking will appear elsewhere BCCF2CF3CF4CF5.

2. Approximation in finite element spaces

The Lipschitz boundary $\Gamma =\partial \Omega$ of the bounded domain $\Omega$ is split into a closed Dirichlet part $\Gamma _D$ with positive surface measure and a remaining, relatively open and possibly empty, Neumann part $\Gamma _N:=\Gamma \setminus \Gamma _D$. Suppose $\mathcal{T}$ be a regular triangulation of the domain $\Omega \subseteq \mathbb{R}^d$, $d=1,2,3$, in the sense of Ciarlet BSCi (no hanging node, domain is matched exactly) with piecewise affine Lipschitz boundary $\Gamma =\partial \Omega =\Gamma _D\cdot \hspace{-8pt}\cup \Gamma _N$, i.e., $\mathcal{T}$ consists of a finite number of closed subsets of $\overline{\Omega }$, that cover $\overline{\Omega }=\cup \mathcal{T}$. Each element $T\in \mathcal{T}$ is either an interval $T=\operatorname {conv}\{a,b\}$ if $d=1$, a triangle $T=\operatorname {conv}\{a,b,c\}$, or a parallelogram $T=\operatorname {conv}\{a,b,c,d\}$ if $d=2$. The extremal points $a,b,c$ are called vertices, the faces $E\subseteq \partial T$, e.g., $E=\operatorname {conv}\{a,b\}$, are called edges. The set of all vertices and all edges appearing for some $T$ in $\mathcal{T}$ are denoted as $\mathcal{N}$ and $\mathcal{E}$. Two distinct and intersecting $T_1$ and $T_2$ share either an entire edge or a vertex. Each edge $E\in \mathcal{E}$ on the boundary $\Gamma$ belongs either to $\Gamma _D$, written $E\in \mathcal{E}_D$, or to $\overline{\Gamma }_N$, written $E\in \mathcal{E}_N$. Therefore the set of edges is partitioned into $\mathcal{E}_\Omega :=\{E\in \mathcal{E}: E\not \subset \Gamma \}$, $\mathcal{E}_D$, and $\mathcal{E}_N$. We stress that $\cup \mathcal{E}$, the union of all egdes, denotes the skeleton of egdes in $\mathcal{T}$, i.e., the set of all points $x$ that belong to some boundary $x\in \partial T$ of some element $T\in \mathcal{T}$. Finally, $\mathcal{K}:=\mathcal{N}\setminus \Gamma _D$ denotes the set of free nodes.

We do not explicitly distinguish between nodes and vertices when we consider conforming finite elements (and avoid these concepts for nonconforming or mixed finite element schemes).

For $T\in \mathcal{T}$, let $P_T^k:=\mathcal{P}_k(T)$ if $T$ is a triangle or $P_T^k:=\mathcal{Q}_k(T)$ if $T$ is a parallelogram. Here, $\mathcal{P}_k(K)$, resp. $\mathcal{Q}_k(K)$, denotes the set of algebraic polynomials in $d$ variables on $K$ of total, resp. partial, degree $\leq k$. The space $\mathcal{L}^k(\mathcal{T})$ of (possibly discontinuous) $\mathcal{T}$-piecewise polynomials of degree $\le k$ is the set of all $U\in L^{\infty }(\Omega )$ with $U|_T\in P_T^k$ for all $T$ in $\mathcal{T}$. Set

$$\begin{equation*} \mathcal{S}^k(\mathcal{T}):=\mathcal{L}^k(\mathcal{T})\cap C(\Omega )\quad \text{and}\quad \mathcal{S}^1_D(\mathcal{T}):=\{u_h\in \mathcal{S}^1(\mathcal{T}):u_h|_{\Gamma _D}=0\}. \end{equation*}$$

Let $(\varphi _z|z\in \mathcal{N})$ denote the nodal basis of $\mathcal{S}^1(\mathcal{T})$, i.e., $\varphi _z\in \mathcal{S}^1(\mathcal{T})$ satisfies $\varphi _z(x)=0$ if $x\in \mathcal{N}\setminus \{z\}$ and $\varphi _z(z)=1$. Note that $(\varphi _z|z\in \mathcal{N})$ is a partition of unity and the open patches

$$\begin{equation} \omega _z:=\{x\in \Omega :0<\varphi _z(x)\} \tag{2.1} \end{equation}$$

form an open cover $(\omega _z:z\in \mathcal{N})$ of $\Omega$ with finite overlap.

In order to define a weak interpolation operator $\mathcal{J}:H^1_D(\Omega )\to \mathcal{S}^1_D(\mathcal{T})$, we modify $(\varphi _z|z\in \mathcal{K})$ to a partition of unity $(\psi _z|z\in \mathcal{K})$. For each fixed node $z\in \mathcal{N}\setminus \mathcal{K}$, we choose a node $\zeta (z)\in \mathcal{K}$ and let $\zeta (z):=z$ if $z\in \mathcal{K}$. In this way, we define a partition of $\mathcal{N}$ into $\operatorname {card}(\mathcal{K})$ classes $I(z):=\{\tilde{z}\in \mathcal{N}: \zeta (\tilde{z})=z\}$, $z\in \mathcal{K}$. For each $z\in \mathcal{K}$ set

$$\begin{equation} \psi _z:=\sum _{\zeta \in I(z)}\varphi _{\zeta } \tag{2.2} \end{equation}$$

and notice that $(\psi _z|z\in \mathcal{K})$ is a partition of unity. It is required that

$$\begin{equation} \Omega _z:=\{x\in \Omega : 0<\psi _z(x)\} \tag{2.3} \end{equation}$$

is connected and that $\psi _z\ne \varphi _z$ implies that $\Gamma _D\cap \partial \Omega _z$ has a positive surface measure.

For $g\in L^1(\Omega )$ and $z\in {\mathcal{K}}$ let $g_z\in \mathbb{R}$ be

$$\begin{equation} g_z := \frac{\int _{\Omega _z}g\psi _z\,dx}{\int _{\Omega _z}\varphi _z\,dx}, {\tag{2.4}} \end{equation}$$

and then define

$$\begin{equation} {\mathcal{J}}g:=\sum _{z\in {\mathcal{K}}}g_z\varphi _z\in \mathcal{S}^1_D(\mathcal{T}). {\tag{2.5}} \end{equation}$$

The local mesh-sizes are denoted by $h_\mathcal{T}$ and $h_\mathcal{E}$, where $h_\mathcal{T}\in \mathcal{L}^0(\mathcal{T})$ denotes the element-size, $h_\mathcal{T}|_T:=h_T:=\operatorname {diam}(T)$ for $T\in \mathcal{T}$, and the edge-size $h_\mathcal{E}\in L^\infty (\cup \mathcal{E})$ is defined on the union or skeleton $\cup \mathcal{E}$ of all edges $E$ in $\mathcal{E}$ by $h_\mathcal{E}|_E:=h_E:=\operatorname {diam}(E)$. The patch-size $h_z:=\operatorname {diam}(\Omega _z)$ is defined for each node $z\in \mathcal{K}$ separately.

Theorem 2.1.

There exist $(h_\mathcal{T},h_\mathcal{E})$-independent constants $c_2,c_3,c_4,c_5>0$ such that for all $g\in H^1_D(\Omega )$ and $f\in L^2(\Omega )$ there holds

$$\begin{align} \|\nabla {\mathcal{J}}g - \nabla g\|_{L^2(\Omega )} &\le c_2\|\nabla g\|_{L^2(\Omega )}, {\tag{2.6}}\\ \int _\Omega f(g-{\mathcal{J}}g) \, dx &\le c_3 \|\nabla g\|_{L^2(\Omega )} \bigl (\sum _{z\in {\mathcal{K}}} h_z^2 \min _{f_z\in \mathbb{R}} \|f-f_z\|_{L^2(\Omega _z)}^2 \bigr )^{1/2}, {\tag{2.7}}\\ \|h_{\mathcal{T}}^{-1} (g-{\mathcal{J}}g)\|_{L^2(\Omega )} &\le c_4\|\nabla g\|_{L^2(\Omega )}, {\tag{2.8}}\\ \|h_{\mathcal{E}}^{-1/2} (g-{\mathcal{J}}g)\|_{L^2(\Gamma _N)} &\le c_5\|\nabla g\|_{L^2(\Omega )}. {\tag{2.9}} \end{align}$$

The constants $c_2,c_3,c_4,c_5$ only depend on $\Omega$, $\Gamma _D$, $\Gamma _N$ and the shape of the elements and patches (not on their sizes).

Remark 2.1.

The assertion of the theorem holds verbatim for three space dimensions where $\mathcal{T}$ consists of tetrahedra or parallelepipeds with the same proof.

Proof.

In this proof and at similar occasions, $\lesssim$ abbreviates an inequality $\le$ up to a constant $(h_\mathcal{T},h_\mathcal{E})$-independent factor. Also, $\|\cdot \|_{p,K}$ abbreviates $\|\cdot \|_{L^p(K)}$ and we neglect $K$ if $\Omega$ is meant, i.e., $\|\cdot \|_2 := \|\cdot \|_{2,\Omega }$. Hence, e.g., 2.6 could be phrased as $\|\nabla {\mathcal{J}}g-\nabla g\|_2\lesssim \|\nabla g\|_2$.

The key estimate for the stability and the approximation property of ${\mathcal{J}}$ will be

$$\begin{equation} \|g_z\varphi _z-g\psi _z\|_{2,\Omega _z}\lesssim h_z\|\nabla g\|_{2,\Omega _z}\qquad (z\in {\mathcal{K}}). {\tag{2.10}} \end{equation}$$

For the proof of 2.10, let $\overline{g}_z$ denote the integral mean of $g$ on $\Omega _z$. Then, using the definition 2.4 for the coefficients $g_z$, Cauchy’s and Young’s inequality, we infer, with $c_6 :=\|1\|_{2,\Omega _z}/\|\varphi _z^{1/2}\|_{2,\Omega _z}$,

$$\begin{equation} \begin{split} \|\varphi _z^{1/2}(g_z-\overline{g}_z)\|_{2,\Omega _z}^2 &= \int _{\Omega _z} \varphi _z (\overline{g}_z -g)(\overline{g}_z -g_z)\,dx + \int _{\Omega _z} (\psi _z - \varphi _z) g (g_z-\overline{g}_z) \,dx\\ &\le \frac{1}{4}\|\varphi _z^{1/2}(g_z-\overline{g}_z)\|_{2,\Omega _z}^2 +\|g-\overline{g}_z\|_{2,\Omega _z}^2 +\frac{1}{4 c_6^2}\| g_z-\overline{g}_z\|_{2,\Omega _z}^2\\ &\quad +c_6^2\|(\psi _z-\varphi _z)g\|_{2,\Omega _z}^2. \end{split} {\tag{2.11}} \end{equation}$$

Absorbing $\frac{1}{4 c_6^2}\| g_z-\overline{g}_z\|_{2,\Omega _z}^2 \le \frac{1}{4} \|\varphi _z^{1/2}(g_z-\overline{g}_z)\|_{2,\Omega _z}^2$, we deduce

$$\begin{equation} \|\varphi _z^{1/2}(g_z-\overline{g}_z)\|_{2,\Omega _z}^2 \lesssim \|g-\overline{g}_z\|_{2,\Omega _z}^2 + \|(\psi _z-\varphi _z)g\|_{2,\Omega _z}^2. {\tag{2.12}} \end{equation}$$

A Poincaré inequality yields

$$\begin{equation*} \|g-\overline{g}_z\|_{2,\Omega _z}\lesssim h_z\|\nabla g\|_{2,\Omega _z}. \end{equation*}$$

Note that $(\psi _z-\varphi _z)g$ is nonzero only if $\Gamma _D\cap (\partial \Omega _z)$ has positive surface measure. Since $g$ vanishes there, Friedrichs’ inequality shows

$$\begin{equation*} \|g\|_{2,\Omega _z} \lesssim h_z \|\nabla g\|_{2,\Omega _z}. \end{equation*}$$

Therefore, (2.12) results in

$$\begin{equation} \|\varphi _z^{1/2}(g_z-\overline{g}_z)\|_{2,\Omega _z} \lesssim h_z\|\nabla g\|_{2,\Omega _z}. {\tag{2.13}} \end{equation}$$

To prove (2.10), we use the triangle inequality, (2.13), and again Cauchy’s and Friedrichs’ inequality to verify

$$\begin{multline} \|g_z\varphi _z-g\psi _z\|_{2,\Omega _z}\le \|(g_z-\overline{g}_z)\varphi _z\|_{2,\Omega _z} +\|(g-\overline{g}_z)\varphi _z\|_{2,\Omega _z}\\ +\|(\psi _z-\varphi _z)g\|_{2,\Omega _z} \lesssim h_z\|\nabla g\|_{2,\Omega _z}. {\tag{2.14}} \end{multline}$$

To prove (2.7), we use that $(\psi _z|z\in {\mathcal{K}})$ is a partition of unity and obtain with (2.10), (2.4) for any $f_z\in \mathbb{R}$ that

$$\begin{multline} \int _\Omega f(g-{\mathcal{J}}g)\,dx=\sum _{z\in {\mathcal{K}}} \int _{\Omega _z}f(g\psi _z-g_z\varphi _z)\,dx=\sum _{z\in \mathcal{K}} \int _{\Omega _z}(f-f_z)(g\psi _z-g_z\varphi _z)\,dx \\ \lesssim \sum _{z\in {\mathcal{K}}}\|f-f_z\|_{2,\Omega _z}h_z\|\nabla g\|_{2,\Omega _z} \lesssim \bigl (\sum _{z\in {\mathcal{K}}}h_z^2\|f-f_z\|_{2,\Omega _z}^2\bigr )^{1/2} \|\nabla g\|_{2,\Omega }. {\tag{2.15}} \end{multline}$$

In the last step we used that $(\varphi _z|z\in {\mathcal{K}})$ has a finite overlap that depends on the shape of the elements only. This concludes the proof of (2.7).

The remaining part of the proof uses standard arguments and is therefore sketched for brevity. To prove (2.8) we let $f:=h_{\mathcal{T}}^{-2}(g-{\mathcal{J}}g)$ and $f_z=0$, $z\in {\mathcal{K}}$, in (2.7). To verify (2.6) we use $\sum _{z\in {\mathcal{K}}} \psi _z =1$ and $\sum _{z\in {\mathcal{K}}} \nabla \psi _z =0$ and repeat the triangle inequality several times for

$$\begin{equation} \|\nabla g - \nabla {\mathcal{J}}g\|_{2}^2 \lesssim \sum _{z\in {\mathcal{K}}} \|\nabla (\psi _z g-\varphi _z g_z)\|_{2}^2. {\tag{2.16}} \end{equation}$$

With Friedrichs’ and Poincaré’s inequality we infer

$$\begin{multline} \|\nabla (g\psi _z-g_z\varphi _z)\|_{2,\Omega _z} \le \|(\psi _z-\varphi _z)\nabla g\|_{2,\Omega _z} + \|\nabla (\varphi _z(g_z-\overline{g}_z))\|_{2,\Omega _z} \\ +\|\nabla (\varphi _z(\overline{g}_z-g))\|_{2,\Omega _z} + \|g\nabla (\psi _z-\varphi _z)\|_{2,\Omega _z} \lesssim \|\nabla g\|_{2,\Omega _z} +\|\nabla (\varphi _z(g_z-\overline{g}_z))\|_{2,\Omega _z}. {\tag{2.17}} \end{multline}$$

It remains to estimate $\|\nabla (\varphi _z(g_z-\overline{g}_z))\|_{2,\Omega _z}$ with 2.10, Friedrichs’ inequality, and the above arguments. A trace inequality BSClCF1 of the form

$$\begin{equation} \Vert w\Vert _{2,E}\lesssim h_E^{-1/2}\Vert w\Vert _{2,T}+ h_E^{1/2}\Vert \nabla w\Vert _{2,T} {\tag{2.18}} \end{equation}$$

for $E\in \mathcal{E}_N$ and $T\in \mathcal{T}$ with $E\subset \partial T\cap \overline{\Gamma }_N$ together with 2.6 and 2.8 implies 2.9.

■

3. Basic estimates

In this section we first derive with the approximation operator $\mathcal{J}$ a global error estimate for a posteriori error control by averaging processes in an abstract setting. We then show the equivalence of local and global averaging techniques. The estimates of this section are then specified, and thereby proved to be substantial, in the subsequent sections to conforming, nonconforming, and mixed finite element methods.

Theorem 3.1.

Suppose $p,q\in H(\operatorname {div};\Omega )$ and $p_h\in \mathcal{L}^k(\mathcal{T})^d$ with $p\cdot n,q\cdot n\in L^2(\Gamma _N)$ and

$$\begin{equation} \int _\Omega (p-p_h)\cdot \nabla w_h\, dx=0\quad \text{for all }w_h\in \mathcal{S}^1_D(\mathcal{T}). {\tag{3.1}} \end{equation}$$

Then there holds

$$\begin{equation} \begin{split} &\sup _{\substack{w\in H^1_D(\Omega )\\ \Vert \nabla w\Vert _{L^2(\Omega )}=1}} \int _\Omega (p-p_h)\cdot \nabla w\, dx \\ &\qquad \le c_2\Vert p_h-q\Vert _{L^2(\Omega )} +c_3\bigl (\sum _{z\in {\mathcal{K}}} h_z^2 \min _{f_z\in \mathbb{R}}\Vert \operatorname {div}(p-q)-f_z\Vert _{L^2(\Omega _z)}^2\bigr )^{1/2} \\ &\qquad \quad +c_5 \Vert h_\mathcal{E}^{1/2}\,(p-q)\cdot n\Vert _{L^2(\Gamma _N)}. \end{split} {\tag{3.2}} \end{equation}$$

Proof.

According to 3.1, 2.6, Cauchy’s inequality, and an integration by parts we have, for each $w\in H^1_D(\Omega )$ with $\Vert \nabla w\Vert _{L^2(\Omega )}=1$, that

$$\begin{multline} \int _\Omega (p-p_h)\cdot \nabla w\,dx=\int _\Omega (p-p_h)\cdot \nabla (w-\mathcal{J}w)\,dx\\ =\int _\Omega (p-q)\cdot \nabla (w-\mathcal{J}w)\,dx +\int _\Omega (q-p_h)\cdot \nabla (w-\mathcal{J}w)\,dx\\ \le \int _{\Gamma _N}(w-\mathcal{J}w)\,(p-q)\cdot n\,dx -\int _\Omega (w-\mathcal{J}w)\,\operatorname {div}(p-q)\,dx +c_2\Vert p_h-q\Vert _{2,\Omega } {\tag{3.3}} \end{multline}$$

since $w$ and $\mathcal{J}w$ vanish on $\partial \Omega \setminus \Gamma _N$. Owing to 2.7 and 2.9 in Theorem 2.1, we conclude 3.2 from 3.3 and Cauchy’s inequality.

■

The second result justifies local averaging. For each edge $E\in \mathcal{E}_\Omega$, let $\omega _E:=\operatorname {int}(T_1\cup T_2)$ and $\mathcal{T}_E:=\{T_1,T_2\}$ for the two distinct elements $T_1,T_2\in \mathcal{T}$ with $E=T_1\cap T_2$ and for each edge $E\in \mathcal{E}_N$, let $\omega _E:=\operatorname {int}(T)$ and $\mathcal{T}_E:=\{T\}$ for the element $T\in \mathcal{T}$ with $E=T\cap \overline{\Gamma }_N$. Let $\mathcal{L}^k(\mathcal{E}_N)$ denote the (possibly discontinuous) $\mathcal{E}_N$-piecewise polynomials of degree $\le k$ on $\Gamma _N$ and let $\mathcal{S}^k(\mathcal{T}_E):=\mathcal{L}^k(\mathcal{T}_E)\cap C(\omega _E)$.

Theorem 3.2.

There exists an $(h_\mathcal{T},h_\mathcal{E})$-independent constant $c_7>0$ which depends on the shape of the elements in $\mathcal{T}$ and on the polynomial degree $k\ge 1$, $c_8= \max _{T\in \mathcal{T}} \operatorname {card}\{E\in \mathcal{E}_\Omega \cup \mathcal{E}_N : E\subseteq \partial T\}$, such that, for all $(p_h,g_h)\in \mathcal{L}^{k-1}(\mathcal{T})^d\times \mathcal{L}^k(\mathcal{E}_N)$, we have

$$\begin{equation} \begin{split} &c_7 \min _{q_h \in \mathcal{S}^k(\mathcal{T})^d}\Bigl (\|p_h-q_h\|_{L^2(\Omega )}^2 +\|h_\mathcal{E}^{1/2}(g_h-q_h\cdot n)\|_{L^2(\Gamma _N)}^2\Bigr )\\ &\qquad \le \sum _{E\in \mathcal{E}_\Omega \cup \mathcal{E}_N}\min _{q_E\in \mathcal{S}^k(\mathcal{T}_E)^d} \Bigl (\|p_h-q_E\|_{L^2(\omega _E)}^2 +h_E\|g_h-q_E\cdot n\|_{L^2(E\cap \Gamma _N)}^2\Bigr ) \\ &\qquad \le c_8 \min _{q_h \in \mathcal{S}^k(\mathcal{T})^d} \Bigl (\|p_h-q_h\|_{L^2(\Omega )}^2 +\|h_\mathcal{E}^{1/2}(g_h-q_h\cdot n)\|_{L^2(\Gamma _N)}^2\Bigr ) . \end{split} {\tag{3.4}} \end{equation}$$

Proof.

The upper estimate follows from $q_E:=q_h|_{\omega _E} \in \mathcal{S}^k(\mathcal{T}_E)^d$ for all $q_h\in \mathcal{S}^k(\mathcal{T})^d$ and a rearrangement of the sums over edges and elements. To verify the lower estimate in 3.4 we consider a subspace $\tilde{\mathcal{S}}^k(\mathcal{T})$ of $\mathcal{S}^k(\mathcal{T})$,

$$\begin{equation*} \tilde{\mathcal{S}}^k(\mathcal{T}):=\{\sum _{z\in \mathcal{N}} q_z \varphi _z: \, q_z \in \mathcal{S}^{k-1}(\mathcal{T}_z) \} \subseteq \mathcal{S}^k(\mathcal{T}), \end{equation*}$$

where $\mathcal{T}_z=\{T\in \mathcal{T}:\,T\subseteq \overline{\omega }_z\}$ denotes the restriction of the triangulation $\mathcal{T}$ to $\omega _z$. Since $\{(q_h,q_h\cdot n|_{\Gamma _N}):q_h\in \tilde{\mathcal{S}}^k(\mathcal{T})^d\}$ is a closed convex subset of $L^2(\Omega )^d \times L^2(\Gamma _N)$, the best-approximation problem

$$\begin{equation} \min _{q_h \in \tilde{\mathcal{S}}^k(\mathcal{T})^d}\Bigl (\|p_h-q_h\|_{2}^2 +\|h_\mathcal{E}^{1/2}(g_h-q_h\cdot n)\|_{2,\Gamma _N}^2\Bigr ) {\tag{3.5}} \end{equation}$$

defines an orthogonal relation, namely, for all ${q}_z\in \mathcal{S}^{k-1}(\mathcal{T}_z)^d$,

$$\begin{equation} \int _\Omega (p_h-\tilde{q_h})\cdot {q}_z\varphi _z\,dx+ \int _{\Gamma _N}h_\mathcal{E}(g_h-\tilde{q}_h\cdot n)\,{q}_z\cdot n\varphi _z\,ds=0, {\tag{3.6}} \end{equation}$$

where $\tilde{q}_h=\sum _{z\in \mathcal{N}}\tilde{q}_z \varphi _z \in \tilde{\mathcal{S}}^k (\mathcal{T})^d$, $\tilde{q}_z\in \mathcal{S}^{k-1}(\mathcal{T}_z)^d$, denotes the minimiser in 3.5. From $\sum _{z\in \mathcal{N}}\varphi _z=1$, (3.6), and Cauchy’s inequality we deduce, for arbitrary ${q}_z\in \mathcal{S}^{k-1}(\mathcal{T}_z)^d$,

$$\begin{equation} \begin{split} \|p_h&-\tilde{q}_h\|_{2}^2 +\|h_\mathcal{E}^{1/2}(g_h-\tilde{q}_h\cdot n)\|_{2,\Gamma _N}^2 \\ =&\sum _{z\in \mathcal{N}}\Bigl (\int _\Omega (p_h-\tilde{q}_h)\varphi _z(p_h-\tilde{q}_z)\,dx +\int _{\Gamma _N}h_\mathcal{E}(g_h-\tilde{q}_h\cdot n) \varphi _z(g_h-\tilde{q}_z\cdot n)ds\Bigr )\\ =&\sum _{z\in \mathcal{N}}\Bigl (\int _\Omega (p_h-\tilde{q}_h)\varphi _z(p_h - q_z)dx +\int _{\Gamma _N}h_\mathcal{E}(g_h-\tilde{q}_h\cdot n) \varphi _z(g_h-q_z\cdot n)\,ds\Bigr )\\ \le &(\|p_h-\tilde{q}_h\|_{2} +\|h_\mathcal{E}^{1/2}(g_h-\tilde{q}_h\cdot n)\|_{2,\Gamma _N}) \\ &\times \Bigl (\sum _{z\in \mathcal{N}} \Bigl (\|\varphi _z^{1/2}(p_h-q_z)\|_{2,\omega _z}^2 + \|h_\mathcal{E}^{1/2}\varphi _z^{1/2}(g_h-q_z\cdot n) \|_{2,\Gamma _N}^2\Bigr )\Bigr )^{1/2}. \end{split} {\tag{3.7}} \end{equation}$$

For each $z\in \mathcal{N}$, we consider the semi-norms on a finite dimensional subspace of $L^2(\omega _z)^d\times L^2((\partial \omega _z) \cap \Gamma _N)$

$$\begin{eqnarray*} |\|(p_h,g_h)\||_{z,1}&:=\displaystyle \min _{q_z\in \mathcal{S}^{k-1}(\mathcal{T}_z)^d}\bigl (\|\varphi _z^{1/2}(p_h-q_z)\|_{2,\omega _z} +\|\varphi _z^{1/2}h_\mathcal{E}^{1/2}(g_h-q_z\cdot n) \|_{2,\Gamma _N}\bigr ), \\ |\|(p_h,g_h)\||_{z,2}&:=\displaystyle \Bigl (\sum _{\substack{E\in \mathcal{E}\\z\in E}} \min _{q_E\in \mathcal{S}^k(\mathcal{T}_E)^d}\Bigl (\|p_h-q_E\|_{2,\omega _E}^2 +h_E\|g_h-q_E\cdot n\|_{2,\Gamma _N\cap E}^2\Bigr )\Bigr )^{1/2}. \end{eqnarray*}$$

Then, 3.7 and $\tilde{\mathcal{S}}^k(\mathcal{T})\subseteq \mathcal{S}^k(\mathcal{T})$ yield

$$\begin{equation} \min _{q_h \in \mathcal{S}^k(\mathcal{T})^d}\bigl (\|p_h-q_h\|_{2}^2 +\|h_\mathcal{E}^{1/2}(g_h-q_h\cdot n)\|_{2,\Gamma _N}^2) \lesssim \sum _{z\in \mathcal{N}} |\|(p_h,g_h)\||_{z,1}^2. {\tag{3.8}} \end{equation}$$

We claim $|\|\cdot \||_{z,1}\lesssim |\|\cdot \||_{z,2}$. For a proof, suppose $|\|(p_h,g_h)\||_{z,2} = 0$. Then, for each $E$ that is an inner edge of $\omega _z$, we have $p_h=q_E$ on the open set $\omega _E$ for some $q_E\in \mathcal{S}^k(\mathcal{T}_E)^d$. Since $p_h\in \mathcal{L}^{k-1}(\mathcal{T})^d$, we find that $p_h|_{\omega _E}\in \mathcal{S}^{k-1}(\mathcal{T}_E)$. The set of all such $\omega _E$ is a cover of $\omega _z$ and there is a sequence $E_1,...,E_m$ of inner edges such that $\omega _{E_j}\cap \omega _{E_{j+1}} \neq \emptyset$, so that we deduce $p_h|_{\omega _z} \in \mathcal{S}^{k-1}(\mathcal{T}_z)$. Moreover, $g_h = p_h\cdot n$ on each edge $E\subseteq \overline{\Gamma }_N$ with $z\in E$, while for edges $E\subseteq \partial \omega _z\cap \overline{\Gamma }_N$ with $z\not \in E$ we have $\varphi _z|_E=0$. Altogether, we deduce $|\|(p_h,g_h)\||_{z,1}= 0$. A compactness and scaling argument then shows our claim

$$\begin{equation} |\|\cdot \||_{z,1}\lesssim |\|\cdot \||_{z,2}\quad \text{on } \mathcal{L}^{k-1}(\mathcal{T}_z)^d\times \mathcal{L}^k(\{E\in \mathcal{E}: E\subseteq \partial \omega _z\}). {\tag{3.9}} \end{equation}$$

Utilizing 3.9 in 3.8, we conclude

$$\begin{equation} \begin{split} &\min _{q_h \in \mathcal{S}^k(\mathcal{T})^d}\bigl (\|p_h-q_h\|_{2}^2 +\|h_\mathcal{E}^{1/2}(g_h-q_h\cdot n)\|_{2,\Gamma _N}^2)\\ &\qquad \lesssim \sum _{z\in \mathcal{N}}|\|(p_h,g_h)\||_{z,1}^2 \lesssim \sum _{z\in \mathcal{N}}|\|(p_h,g_h)\||_{z,2}^2\\ &\qquad \lesssim \sum _{\substack{E\in \mathcal{E}}}\min _{q_E\in \mathcal{S}^k(\mathcal{T}_E)^d} \bigl (\|p_h-q_E\|_{2,\omega _E}^2 +h_E\|g_h-q_E\cdot n\|_{2,\Gamma _N\cap E}^2\bigr ).\end{split} {\tag{3.10}} \end{equation}$$■

Remark 3.1.

The assertions of Theorems 3.1 and 3.2 hold verbatim for three space dimensions where $\mathcal{T}$ consists of tetrahedra or parallelepipeds with the same proofs.

4. Applications to conforming finite element schemes

Given right-hand sides $f\in L^2(\Omega )$, $g\in L^2(\Gamma _N)$, and $u_D\in H^1(\Gamma _D)$, let $u\in H^1(\Omega )$ denote the unique weak solution to

$$\begin{alignat}{2} -\Delta u&=f &\qquad &\text{in }\Omega , {\tag{4.1}}\\ u&= u_D &\qquad &\text{on }\Gamma _D, {\tag{4.2}}\\ \partial u/\partial n&=g&\qquad &\text{on }\Gamma _N. {\tag{4.3}} \end{alignat}$$

Suppose a finite element scheme, based on a regular triangulation $\mathcal{T}$, provided a discrete flux $p_h:=\nabla u_h$ to the exact flux $p:=\nabla u\in H(\operatorname {div};\Omega )$ such that $u_h\in \mathcal{S}^1(\mathcal{T})$, $u_h(z) =u_D(z)$ for all $z\in \mathcal{N}\cap \Gamma _D$ and

$$\begin{equation} \int _\Omega \nabla u_h\cdot \nabla w_h\, dx= \int _\Omega f w_h\,dx+\int _{\Gamma _N} g w_h\,ds \quad \text{for all }w_h\in \mathcal{S}^1_D(\mathcal{T}). {\tag{4.4}} \end{equation}$$

Theorem 4.1.

There exists an $(h_\mathcal{T},h_\mathcal{E})$-independent constant $c_9>0$ (that depends on $k$ and the shape of the elements and patches) such that

$$\begin{multline} \Vert \nabla (u-u_h)\Vert _{L^2(\Omega )}\le \min _{q_h\in \mathcal{S}^k(\mathcal{T})^d}\bigl (c_9\Vert \nabla u_h-q_h\Vert _{L^2(\Omega )} +2c_5\Vert h_\mathcal{E}^{1/2}(g-q_h\cdot n)\Vert _{L^2(\Gamma _N)}\bigr ) \\ +\inf _{v|_{\Gamma _D}=u_D} \Vert \nabla (u_h-v)\Vert _{L^2(\Omega )} +2c_3\bigl (\sum _{z\in {\mathcal{K}}} h_z^2 \min _{f_z\in \mathbb{R}}\Vert f-f_z\Vert _{L^2(\Omega _z)}^2\bigr )^{1/2}. {\tag{4.5}} \end{multline}$$

In the infimum, “${v|_{\Gamma _D}=u_D}$” stands for all $v\in H^1(\Omega )$ with $v=u_D$ on $\Gamma _D$.

Proof.

Abbreviate $e:=u-u_h$ and let $q_h\in \mathcal{S}^k(\mathcal{T})^d$. Assume that $v\in H^1(\Omega )$ satisfies $v=u_D$ on $\Gamma _D$ and $\Vert \nabla (u_h-v)\Vert _{2}\le \Vert \nabla e\Vert _{2}$. Recall $p=\nabla u$ and $p_h=\nabla u_h$. Then 4.1-4.4 imply 3.1. Hence, we may choose $q=q_h$ and $w=u-v$ in Theorem 3.1 to obtain with Cauchy’s inequality for the second term that

$$\begin{equation} \begin{split} \Vert \nabla e\Vert _{2}^2&=\int _\Omega \nabla e\cdot \nabla w\,dx +\int _\Omega \nabla e\cdot \nabla (v-u_h)\,dx\\ &\le \Vert \nabla w\Vert _{2}\,\Bigl ( c_2\Vert p_h-q_h\Vert _{2} +c_5 \Vert h_\mathcal{E}^{1/2}\,(g-q_h\cdot n)\Vert _{2,\Gamma _N}\\ &\quad \hphantom {\Vert \nabla w\Vert _{2}\,\Bigl ( c_2} +c_3\bigl (\sum _{z\in {\mathcal{K}}} h_z^2 \min _{f_z\in \mathbb{R}} \Vert f+\operatorname {div}q_h-f_z\Vert _{2,\Omega _z}^2\bigr )^{1/2}\Bigr ) \\ &\quad +\Vert \nabla (u_h-v)\Vert _{2}\,\Vert \nabla e\Vert _{2}. \end{split} {\tag{4.6}} \end{equation}$$

Since $\Vert \nabla w\Vert _{2}\le \Vert \nabla e\Vert _{2} +\Vert \nabla (u_h-v)\Vert _{2}\le 2\Vert \nabla e\Vert _{2}$, we can divide 4.6 by $\Vert \nabla e\Vert _{2}$ to verify

$$\begin{multline} \Vert \nabla e\Vert _{2}\le 2c_2\Vert p_h-q_h\Vert _{2} +2c_5\Vert h_\mathcal{E}^{1/2}\,(g-q_h\cdot n)\Vert _{2,\Gamma _N} +\Vert \nabla (u_h-v)\Vert _{2}\\ +2c_3\bigl (\sum _{z\in {\mathcal{K}}} h_z^2 \min _{f_z\in \mathbb{R}} \Vert f+\operatorname {div}q_h-f_z\Vert _{2,\Omega _z}^2\bigr )^{1/2}. {\tag{4.7}} \end{multline}$$

Let $\operatorname {div}_\mathcal{T}$ denote the $\mathcal{T}$-piecewise action of the $\operatorname {div}$-operator. The triangle inequality in the last summand in 4.7 and $h_z\lesssim h_T$ for $z\in T\cap \mathcal{N}$ and $T\in \mathcal{T}$ and a summation over elements show

$$\begin{equation} \begin{split} &\sum _{z\in {\mathcal{K}}} h_z^2 \min _{f_z\in \mathbb{R}} \Vert f+\operatorname {div}q_h-f_z\Vert _{2,\Omega _z}^2\\ &\qquad \lesssim \Vert h_\mathcal{T}\,\operatorname {div}_\mathcal{T}(p_h-q_h)\Vert _{2}^2 + \sum _{z\in {\mathcal{K}}} h_z^2 \min _{f_z\in \mathbb{R}}\Vert f+\operatorname {div}_\mathcal{T}p_h-f_z\Vert _{2,\Omega _z}^2. \end{split} {\tag{4.8}} \end{equation}$$

Note that $\operatorname {div}_\mathcal{T}p_h=\Delta _\mathcal{T}u_h=0$ for our choices of $u_h\in \mathcal{S}^1(\mathcal{T})$. A $\mathcal{T}$-elementwise inverse estimate shows $\Vert h_\mathcal{T}\,\operatorname {div}_\mathcal{T}(p_h-q_h)\Vert _{2}\lesssim \Vert p_h-q_h\Vert _{2}$ (with a constant that depends on the shape of the finite elements only). Utilising this in 4.7–4.8, we deduce 4.5.

■

Remark 4.1.

In the proof of Theorem 4.1 we used the assumption that $u_h$ is of lowest order, i.e., $\nabla u_h \in \mathcal{L}^0(\mathcal{T})^d$, for the purpose of estimating $\|h_\mathcal{T}\operatorname {div}q_h \|_{L^2(\Omega )}$ by $\|q_h - \nabla u_h\|_{L^2(\Omega )}$. We refer to BC for related error estimates for higher order methods.

The subsequent lemma shows that $\inf _{v|_{\Gamma _D}=u_D} \Vert \nabla (u_h-v)\Vert _{L^2(\Omega )}$ is a higher order term.

Lemma 4.1.

Suppose that $u_h(z)=u_D(z)$ for all $z\in \mathcal{N}\cap \Gamma _D$. Then there exists an $h_\mathcal{E}$-independent constant $c_{10}>0$ (that depends on the shapes of the elements only) such that

$$\begin{equation} \inf _{v|_{\Gamma _D}=u_D} \Vert \nabla (u_h-v)\Vert _{L^2(\Omega )} \le c_{10} \Vert h_\mathcal{E}^{1/2}\partial (u_h-u_D)/\partial s\Vert _{L^2(\Gamma _D)}. {\tag{4.9}} \end{equation}$$

If $u_D \in H^2(\mathcal{E}_D):= \{ v\in L^2(\Gamma _D): \forall E\in \mathcal{E}_D,\; v|_E \in H^2(E) \}$, we have

$$\begin{equation} \inf _{v|_{\Gamma _D}=u_D} \Vert \nabla (u_h-v)\Vert _{L^2(\Omega )} \le c_{10} \Vert h_\mathcal{E}^{3/2}\partial ^2_{\mathcal{E}} u_D/\partial s^2\Vert _{L^2(\Gamma _D)}. {\tag{4.10}} \end{equation}$$

Proof.

Let $E\in \mathcal{E}_D$ belong to some $T\in \mathcal{T}$ and denote $\gamma :=\Gamma _D\cap \partial T$. We determine $w\in H^1(T)$ by extending the boundary values $w|_\gamma =u_h-u_D$ and $w|_{\partial T\setminus \gamma }=0$. Note that $w$ is continuous on $\partial T$ since $u_h$ interpolates $u_D=v$ at each node on $\Gamma _D$. An harmonic extension of $w|_{\partial T}$ to $w\in H^1(\mathcal{T})$ yields

$$\begin{equation} \Vert \nabla w\Vert _{2,T}\lesssim \Vert w\Vert _{H^{1/2}(\partial T)} \lesssim \Vert w\Vert _{2,\partial T}^{1/2}\Vert \partial w/\partial s\Vert _{2,\partial T}^{1/2}, {\tag{4.11}} \end{equation}$$

where we applied an interpolation estimate. A one-dimensional integration argument shows $\Vert w\Vert _{2,\partial T}\le h_T\,\Vert \partial w/\partial s\Vert _{2,\partial T}$. Consequently,

$$\begin{equation} \Vert \nabla w\Vert _{2,T}\lesssim h_T^{1/2}\, \Vert \partial w/\partial s\Vert _{2,\partial T} = h_T^{1/2}\,\Vert \partial (u_h-u_D)/\partial s\Vert _{2,\gamma }. {\tag{4.12}} \end{equation}$$

A scaling argument guarantees that the constant in 4.12 is $h_T$-independent. Defining $v$ by $u_h-w$ on elements on $\Gamma _D$ and by zero on other elements then shows the lemma. The second estimate follows from $\Vert w\Vert _{2,\partial T}\le h_T^{2}\,\Vert \partial ^2 w/\partial s^2\Vert _{2,\partial T}$.

■

Lemma 4.2.

Suppose $g\in H^1(\mathcal{E}_N)$ and, for each node $z\in \mathcal{N}\cap \overline{\Gamma }_N$ where the outer unit normal $n$ on $\Gamma _N$ is continuous (hence constant in a neighbourhood of $z$ as $\Gamma _N$ is a polygon), let $g$ be continuous. Then, the set

$$\begin{equation} \mathcal{S}_N^1(\mathcal{T},g):=\{q_h\in \mathcal{S}^1(\mathcal{T})^d:\forall E\in \mathcal{E}_N\; \forall z\in E\cap \mathcal{N},\; q_h(z)\cdot n_E=g(z)\} {\tag{4.13}} \end{equation}$$

is nonvoid and, for each $q_h \in \mathcal{S}_N^1(\mathcal{T},g)$,

$$\begin{equation} \Vert h_\mathcal{E}^{1/2}(g-q_h\cdot n)\Vert _{L^2(\Gamma _N)} \le \Vert h_\mathcal{E}^{3/2}\partial _\mathcal{E}g/\partial s\Vert _{L^2(\Gamma _N)}. {\tag{4.14}} \end{equation}$$

Proof.

Elementary estimates on each edge on $\Gamma _N$ verify 4.14; the proof of $\mathcal{S}_N^1(\mathcal{T},g) \neq \emptyset$ follows from an explicit construction in Example 4.1.

■

Example 4.1.

We define an operator ${\mathcal{A}}:L^2(\Omega )^2 \to {\mathcal{S}}_N^1(\mathcal{T},g)$ by

$$\begin{equation} {\mathcal{A}} p := \sum _{z\in {\mathcal{N}}} p_z \varphi _z, {\tag{4.15}} \end{equation}$$

where $p_z := -\kern -1.00em\int _{\omega _z} p \,dx := \frac{1}{|\omega _z|} \int _{\omega _z} p \,dx \in \mathbb{R}^2$ for $z\in {\mathcal{N}}\setminus \overline{\Gamma }_N$ while we incorporate ${\mathcal{A}} p(z) \cdot n_E = g(z)$ for $z \in \mathcal{N}\cap \overline{\Gamma }_N$. In case $z = E_1\cap E_2$ for two distinct edges $E_1,E_2 \in \mathcal{E}_N$ with distinct outer unit normals $n_{E_1}$, $n_{E_2}$ on $E_1$, $E_2$ at a corner $z$ we choose $p_z \in \mathbb{R}^2$ to be the unique solution of the $2\times 2$ linear system

$$\begin{equation} n_{E_1} \cdot p_z = g|_{E_1}(z) \quad \text{and}\quad n_{E_2} \cdot p_z = g|_{E_2}(z). {\tag{4.16a}} \end{equation}$$In the remaining cases $z\in E_1 \cap \Gamma _D$ for $E_1 \in \mathcal{E}_N$ or $z=E_1\cap E_2$ with two parallel edges $E_1,E_2 \in \mathcal{E}_N$ with the unit tangent vector $t_{E_1}$ let $p_z\in \mathbb{R}^2$ solve$$\begin{equation} n_{E_1} \cdot p_z = g|_{E_1}(z) \quad \text{and}\quad t_{E_1} \cdot p_z = -\kern -1.00em\int _{\omega _z} t_{E_1} \cdot p \,dx. {\tag{4.16b}} \end{equation}$$■

The following corollary is 1.5 with a constant $c_1=c_9$ as in Theorem 4.1 and with specified higher order terms from Lemma 4.1 and 4.2 and a Poincaré inequality.

Corollary 4.1.

Under the conditions of Theorem 4.1 and Lemmas 4.1 and 4.2 we have for $f\in H^1(\Omega )$ that

$$\begin{multline} \Vert \nabla (u-u_h)\Vert _{L^2(\Omega )}\le c_9\, \min _{q_h\in \mathcal{S}_N^1(\mathcal{T},g)}\Vert \nabla u_h-q_h\Vert _{L^2(\Omega )} \\ +c_{11} \Bigl (\Vert h_\mathcal{E}^{3/2}\,\partial _{\mathcal{E}}^2 u_D/\partial s^2\Vert _{L^2(\Gamma _D)} +\Vert h_\mathcal{E}^{3/2}\,\partial _{\mathcal{E}} g/\partial s\Vert _{L^2(\Gamma _N)} +\Vert h_\mathcal{T}^2\,\nabla f\Vert _{L^2(\Omega )}\Bigr ). {\tag{4.17}} \end{multline}$$

The $(h_\mathcal{T},h_\mathcal{E})$-independent constant $c_{11} >0$ depends on the shape of the elements and patches only.■

Remark 4.2.

Let us emphasise that the derivatives along $\Gamma$ are required only $\mathcal{E}$-piecewisely while $f$ needs to be patch-wise (not only elementwise) in $H^1$ and so $f\in H^1(\Omega )$. For a nonsmooth right-hand side $f$, $\Vert h_\mathcal{T}^2\,\nabla f\Vert _{L^2(\Omega )}$ may be replaced by a patch-wise $L^2$–best approximation error in the approximation through constants of $f$ (cf. 2.7).

The global averaging process might be too expensive or its approximation may be inefficient and hence a local averaging process of interest. Recall that $\omega _E$ is the (interior of the) union of all elements in $\mathcal{T}$ that share the edge $E\in \mathcal{E}$.

Corollary 4.2.

Under the conditions of Theorem 4.1 and Lemmas 4.1 and 4.2 we have for $f\in H^1(\Omega )$ that

$$\begin{equation} \begin{split} \Vert \nabla (u-u_h)\Vert _{L^2(\Omega )}&\le c_{12}\, \Bigl (\sum _{E\in \mathcal{E}}\min _{q_E\in \mathcal{S}^1(\mathcal{T}_E)^d} \Bigl (\| \nabla u_h-q_E\|_{L^2(\omega _E)}^2\\ &\hphantom {\le c_{12}\, \Bigl (\sum _{E\in \mathcal{E}}\min _{q_E\in \mathcal{S}^1(\mathcal{T}_E)^d} \Bigl (\|} +h_E\|g_h-q_E\cdot n\|_{L^2(E\cap \Gamma _N)}^2\Bigr )\Bigr )^{1/2} \\ &\quad +c_{11}\Bigl (\Vert h_\mathcal{E}^{3/2}\,\partial _{\mathcal{E}}^2 u_D/\partial s^2\Vert _{L^2(\Gamma _D)}\\ &\quad \hphantom {+c_{11}\Bigl (\|} +\Vert h_\mathcal{E}^{3/2}\,\partial _{\mathcal{E}} g/\partial s\Vert _{L^2(\Gamma _N)} +\Vert h_\mathcal{T}^2\,\nabla f\Vert _{L^2(\Omega )}\Bigr ). \end{split} {\tag{4.18}} \end{equation}$$

The $(h_\mathcal{T},h_\mathcal{E})$-independent constant $c_{12}=\max \{c_9,2c_5\}/c_7$ depends on the shape of the elements and patches only.

Proof.

Theorem 4.1, Lemma 4.1, an approximation $g_h$ of $g$ as in Lemma 4.2, and a Poincaré inequality show

$$\begin{equation} \begin{split} &\Vert \nabla (u-u_h)\Vert _{2}\\ &\qquad \lesssim \min _{q_h\in \mathcal{S}^1(\mathcal{T})^d} (\Vert \nabla u_h-q_h\Vert _{2} + \Vert h_\mathcal{E}^{1/2}(g_h-q_h \cdot n)\Vert _{2,\Gamma _N}) +\Vert h_\mathcal{E}^{3/2}\partial _{\mathcal{E}} g/\partial s\Vert _{2,\Gamma _N}\\ &\qquad \hphantom {\lesssim \min _{q_h\in \mathcal{S}^1(\mathcal{T})^d}} +\Vert h_\mathcal{E}^{3/2}\partial _{\mathcal{E}}^2 u_D/\partial s^2\Vert _{2,\Gamma _D} +\Vert h_\mathcal{T}^2 \nabla f\Vert _{2}. \end{split} {\tag{4.19}} \end{equation}$$

This and the first inequality of Theorem 3.2 imply the assertion.

■

Remark 4.3.

The results of this section hold also in three dimensions where $\mathcal{T}$ consists of tetrahedra or parallelepipeds. The proofs of some details as Lemma 4.1 or Lemma 4.2 require much more technical preparations and so are omitted in this overview.

Remark 4.4.

It is shown in CVC2 that the edge-contributions (jump differences in the normal fluxes components across edges) dominate in standard residual a posteriori error estimates BaRBBSCF1EEHJV. Arguing as in R1R2DMR, one can hence derive alternative proofs of 4.18 and then of 4.17.

Remark 4.5.

In an $L^\infty$-estimate of HSWW it is suggested to average over a domain of size $O(h \log (1/h))$ instead of merely over patches or the entire domain to obtain asymptotic exact results.

5. Applications to nonconforming finite element schemes

In the Laplace problem with mixed boundary conditions 4.1–4.3, we suppose that the discrete flux $p_h:=\nabla _\mathcal{T}u_h\in \mathcal{L}^0(\mathcal{T})^d$, where $\nabla _\mathcal{T}$ denotes the $\mathcal{T}$-piecewise application of the gradient, satisfies

$$\begin{equation} \int _\Omega \nabla _\mathcal{T}u_h\cdot \nabla w_h\, dx= \int _\Omega f w_h\,dx+\int _{\Gamma _N} g w_h\,ds \quad \text{for all }w_h\in \mathcal{S}^1_D(\mathcal{T}). {\tag{5.1}} \end{equation}$$

The usual conformity conditions read for all $E\in \mathcal{E}_\Omega \cup \mathcal{E}_D$,

$$\begin{equation} \int _E [u_h]\, ds=0, {\tag{5.2}} \end{equation}$$

where $[u_h]|_E$ denotes the jump of $u_h$ across $E\in \mathcal{E}_\Omega$ and denotes $u_D-u_h$ on $\Gamma _D$. Those conditions are satisfied by construction for Crouzeix–Raviart finite elements of lowest order.

Remark 5.1.

It is stressed that $\mathcal{S}^1_D(\mathcal{T})$ is a conforming test function space which is included in the nonconforming finite element spaces for triangles or tetrahedra. For parallelograms, 5.1 means that the polynomial degrees are at least of second order to include the conforming term $x_1\, x_2$. This technical detail could actually be dropped since the contribution from an enhanced finite element space leads to a higher order term KS. We restrict our analysis to triangles or tetrahedra for simplicity.

Theorem 5.1.

Suppose that $\Gamma _N$ is connected and that $\Gamma _D$ belongs to only one connectivity component of $\partial \Omega$. Then, there exists an $(h_\mathcal{T},h_\mathcal{E})$-independent constant $c_{13}>0$ (that depends on $k\ge 1$ and the shape of the elements and patches) such that

$$\begin{equation} \begin{split} &\Vert \nabla _\mathcal{T}(u-u_h)\Vert _{L^2(\Omega )}\\ &\qquad \le \min _{q_h\in \mathcal{S}^k(\mathcal{T})^d}\bigl (c_{13}\Vert \nabla _\mathcal{T}u_h-q_h\Vert _{L^2(\Omega )} +c_5\Vert h_\mathcal{E}^{1/2}(g-q_h\cdot n)\Vert _{L^2(\Gamma _N)}\\ &\hphantom {\qquad \le \min _{q_h\in \mathcal{S}^k(\mathcal{T})^d}\bigl (c_{13}\Vert \nabla _\mathcal{T}u_h-q_h\Vert _{L^2(\Omega )} } \quad +c_5 \Vert h_\mathcal{E}^{1/2}(q_h\cdot t-\partial u_D/\partial s)\Vert _{L^2(\Gamma _D)}\bigr )\\ &\qquad \quad +c_3\bigl (\sum _{z\in {\mathcal{K}}} h_z^2 \min _{f_z\in \mathbb{R}}\Vert f-f_z\Vert _{L^2(\Omega _z)}^2\bigr )^{1/2}. \end{split} {\tag{5.3}} \end{equation}$$

Here, $t \in \mathcal{L}^0(\mathcal{E}_D)^d$ denotes the unit tangent vector on $\Gamma _D$.

Remark 5.2.

The following lemma is based on the Helmholtz decomposition of a vector field. The decomposition is available in three dimensions as well (cf., e.g., [GR]) but the notation is more involved so we restrict the discussion to the two-dimensional setting for brevity.

Lemma 5.1.

For all $p-p_h\in L^2(\Omega )^2$, there exist $\alpha ,\beta \in H^1(\Omega )$ that satisfy boundary conditions $\alpha |_{\Gamma _D}=0$ and $\beta |_{\Gamma _N}$ is constant such that

$$\begin{equation} p-p_h=\nabla \alpha +\operatorname {Curl}\beta \quad \text{and}\quad \Vert p-p_h\Vert _{L^2(\Omega )}^2=\Vert \nabla \alpha \Vert _{L^2(\Omega )}^2+ \Vert \nabla \beta \Vert _{L^2(\Omega )}^2. {\tag{5.4}} \end{equation}$$

Proof.

The lemma follows from the Helmholtz decomposition where $\alpha \in H^1_D(\Omega )$ solves $\Delta \alpha =\operatorname {div}(p-p_h)$ and $\Delta \beta =\operatorname {curl}(p-p_h)$ with proper boundary conditions (cf., e.g., GR).

■

Proof of Theorem 5.1.

For $p=\nabla u$ and $p_h=\nabla _\mathcal{T}u_h$, Lemma 5.1 yields

$$\begin{equation} \Vert p-p_h\Vert _{2}^2=\int _\Omega (p-p_h)\cdot \nabla \alpha \,dx +\int _\Omega (p-p_h)\cdot \operatorname {Curl}\beta \,dx. {\tag{5.5}} \end{equation}$$

Since $\Gamma _N$ is connected, we may and will assume without loss of generality that $\beta =0$ on $\Gamma _N$. According to 4.1–4.3 and 5.1, we infer 3.1 and hence may choose $q=q_h\in \mathcal{S}^k(\mathcal{T})^d$ and, in case $\alpha \not \equiv 0$, $w=\alpha /\|\alpha \|_{2,\Omega }$ in Theorem 3.1 to obtain

$$\begin{multline} \int _\Omega (p-p_h)\cdot \nabla \alpha \,dx\le \Vert \nabla \alpha \Vert _{2}\,\Bigl ( c_2\Vert p_h-q_h\Vert _{2}+c_5 \Vert h_\mathcal{E}^{1/2}\,(g-q_h\cdot n)\Vert _{2,\Gamma _N}\\ +c_3\bigl (\sum _{z\in {\mathcal{K}}} h_z^2 \min _{f_z\in \mathbb{R}} \Vert f+\operatorname {div}q_h-f_z\Vert _{2,\Omega _z}^2\bigr )^{1/2}\Bigr ). {\tag{5.6}} \end{multline}$$

The estimate of the last term in 5.5 will follow from Theorem 3.1 as well once we establish an analogy to 3.1, namely

$$\begin{equation} \int _\Omega (p-p_h)\cdot \operatorname {Curl}\ w_h\, dx=0 \quad \text{for all }w_h\in \mathcal{S}_N^1(\mathcal{T}), {\tag{5.7}} \end{equation}$$

where $\mathcal{S}^1_N(\mathcal{T}):=\{v_h\in \mathcal{S}^1(\mathcal{T}): v_h=0\text{ on }\Gamma _N\}$. It is essential to notice that $\partial w_h/\partial s$ is constant and $[u_h]$ has a vanishing integral on any edge. An elementwise integration by parts on the left-hand side of 5.7 yields volume terms $(u-u_h)\operatorname {div}_\mathcal{T}\operatorname {Curl}\, w_h=0$ and edge terms $[(u-u_h)\partial w_h/\partial s]=[u_h]\partial w_h/\partial s$ whose integral vanishes on any $E$ (the case $E\in \mathcal{E}_\Omega$ is indicated and the assertion is true for $E\in \mathcal{E}_D$ as well; $w_h=0$ on $\Gamma _N$ shows it for $E\in \mathcal{E}_N$). In this way we establish 5.7.

To employ Theorem 3.1, we interchange components, writing in this proof $Q(a_1,a_2):=(-a_2,a_1)$ for vectors, and we interchange the role of the boundaries and adopt Theorem 2.1 and 3.1 where $\tilde{\Gamma }_D=\Gamma _N$ acts as the Dirichlet boundary and $\tilde{\Gamma }_N=\Gamma _D$ acts as the Neumann boundary. Writing $\tilde{p}=Qp$ and $\tilde{p}_h=Q p_h$, 5.7 reads $\int _\Omega (\tilde{p}-\tilde{p}_h)\cdot \nabla w_h\, dx=0$ for all $w_h \in \tilde{\mathcal{S}}^1_D(\mathcal{T})=\mathcal{S}^1_N(\mathcal{T})$ and this is 3.1. Reading Theorem 3.1 in the present notation, we obtain

$$\begin{equation} \begin{split} &\int _\Omega (p-p_h)\cdot \operatorname {Curl}\beta \,dx\\ &\qquad \le \Vert \nabla \beta \Vert _{2}\,\Bigl ( c_2\Vert p_h-q_h\Vert _{2}+c_5 \Vert h_\mathcal{E}^{1/2}\,(\partial u_D/\partial s-q_h\cdot t)\Vert _{2,\Gamma _D}\\ &\qquad \hphantom {\le \Vert \nabla \beta \Vert _{2}\,\Bigl ( c_2\Vert p_h-q_h\Vert _{2}}\, +c_3\bigl (\sum _{z\in {\mathcal{K}}} h_z^2 \Vert \operatorname {curl}_\mathcal{T}(q_h-p_h)\Vert _{2,\Omega _z}^2\bigr )^{1/2}\Bigr ) \end{split} {\tag{5.8}} \end{equation}$$

with $\operatorname {curl}_\mathcal{T}(q_h-p_h) := \operatorname {div}_\mathcal{T}Q(q_h-p_h)$. In the second last term, $t=Qn$ denotes the unit tangent vector and in the last term we used that $\operatorname {curl}_\mathcal{T}p=0=\operatorname {curl}_\mathcal{T}p_h$.

The remaining arguments are similar to those in the proof of Theorem 4.1 and hence are omitted.

■

In contrast to the conforming situation, Theorem 5.1 demands averaging functions to satisfy some conditions on the Dirichlet boundary.

Lemma 5.2.

Suppose $u_D\in H^2(\mathcal{E}_D)$ and, for each node $z\in \mathcal{N}\cap \Gamma _D$ where the outer unit normal $n$ on $\Gamma$ is continuous, let $\partial u_D/\partial s$ be continuous. Then, the set

$$\begin{equation} \mathcal{S}^1_D(\mathcal{T},u_D):=\{ q_h\in \mathcal{S}^1(\mathcal{T})^d:\, \forall E\in \mathcal{E}_D \, \forall z\in E\cap \mathcal{N}, \, q_h(z)\cdot t_E = \partial u_D/\partial s(z)\} {\tag{5.9}} \end{equation}$$

is nonvoid and, for each $q_h\in \mathcal{S}^1_D(\mathcal{T},u_D)$,

$$\begin{equation} \|h_\mathcal{E}^{1/2}(q_h\cdot t_\mathcal{E}- \partial u_D / \partial s) \|_{L^2(\Gamma _D)} \le \|h_\mathcal{E}^{3/2} \partial _\mathcal{E}^2 u_D/\partial s^2 \|_{L^2(\Gamma _D)}. {\tag{5.10}} \end{equation}$$

Proof.

Similar to 4.10 in Lemma 4.1 or 4.2.

■

Example 5.1.

Assume the conditions of Lemmas 4.2 and 5.2 on the data $g$ and $u_D$. We define an operator $\mathcal{A}:L^2(\Omega )^2 \to \mathcal{S}^1_N(\mathcal{T},g) \cap \mathcal{S}^1_D(\mathcal{T},u_D)$ by 4.15 and $p_z:=-\kern -1.00em\int _{\omega _z}p\,dx$ for $z\in \mathcal{N}\setminus \Gamma$. In case $z\in (\mathcal{N}\cap \overline{\Gamma }_N)\setminus \Gamma _D$ we preceed as in 4.16a, resp. 4.16b. In case $z\in (\mathcal{N}\cap \Gamma _D)\setminus \overline{\Gamma }_N$ we consider the analogous $2\times 2$ systems

$$\begin{equation} t_{E_1}\cdot p_z = \partial u_D |_{E_1}/\partial s (z) \quad \text{and}\quad t_{E_2}\cdot p_z = \partial u_D |_{E_2}/\partial s (z) {\tag{5.11a}} \end{equation}$$(cf. 4.16a and notation from Example 4.1), resp., as an analog to 4.16b,$$\begin{equation} t_{E_1}\cdot p_z = \partial u_D |_{E_1}/\partial s (z) \quad \text{and}\quad n_{E_1}\cdot p_z = -\kern -1.00em\int _{\omega _z} n_{E_1} \cdot p \,dx. {\tag{5.11b}} \end{equation}$$

For $z \in \Gamma _D\cap \overline{\Gamma }_N \subseteq \mathcal{N}$ we require a compatibility condition if $n_{E_1}=t_{E_2}$, namely $g(z)= \partial u_D/ \partial s (z)$. Then, we define $p_z\in \mathbb{R}^2$ as in 5.11b when $E_1\subseteq \Gamma _D$ (the case $E_2 \subseteq \Gamma _D$ is analogous). For $n_{E_1}\neq \pm t_{E_2}$, we need no further compatibility of the data and solve the $2\times 2$ linear system

$$\begin{equation} t_{E_1} \cdot p_z = \partial u_D|_{E_1}/\partial s (z) \quad \text{and}\quad n_{E_2}(z)\cdot p_z = g(z). {\tag{5.12}} \end{equation}$$

(Here $E_1\subseteq \Gamma _D$; the case $E_2\subseteq \Gamma _D$ is analogous.)■

The modification of 1.5 in the nonconforming setting is a direct consequence of Theorem 5.1, Lemma 5.2 and Example 5.1. Note that Corollary 4.1 is a special case apart from the different treatment of the Dirichlet boundary conditions.

Corollary 5.1.

Under the conditions of Theorem 5.1 and Lemmas 4.2 and 5.2, we have for $f\in H^1(\Omega )$ that

$$\begin{multline} \|\nabla _\mathcal{T}(u-u_h)\|_{L^2(\Omega )} \le c_{13} \min _{q_h\in \mathcal{S}^1_N(\mathcal{T},g)\cap \mathcal{S}^1_D(\mathcal{T},u_D)} \|\nabla _\mathcal{T}(u_h-q_h)\|_{L^2(\Omega )} \\ +c_{14} \bigl ( \|h_\mathcal{E}^{3/2}\partial _\mathcal{E}^2 u_D/\partial s^2\|_{L^2(\Gamma _D)} + \|h_\mathcal{E}^{3/2}\partial _\mathcal{E}g/\partial s\|_{L^2(\Gamma _N)} + \|h_\mathcal{T}^2\nabla f\|_{L^2(\Omega )} \bigr ). {\tag{5.13}} \end{multline}$$■

The analog to Corollary 4.2 concludes this section on lowest order Raviart–Crouzeix finite elements.

Corollary 5.2.

Under the conditions of Theorem 5.1 and Lemmas 4.2 and 5.2, there exists a constant $c_{15}>0$ such that we have for $f\in H^1(\Omega )$

$$\begin{equation} \begin{split} &\|\nabla _\mathcal{T}(u-u_h)\|_{L^2(\Omega )} \\ &\qquad \le c_{15} \Bigl ( \sum _{E\in \mathcal{E}} \min _{q_E\in \mathcal{S}^1(\mathcal{T}_E)^d} \bigl ( \|\nabla _\mathcal{T}u_h - q_E \|_{L^2(\omega _E)}^2 \\ &\qquad \hphantom {\le c_{15} \Bigl ( \sum _{E\in \mathcal{E}} } + h_E \|g_h-q_E\cdot n\|_{L^2(E\cap \Gamma _N)}^2 + h_E \|q_E\cdot t - u_{D,h}'\|_{L^2(E\cap \Gamma _D)}^2 \bigr )\Bigr )^{1/2} \\ &\qquad \quad +c_{14} \bigl ( \|h_\mathcal{E}^{3/2}\partial _\mathcal{E}^2 u_D/\partial s^2\|_{L^2(\Gamma _D)} + \|h_\mathcal{E}^{3/2}\partial _\mathcal{E}g/\partial s\|_{L^2(\Gamma _N)} +\|h_\mathcal{T}^2\nabla f\|_{L^2(\Omega )} \bigr ). \end{split} {\tag{5.14}} \end{equation}$$

Here, $u_{D,h}'$ denotes an approximation of $\partial u_D/\partial s$ as in Lemma 5.2, i.e., $u_{D,h}'= q_h \cdot t_\mathcal{E}$ on $\Gamma _D$ for some $q_h \in \mathcal{S}^1_D(\mathcal{T},u_D)$.■

Remark 5.3.

The results of this section can be generalised to three space dimensions as all the required tools such as a Helmholtz decomposition are available then as well. Details on the three-dimensional case are omitted for notational simplicity.

Remark 5.4.

Arguing as in CVC2, one can prove that edge contributions (jumps in the fluxes across edges) dominate the residual based error estimates from DDPVC2KS. Arguing in the spirit of R1R2DMR, one can hence derive alternative proofs of 5.14 and then of 5.13.

6. Applications to mixed finite element schemes

In the Laplace problem with mixed boundary conditions 4.1–4.3, we suppose that the discrete flux $p_h\in H(\operatorname {div},\Omega )\cap \mathcal{L}^k(\mathcal{T})^2$ and the displacement approximation $u_h\in \mathcal{L}^k(\mathcal{T})$ satisfy, for all $q_h\in \operatorname {Curl}\mathcal{S}^1(\mathcal{T})$ with $q_h\cdot n=0$ on $\Gamma _N$, and for all $T\in \mathcal{T}$ and $E\in \mathcal{E}_N$ that

$$\begin{eqnarray} &\int _\Omega (p_h\cdot q_h+u_h\,\operatorname {div}q_h)dx =\int _{\Gamma _D} u_D\,q_h\cdot n\,ds,&{\tag{6.1}}\\ &\int _T (f+\operatorname {div}p_h) dx =0,& {\tag{6.2}}\\ &\int _E (g-p_h\cdot n)ds =0.& {\tag{6.3}} \end{eqnarray}$$

Remark 6.1.

Standard mixed finite element methods of any order, such as Raviart–Thomas (RT), Brezzi–Douglas–Marini (BDM), or Brezzi–Douglas–Fortin–Marini (BDFM) elements (cf. BF for details), provide 6.1–6.3 C1.

Theorem 6.1.

Suppose that $\Gamma _N$ is connected and that $\Gamma _D$ belongs to only one connectivity component of $\partial \Omega$ and let $f\in H^1(\mathcal{T})$, i.e., $f|_T\in H^1(T)$ for all $T\in \mathcal{T}$. Then, there exists an $(h_\mathcal{T},h_\mathcal{E})$-independent constant $c_{16}>0$ (that depends on $k\ge 1$ and the shape of the elements and patches) such that

$$\begin{equation} \begin{split} \Vert p-p_h\Vert _{L^2(\Omega )}&\le \min _{q_h\in \mathcal{S}^k(\mathcal{T})^2}\Bigl ( c_2\Vert p_h-q_h\Vert _{L^2(\Omega )}+c_5 \Vert h_\mathcal{E}^{1/2}\,(\partial u_D/\partial s-q_h\cdot t)\Vert _{L^2(\Gamma _D)}\\ &\hphantom {\le \min _{q_h\in \mathcal{S}^k(\mathcal{T})^2}\Bigl ( c_2\|p_h-} +c_3\bigl (\sum _{z\in {\mathcal{K}}} h_z^2 \min _{f_z\in \mathbb{R}} \Vert \operatorname {curl}_\mathcal{T}q_h-f_z\Vert _{L^2(\Omega _z)}^2\bigr )^{1/2}\Bigr )\\ &\quad +c_{16}\Vert h_\mathcal{T}^2\nabla _\mathcal{T}(f+\operatorname {div}p_h)\Vert _{2,\Omega } +c_{16}\Vert h^{1/2}_\mathcal{E}(g-p_h\cdot n)\Vert _{2,\Gamma _N}. \end{split} {\tag{6.4}} \end{equation}$$

Proof.

Lemma 5.1 provides 5.5, and we may and will assume without loss of generality that $\beta =0$ on $\Gamma _N$. An integration by parts and 6.2–6.3 show for the $\mathcal{T}$-piecewise integral mean $\alpha _\mathcal{T}\in \mathcal{L}^0(\mathcal{T})$ of $\alpha \in H^1_D(\Omega )$ that

$$\begin{multline} \int _\Omega (p-p_h)\cdot \nabla \alpha \,dx =\int _\Omega (f+\operatorname {div}p_h)\,\alpha \,dx +\int _{\Gamma _N}(g-p_h\cdot n)\,\alpha \,dx\\ =\int _\Omega (f+\operatorname {div}p_h)\,(\alpha -\alpha _\mathcal{T})\,dx +\int _{\Gamma _N}(g-p_h\cdot n)\,(\alpha -\alpha _\mathcal{T})\,dx. {\tag{6.5}} \end{multline}$$

The second last term is estimated with an elementwise Poincaré inequality while the last term in 6.5 involves a trace theorem BSCF1Cl, namely

$$\begin{equation} \Vert a\Vert _{2,E}\lesssim h_E^{-1/2}\,\Vert a\Vert _{2,T_E} + h_E^{1/2}\,\Vert \nabla a\Vert _{2,T_E} {\tag{6.6}} \end{equation}$$

for $a=\alpha -\alpha _\mathcal{T}\in H^1(T_E)$ on the triangle $T_E\in \mathcal{T}$ and the edge $E\in \mathcal{E}_N$, $E\subset \partial T_E$. With a second application of Poincaré’s inequality, 6.6, and Cauchy’s inequality we show that

$$\begin{equation} \begin{split} \int _\Omega (p-p_h)\cdot \nabla \alpha \,dx&\lesssim \Vert h_\mathcal{T}(f+\operatorname {div}p_h)\Vert _{2}\,\Vert \nabla \alpha \Vert _{2}\\ &\qquad +\Vert h^{1/2}_\mathcal{E}(g-p_h\cdot n)\Vert _{2,\Gamma _N} \Vert h^{-1/2}_\mathcal{E}(\alpha -\alpha _\mathcal{T})\Vert _{2,\Gamma _N} \\ &\lesssim \Vert \nabla \alpha \Vert _{2}\Bigl (\Vert h_\mathcal{T}^2\nabla _\mathcal{T}(f+\operatorname {div}p_h)\Vert _{2}+\Vert h^{1/2}_\mathcal{E}(g-p_h\cdot n)\Vert _{2,\Gamma _N}\Bigr ). \end{split} {\tag{6.7}} \end{equation}$$

The second contribution on the right-hand side of 5.5 is analysed with Theorem 3.1, where, as in the proof of Theorem 5.1, we interchange components and the role of the boundary conditions. As already employed in C1C2, $\operatorname {Curl}w_h\in H(\operatorname {div},\Omega )$ for all $w_h\in \mathcal{S}^1_N(\mathcal{T})$. Moreover, $\operatorname {Curl}w_h\cdot n=\partial w_h/\partial s=0$ on $\Gamma _N$. Hence, 6.1 and an integration by parts for $p$ yield 5.7 because of 4.2. Arguing as in the proof of Theorem 5.1, we deduce for arbitrary $q_h \in \mathcal{S}^k(\mathcal{T})^2$ that

$$\begin{multline} \int _\Omega (p-p_h)\cdot \operatorname {Curl}\beta \,dx\le \Vert \nabla \beta \Vert _{2}\,\Bigl ( c_2\Vert p_h-q_h\Vert _{2}+c_5 \Vert h_\mathcal{E}^{1/2}\,(\partial u_D/\partial s-q_h\cdot t)\Vert _{2,\Gamma _D}\\ +c_3\bigl (\sum _{z\in {\mathcal{K}}} h_z^2 \min _{f_z\in \mathbb{R}} \Vert \operatorname {curl}_\mathcal{T}q_h-f_z\Vert _{2,\Omega _z}^2\bigr )^{1/2}\Bigr ). {\tag{6.8}} \end{multline}$$

The remaining details are analogous to the proof of Theorem 5.1 and hence are omitted.

■

The precise version of 1.5 for lowest order mixed finite element methods is summarised as follows.

Corollary 6.1.

Suppose that the discrete flux $p_h$ satisfies $\operatorname {curl}_\mathcal{T}p_h =0$, $\operatorname {div}_\mathcal{T}p_h \in \mathcal{L}^0(\mathcal{T})$ and $p_h \cdot n \in \mathcal{L}^0(\mathcal{E}_N)$. Then,

$$\begin{multline} \|p-p_h\|_{L^2(\Omega )} \le \min _{q_h\in \mathcal{S}^1_D(\mathcal{T},u_D)} c_{17}\|p_h-q_h\|_{L^2(\Omega )} \\ +c_{18} \bigl (\|h_\mathcal{E}^{3/2}\partial _\mathcal{E}^2 u_D/\partial s^2 \|_{L^2(\Gamma _D)} +\|h_\mathcal{E}^{3/2}\partial _\mathcal{E}g/\partial s\|_{L^2(\Gamma _N)} +\|h_\mathcal{T}^2 \nabla _\mathcal{T}f\|_{L^2(\Omega )} \bigr ). {\tag{6.9}} \end{multline}$$

Proof

Combine Theorem 6.1 and Lemmas 4.1 and 5.2, and use an inverse estimate to prove

$$\begin{align*} \sum _{z\in \mathcal{K}} h_z^2\|\operatorname {curl}_\mathcal{T}q_h \|_{2,\Omega _z}^2 &= \sum _{z\in \mathcal{K}}h_z^2 \|\operatorname {curl}_\mathcal{T}(q_h-p_h)\|_{2,\Omega _z}^2 \\ &\lesssim \sum _{z\in \mathcal{K}}\|q_h-p_h\|_{2,\Omega _z}^2 \lesssim \|q_h-p_h\|_2. \end{align*}$$■

Remark 6.2.

The assumptions in Corollary 6.1 are satisfied for lowest order Raviart–Thomas and Brezzi–Douglas–Fortin–Marini finite elements.

Example 6.1.

Assume the conditions of Lemma 5.2 on the data $u_D$. We define an operator $\mathcal{A}:L^2(\Omega )^2 \to \mathcal{S}^1_D(\mathcal{T},u_D)$ by 4.15 and $p_z:=-\kern -1.00em\int _{\omega _z}p\,dx$ for $z\in \mathcal{N}\setminus \Gamma _D$. In case $z\in \mathcal{N}\cap \Gamma _D$ we consider $2\times 2$ systems

$$\begin{equation} t_{E_1}\cdot p_z = \partial u_D |_{E_1}/\partial s (z) \quad \text{and}\quad t_{E_2}\cdot p_z = \partial u_D |_{E_2}/\partial s (z) {\tag{6.10a}} \end{equation}$$(cf. 4.16a and notation from Example 4.1), resp., as an analog to 4.16b,$$\begin{equation} t_{E_1}\cdot p_z = \partial u_D |_{E_1}/\partial s (z) \quad \text{and}\quad n_{E_1}\cdot p_z = -\kern -1.00em\int _{\omega _z} n_{E_1} \cdot p \,dx. {\tag{6.10b}} \end{equation}$$■

A local version follows from Theorem 3.2 and concludes this section on mixed finite element methods.

Corollary 6.2.

Under the conditions of Theorem 6.1 and Corollary 6.1 we have

$$\begin{multline*} \|p-p_h\|_{L^2(\Omega )} \le c_{19} \bigl ( \sum _{E\in \mathcal{E}} \min _{q_E\in \mathcal{S}^1(\mathcal{T}_E)^2} \|p_h-q_E\|_{L^2(\omega _E)}^2 +h_E\|u_{D,h}' - q_E\cdot t\|^2_{L^2(\Gamma _D)}\bigr )^{1/2} \\ + c_{18} \bigl (\|h_\mathcal{E}^{3/2}\partial _\mathcal{E}^2 u_D/\partial s^2 \|_{L^2(\Gamma _D)} +\|h_\mathcal{E}^{3/2}\partial _\mathcal{E}g /\partial s\|_{L^2(\Gamma _N)} +\|h_\mathcal{T}^2\nabla _\mathcal{T}f\|_{L^2(\Omega )} \bigr ). \end{multline*}$$■

Remark 6.3.

The results of this section could be generalised to three space dimensions. Details are omitted for brevity.

Remark 6.4.

For related residual based a posteriori error estimates we refer to ABVC1C2HW.

7. Numerical experiments

The theoretical results of this paper are supported by numerical experiments. In this section, we report on two examples of the problem 4.1–4.3 on uniform, adapted, and perturbed meshes for conforming, nonconforming, and mixed finite element methods.

Example 7.1.

Let $f:=0$ on the L-shaped domain $\Omega :=(-1,1)^2\setminus [0,1]\times [-1,0]$, $u_D:=0$ on the Dirichlet boundary $\Gamma _D:=\{0\}\times [-1,0]\cup [0,1] \times \{0\}$, and on the Neumann boundary $\Gamma _N:=\partial \Omega \setminus \Gamma _D$,

$$\begin{equation*} g(r,\varphi ):=2/3 \, r^{-1/3}(-\sin (\varphi /3), \cos (\varphi /3)) \cdot n \end{equation*}$$

using polar coordinates $(r,\varphi )$. The exact solution $u(r,\varphi ):=r^{2/3}\sin (2\varphi /3)$ of 4.1–4.3 has a typical corner singularity at the origin. In this example, the right-hand sides are smooth, but the solution is not. The coarsest triangulation $\mathcal{T}_0$ consists of three squares halved by diagonals parallel to the vector $(1,1)$ (cf. Figure 1).

Example 7.2.

Let $f:=-\Delta u$ for the function

$$\begin{equation*} u(x,y):=x(1-x)y(1-y) \arctan (60(r-1)), \quad r^2:=(x-1.25)^2+(y+0.25)^2 \end{equation*}$$

on the unit square $\Omega := (0,1)^2$ and set $u_D:=0$ on the entire boundary $\Gamma _D:=\partial \Omega$ ($\Gamma _N=\emptyset$). The solution $u$ to 4.1–4.3 is $H^2$-regular but $f$ (although theoretically smooth) has huge gradients on the circle with radius $1$ around $(1.25,-0.25)$. The coarsest triangulation $\mathcal{T}_0$ consists of four squares halved by diagonals parallel to the vector $(1,1)$ (cf. Figure 3).

The following adaptive algorithm generates all the sequences of meshes $\mathcal{T}_0,\mathcal{T}_1,\mathcal{T}_2,...$ in this paper which are uniform for $\Theta =0$ or adapted for $\Theta =1/2$ in 7.2. Since the resulting meshes might show local symmetries, we considered meshes that are either unperturbed (relative to $\mathcal{T}_0$) for $\vartheta =0$ and randomly perturbed for $\vartheta =1$ in step (e). The implementation was performed in Matlab in the spirit of ACF with a direct solution of linear systems of equations. For details on the red-blue-green refinements we refer to V.

Algorithm $(A^\vartheta _\Theta )$

(a)

Start with a coarse mesh $\mathcal{T}_0$, $k=0$.

(b)

Compute the discrete solution $p_h$ on the actual mesh $\mathcal{T}_k$.

(c)

Compute error indicators$$\begin{equation} \eta _{Z,T}:=\|p_h-\mathcal{A}p_h\|_{L^2(T)} {\tag{7.1}} \end{equation}$$

for all $T\in \mathcal{T}_k$ and plot energy error $e_N := \|p - p_h\|_{L^2(\Omega )}$ and its estimator $\eta _N^2 := \sum _{T\in \mathcal{T}} \eta _{Z,T}^2$ versus the degree of freedom $N$ of the triangulation $\mathcal{T}_k$.

(d)

Mark the element $T$ for red-refinement provided$$\begin{equation} \eta _{Z,T} \ge \Theta \,\max _{T'\in \mathcal{T}_k} \eta _{Z,T'}. {\tag{7.2}} \end{equation}$$

(e)

Mark further elements (red-blue-green refinement) to avoid hanging nodes. Generate a new triangulation $\tilde{\mathcal{T}}_{k+1}$ using edge-midpoints if $\vartheta =0$ and points on the edges at a random distance at most $0.3\,h_E$ from the edge-midpoints if $\vartheta =1$. Perturbe the nodes $z\in \mathcal{N}_{k+1}$ of the mesh $\tilde{\mathcal{T}}_{k+1}$ at random with values taken uniformly from a ball around $z$ of radius $\vartheta \,2^{-k}/15$. Correct boundary nodes by orthogonal projection onto that boundary piece they are expected such that $\Omega ,\Gamma _D,\Gamma _N$ are matched by the resulting mesh $\mathcal{T}_{k+1}$ exactly. Update $k$ and go to (b).

7.1. Results for conforming finite element methods

In the conforming finite element scheme, we use operator $\mathcal{A}$ from Example 4.1 in 7.1 of Algorithm $(A^\vartheta _\Theta )$ and report on results obtained for $(\Theta ,\vartheta )=(0,0)$ (uniform), $(1/2,0)$ (adaptive), and $(1/2,1)$ (adaptive, perturbed).

Some meshes obtained for Example 7.1 are shown in Figure 1 and illustrate a high automatic mesh-refinement of the adapted meshes towards the origin, which is expected to improve the convergence rate of $2/3$ possibly to the optimal value $1$. The result of the perturbation in step (e) of Algorithm $(A^1_{1/2})$ is seen in the right half of Figure 1. We believe that the meshes generated by Algorithm $(A^1_{1/2})$ have less local symmetry than that by $(A^0_{1/2})$. According to local extrapolation, symmetry could cause superconvergence phenomena. To check the practical convergence behaviour, we plotted in Figure 2 for each mesh $\mathcal{T}_k$ an entry $(N,e_N)$ and $(N,\eta _N)$. A log-scaling on both axes allows a slope $-\alpha$ of a straight line in the plot that connects two subsequent entries for a series of meshes $\mathcal{T}_0,\mathcal{T}_1,\mathcal{T}_2,...$ generated by Algorithm $(A^\vartheta _\Theta )$, to be interpreted as an experimental convergence rate $2\alpha$ (owing to $N\propto h^{-2}$ in two dimensions). We observe experimental convergence rates $2/3$, resp. $1$, for uniform, resp. adapted, meshes (generated by Algorithm $(A^\vartheta _\Theta )$ for $\Theta =0$, resp. $\Theta =1/2$). Furthermore, even for coarse meshes, $\eta _N$ appears to be a very good approximation to $e_N$; corresponding entries almost coincide for $(\Theta ,\vartheta ) = (1/2,0)$. If these meshes are perturbed (cf. Figure 1), the quotient $\eta _N/e_N$ is almost a constant very close to $1$. Numerical checks with different numerical quadrature rules (used to evaluate $e_N$) convinced us that, in general, $\eta _N$ behaves not asymptotically exact in practise but is very accurate.

In Example 7.2 we obtained meshes and experimental convergence rates displayed in Figures 3 and 4. Although $u$ belongs to $H^2(\Omega )$ and we expect linear convergence, $u$ has huge second order derivatives along a circular arc where $f$ is steep. We observe high refinements in the adapted meshes towards this arc. In this example the preasymptotic range is very large, an experimental convergence rate $1$ can be observed only for $N\ge 300$ for all refinement strategies. In this regime the estimator $\eta _Z$ appears as a good approximate for $e_N$ and the entries $(N,e_N)$ and $(N,\eta _N)$ almost conincide for $(\Theta ,\vartheta ) = (0,0),(1/2,0)$. This is not the case for $(\Theta ,\vartheta ) = (1/2,1)$, but the quotient $\eta _N/e_N\approx 0.9$ is still close to $1$.

7.2. Results for nonconforming finite element methods

The operator $\mathcal{A}$ from Example 5.1 serves in 7.1 to define $\eta _{Z,T}$ for first order Crouzeix–Raviart finite elements (cf., e.g., BSCi) in Algorithm $(A^\vartheta _\Theta )$. The generated meshes look similar to those shown in Figure 1, resp. Figure 3, and therefore are not displayed in this paper. The experimental convergence rates for Example 7.1, resp. 7.2, are illustrated in Figure 5, resp. Figure 6. The overall picture appears similarly to the above discussions and we draw the same conclusions. For uniform meshes, the quotient $\eta _N/e_N$ is nearly constant $1.2$ but significantly larger than $1.1$ in Figure 2.

7.3. Results for mixed finite element methods

For the Raviart–Thomas finite element method (cf. BBFBS) we use the operator $\mathcal{A}$ from Example 6.1 to define $\eta _{Z,T}$. The adapted meshes look similar to those shown in Figures 1 and 3 and therefore are not displayed in this paper. Figures 7 and 8 display the error and the estimator for the mixed finite element scheme in Examples 7.1 and 7.2 obtained from Algorithm $(A^\vartheta _\Theta )$ and $(\Theta ,\vartheta ) = (0,0),(1/2,0),(1/2,1)$. We obtain the same experimental convergence rates as in the previous methods.

7.4. Remarks

(i)

Our overall experience with Algorithm $(A^\vartheta _\Theta )$ and other (e.g., residual-based) adaptive algorithms supports that all such adaptive algorithms yield a considerable convergence improvement.

(ii)

Although asymptotic exactness of $\eta _N$ is not observed, the reliability constant $c_1$ in 1.5 and the efficiency constant are experimentally very close to $1$ since $\eta _N$ is a very good approximation to $e_N$ for very fine meshes (i.e., when h.o.t. is neglegible, say, for $N\ge 100$).

(iii)

Note that the efficiency constant is not known to be one as $\eta _N$ (based on the averaging operator $\mathcal{A}$) is different from$$\begin{equation} \eta _Z:=\min _{q_h \in \mathcal{S}^1_N(\mathcal{T},g)} \|p_h-q_h\|_{L^2(\Omega )}. {\tag{7.3}} \end{equation}$$

For conforming linear triangular finite elements, the efficiency of $\eta _N$ with an averaging operator follows from R1R2.

(iv)

Instead of the averaging operator $\mathcal{A}$, we tested the error estimator $\eta _Z \le \eta _N$ from 7.3 and found that sometimes the performance is poorer than that of $\eta _N$: In Figure 4, for instance, the results of $\eta _Z$ are much smaller than those of $\eta _N \approx e_N$. The averaging technique suggested in HSWW; average over a domain of size $O(h\log (1/h))$, is expected to give values between $\eta _N$ and $e_N$.

(v)

The error estimation in Example 7.1 is very accurate even for very coarse meshes. Hence the higher order terms do not seem to be important here although $u$ is nonsmooth. This agrees with our theoretical prediction in 1.5 since $f$ and $u_D$ are zero and $g$ is piecewise analytic. For a generic corner singularity of $u$, we expect$$\begin{equation*} \min _{q_h\in \mathcal{S}^1_N(\mathcal{T},g)} \|p-q_h\|_{L^2(\Omega )} \ll \min _{q_h\in \mathcal{L}^0(\mathcal{T})} \|p-q_h\|_{L^2(\Omega )} \le \|p-p_h\|_{L^2(\Omega )} \end{equation*}$$

and so the h.o.t. in 1.4 are not expected to be dominant even for coarse meshes.

(vi)

For coarse meshes in Example 7.2, higher order terms may cause the overall observation that $\eta _N$ is much smaller than $e_N$. Assuming $c_1\approx 1$ (which is seen for fine meshes, whence for neglegible h.o.t.) nonsmooth data ($\nabla f$ is large) indicate that $e_N \le c_1\eta _N +\text{h.o.t}$ cannot be improved to $e_N \le c_1 \eta _N$.

Acknowledgments

The second author (S.B.) thankfully acknowledges the partial support by the German research foundation at the Graduiertenkolleg 357 “Effiziente Algorithmen und Mehrskalenmethoden”.

Figure 1.

Adaptively refined meshes $\mathcal{T}_0$ (left upper) to $\mathcal{T}_8$ (right lower) (left) and perturbed triangulation $\mathcal{T}_{16}$ with $1157$ free nodes (right) in Example 7.1 for the conforming finite element scheme.

Figure 2.

Error and error estimator for uniform, adaptive, and perturbed adaptive mesh-refinement in the conforming finite element scheme in Example 7.1.

Figure 3.

Adaptively refined meshes $\mathcal{T}_0$ (left upper) to $\mathcal{T}_8$ (right lower) (left) and perturbed triangulation $\mathcal{T}_{12}$ with $1909$ free nodes (right) for the conforming finite element scheme in Example 7.2.

Figure 4.

Error and error estimator for uniform, adaptive, and perturbed adaptive mesh-refinement in the conforming finite element scheme in Example 7.2.

Figure 5.

Error and error estimator for uniform, adaptive, and perturbed adaptive mesh-refinement in the nonconforming finite element scheme in Example 7.1.

Figure 6.

Error and error estimator for uniform, adaptive, and perturbed adaptive mesh-refinement in the nonconforming finite element scheme in Example 7.2.

Figure 7.

Error and error estimator for uniform, adaptive, and perturbed adaptive mesh-refinement in the mixed finite element scheme in Example 7.1.

Figure 8.

Error and error estimator for uniform, adaptive, and perturbed adaptive mesh-refinement in the mixed finite element scheme in Example 7.2.