Local problems on stars: A posteriori error estimators, convergence, and performance

Abstract

A new computable a posteriori error estimator is introduced, which relies on the solution of small discrete problems on stars. It exhibits built-in flux equilibration and is equivalent to the energy error up to data oscillation without any saturation assumption. A simple adaptive strategy is designed, which simultaneously reduces error and data oscillation, and is shown to converge without mesh pre-adaptation nor explicit knowledge of constants. Numerical experiments reveal a competitive performance, show extremely good effectivity indices, and yield quasi-optimal meshes.

1. Introduction and main results

Adaptive finite element methods (FEM) have become essential tools in science and engineering for the numerical solution of multiscale phenomena governed by partial differential equations (PDE). In order to equidistribute the approximation errors, and thus the computational effort, adaptive FEM for elliptic PDE give rise to a sequence of graded meshes. This in turn leads to iterations of the form

$$\begin{equation} \text{solve }\to \text{ estimate } \to \text{ refine}. {\tag{1.1}} \end{equation}$$

A posteriori error estimators are the essence of estimate in 1.1. They are computable quantities, depending on the computed solution(s) and data, that provide information about the quality of approximation. They are problem-dependent and may be used to make judicious mesh refinements; coarsening is also important for evolution PDE. An efficient tool for local mesh refinement (coarsening) is the key component of refine in 1.1 and is typically problem-independent (see, e.g., 1516). We refer to 18 for references on adaptivity for elliptic PDE, and restrict the list of papers to those strictly related to our work.

The solution of local Dirichlet problems on stars was first proposed by Babuška and Miller in 3; a star is the support of a piecewise linear nodal basis function. The resulting error estimators are noncomputable since they require the solution of infinite dimensional problems of the same complexity as the original PDE. However, they are used in 3 to derive the first element-wise residual estimators in 2d with rigorous upper and lower bounds up to unknown interpolation constants. Bank and Weiser introduced in 4 a posteriori error estimators based on the solution of local Neumann problems on elements, which seem to allow for cancellation and thus lead to better results than the residual estimators. Their proofs of equivalence with the energy error require the saturation assumption though. The concept of flux equilibration was subsequently studied by Ainsworth and Oden 2 to enforce further cancellation, and thereby obtain better effectivity indices than 4. Equilibration consists of first solving linear systems on stars to compute suitable corrections to the weights of element jump residuals, and next of solving local problems on elements. We refer to 1 for the simplest proofs of upper and lower bounds, the former still requiring the solution of infinite dimensional local problems but being a constant free estimate. In 5 Carstensen and Funken also showed a constant free upper bound for an estimator based on the exact solution of weighted local problems on stars. Strouboulis et al. 17 derived guaranteed computable error bounds. In order to compute the estimators, they first solve equilibrated local Neumann problems with higher order polynomial degree, and then correct these local solutions via residual-type estimators. This yields an improved approximation to the exact solution of the local problems.

Computational experience strongly suggests that, starting from a coarse mesh, 1.1 converges within any prescribed tolerance in a finite number of steps. This issue has been recently tackled by Morin et al. 11 in the multidimensional setting exploiting an idea of Dörfler 6. The crucial role of data oscillation is disclosed in 11, which is intrinsic information missed by the averaging process associated with FEM regardless of quadrature. Ensuring a reduction rate of data oscillation in every step of 1.1, together with an error reduction, a linear rate of convergence of 1.1 for piecewise linear FEM is proven in 11 without any mesh pre-adaptation nor explicit knowledge of constants. This is achieved via a simple modification of the marking strategy of 6 which accounts for data oscillation and ensures an interior node in the refined mesh for every marked element. This theory is extended to any polynomial degree in 12.

The purpose of this paper is twofold. We first go back to the basis and propose a new computable a posteriori error estimator, with built-in flux equilibration, based on the solution of small discrete weighted problems on stars. Secondly, we devise an adaptive FEM which accounts for data oscillation, and prove that the resulting iteration 1.1 converges with a linear rate.

To describe the results in more detail, let $\Omega$ be a polygonal and bounded domain of $\mathbb{R}^2$, and let $u\in H^1(\Omega )$ be the solution of the following model problem

$$\begin{equation} \left\{ \begin{aligned} - \operatorname {div}(A \, \nabla u) &= f \qquad \text{in } \Omega , \\ u &= g \qquad \text{on } \partial \Omega , \end{aligned} \right. {\tag{1.2}} \end{equation}$$

where $f\in L^2(\Omega )$ and $g\in H^{1}(\partial \Omega )$. Furthermore, let $A$ be a piecewise constant, positive definite and symmetric matrix. Let ${\mathcal{T}}_h$ be a graded mesh with set of nodes ${\mathcal{N}}_h$. Since stars overlap, our current notion of data oscillation $\operatorname {osc}(f,g;{{\mathcal{N}}_h})$ differs from that in 11, and accounts for fine structure of both $f$ and $g$ which is averaged out by the FEM; see §3. Our first result reads as follows.

Main Result 1.

Let $u_h$ be the piecewise linear finite element solution over a triangulation ${\mathcal{T}}_h$. Then there exist three positive constants $C_1$, $C_2$, $C_3$, depending only on the minimum angle of ${\mathcal{T}}_h$ and $A$ such that the a posteriori error estimator $\mathcal{ E}_h$ based on the solution of discrete local problems on stars of §2 satisfies

$$\begin{equation} C_1 \mathcal{ E}_h \le \left|\negthinspace \left|\negthinspace \left|{u - u_h} \right|\negthinspace \right|\negthinspace \right| \le C_2 \mathcal{ E}_h + C_3 \operatorname {osc}(f,g;{{\mathcal{N}}_h}), {\tag{1.3}} \end{equation}$$

where $\left|\negthinspace \left|\negthinspace \left|{u-u_h} \right|\negthinspace \right|\negthinspace \right|$ denotes the energy norm of the error $u-u_h$.

Several remarks and comparisons with the existing literature are now in order:

•

The error estimator $\mathcal{ E}_h$ requires the solution of finite dimensional local problems on stars, and is thus computable. However, our proof of 1.3 does not make use of the saturation assumption as in 24. Such an assumption is shown to be superfluous in 14 for the estimator of 4, and indeed related to the relative size of data oscillation in 7. However, removing the saturation assumption requires comparison with residual estimators, which exhibit a notorious worse performance, thereby weakening the result. An important new feature of our proof of 1.3 is that it applies directly without reference to residual estimators.

•

Our estimators possess a self-equilibration property which result from direct exploitation of Galerkin orthogonality on stars. Therefore, we do not need to equilibrate jump residuals as in 12.

•

Our numerical experiments of §8 for the case $A=I$ show effectivity indices $\mathcal{ E}_h/\left|\negthinspace \left|\negthinspace \left|{u-u_h} \right|\negthinspace \right|\negthinspace \right|$ very close to unity as the mesh ${\mathcal{T}}_h$ is refined adaptively. In case $A$ has discontinuous coefficients with jumps of order $160$, and a solution barely in $H^1(\Omega )$, effectivity indices are around $0.7$ which means an error underestimation of about 30%.

•

The cost of error estimation is slightly higher than solving local problems on elements but the estimators are better than those of 4; see §8. However, when added, the computing time of flux equilibration, which entails solving a linear system for every interior star and a local problem for every element as in 4, the performance of our adaptive FEM is quite competitive.

•

The error estimator $\mathcal{ E}_h$ is computed with contributions from all the stars of ${\mathcal{T}}_h$, including boundary stars. As usual, the upper bound in (1.3) is global, whereas the lower bound is local and valid for every star.

•

Since there is no data oscillation in the lower bound of (1.3), we end up with a strong concept of reliability: the relative size of star indicators dictates mesh refinement regardless of fine structure of both $f$ and $g$, and thus without resorting to asymptotics.

Our second main result deals with convergence of 1.1 as the adaptive procedure Algorithm C of §4 is iterated.

Main Result 2.

Let $u_k$ be a sequence of finite element solutions over nested triangulations ${\mathcal{T}}_k$ produced by Algorithm C. There exist positive constants $C_0$ and $\alpha <1$, depending only on given data and the initial grid, such that

$$\begin{equation} \left|\negthinspace \left|\negthinspace \left|{u - u_k} \right|\negthinspace \right|\negthinspace \right| \le C_0 \,\alpha ^k. {\tag{1.4}} \end{equation}$$

The initial coarse mesh need not be adjusted to resolve data $f$ and $g$ to any tolerance, and no explicit constants are needed for Algorithm C to work.

Several remarks are now in order:

•

The basic ingredients for this result to hold are as follows: a global upper bound of the error in terms of the estimators (such as Theorem 3.6); a local lower bound for the difference between two consecutive discrete solutions (such as Lemma 5.2); nested meshes (which yield an energy error decrease such as Lemma 5.1); and a suitable marking strategy (such as that in Proposition 4.1). The linear decay rate (1.4) is thus valid for every adaptive procedure with these properties. This is the case of the procedure of 11 based on residual type error estimators.

•

Depending on the flatness of $u$, the mesh size of ${\mathcal{T}}_k$ may not necessarily tend to zero as $k\uparrow \infty$, which makes 1.4 a nonstandard finite element asymptotic statement.

•

Even though no stability constants are required for Algorithm C, nor for convergence, constants $C_2$ and $C_3$ in 1.3 are needed to stop the iterations; this is customary in adaptivity 18.

•

Main Result 2 is stronger than that of 11 in regard to data oscillation. We no longer impose a data oscillation reduction rate but rather that $\operatorname {osc}(f,g;{{\mathcal{N}}_k})$ stays bounded by a decreasing exponential of $k$. This simple modification leads to a dramatic improvement in that data oscillation plays now a much weaker computational role than in 11; see §8 for details.

The rest of the paper is organized as follows. In §2 we introduce some notation, motivate the definition of the a posteriori error estimator, and define the estimator $\mathcal{ E}_h$ based on star indicators $\mathcal{ E}_i$. In §3 we prove the equivalence between error and estimator (Main Result 1). In §4 we prove the convergence of 1.1 (Main Result 2). The proof hinges on a crucial error reduction property and an analogous data oscillation reduction property, which are then shown in §§5 and 6, respectively. In §7 we demonstrate an auxiliary weighted Poincaré inequality, which is a crucial technical devise. We conclude with several numerical experiments in §8. They reveal a competitive performance of Algorithm C in terms of both computing time and effectivity indices, and show that the graded meshes possess quasi-optimal number of degrees of freedom.

2. A posteriori error estimators

We start this section with some useful notation. For an open set $G \subset \mathbb{R}^2$ we denote by $H^1(G)$ the usual Sobolev space of functions in $L^2(G)$ whose first derivatives are also in $L^2(G)$, and by $H^1_{\!g}(G)$ the subspace of functions in $H^1(G)$ which are equal to $g$ on the boundary, both spaces endowed with the norm

$$\begin{equation*} \left\| u \right\|_{H^1(G)} := \left( \|u\|_G^2 + \|\nabla u\|_G^2 \right)^{1/2} \end{equation*}$$

where $\|\cdot \|_G$ stands for the $L^2(G)$-norm.

Since $A$ is piecewise constant, symmetric and positive definite, the bilinear form $\left(\hskip -1pt\left({\cdot }\, , \,{\cdot }\right)\hskip -1pt\right)_G$ defined for any open set $G \subset \Omega$ by

$$\begin{equation*} \left(\hskip -1pt\left({u}\, , \,{v}\right)\hskip -1pt\right)_G := \int _G A \nabla u \cdot \nabla v, \end{equation*}$$

is bounded in $H^1(G)$, and $\left(\hskip -1pt\left({\cdot }\, , \,{\cdot }\right)\hskip -1pt\right)_\Omega$ is coercive on $H_0^1(\Omega )$, i.e., there exist constants $0 < c_a \le C_a < \infty$ such that

$$\begin{gather} \left(\hskip -1pt\left({v}\, , \,{w}\right)\hskip -1pt\right)_G \le C_a \left\| v \right\|_{H^1(G)} \left\| w \right\|_{H^1(G)}, \qquad \forall\nobreakspace v, w \in H^1(G), {\tag{2.1}}\\ c_a \left\| v \right\|_{H^1(\Omega )}^2 \le \left(\hskip -1pt\left({v}\, , \,{v}\right)\hskip -1pt\right)_\Omega , \qquad \forall\nobreakspace v \in H_0^1(\Omega ). {\tag{2.2}} \end{gather}$$

This implies, in particular, that the seminorm $\left|\negthinspace \left|\negthinspace \left|{\cdot } \right|\negthinspace \right|\negthinspace \right|_G$ defined by $\left|\negthinspace \left|\negthinspace \left|{v} \right|\negthinspace \right|\negthinspace \right|_G^2 := \left(\hskip -1pt\left({v}\, , \,{v}\right)\hskip -1pt\right)_G$ is equivalent to the $H_0^1(\Omega )$-norm when $G = \Omega$.

In view of (2.1) and (2.2), problem (1.2) admits a unique weak solution $u$ for any $f\in L^2(\Omega )$ and $g \in H^1(\partial \Omega )$, i.e., there exists a unique $u \in H_g^1(\Omega )$ satisfying

$$\begin{equation} \left(\hskip -1pt\left({u}\, , \,{v}\right)\hskip -1pt\right)_\Omega = \int _\Omega f v, \quad \text{for all } v \in H_0^1(\Omega ). {\tag{2.3}} \end{equation}$$

From now on, whenever it is clear from the context, we will omit the subscript $\Omega$ from $\left(\hskip -1pt\left({\cdot }\, , \,{\cdot }\right)\hskip -1pt\right)_\Omega$ and $\left|\negthinspace \left|\negthinspace \left|{\cdot } \right|\negthinspace \right|\negthinspace \right|_\Omega$.

Let ${\mathcal{T}}_h$ be a conforming shape-regular triangulation of $\Omega$. Let ${\mathbb{V}}_{\!h}$ be the space of continuous piecewise linear functions over ${\mathcal{T}}_h$, and let ${\mathbb{V}}_{\!h}^{g_h}$ be the subspace of functions of ${\mathbb{V}}_{\!h}$ that are equal to $g_h$ at the boundary, where $g_h$ is the piecewise linear interpolant of $g$. The finite element approximation $u_h$ to $u$ is defined by

$$\begin{equation} u_h \in {\mathbb{V}}_{\!h}^{g_h}: \quad \int _\Omega A \nabla {u_h}\nabla {\phi } = \int _\Omega f \phi , \quad \forall\nobreakspace \phi \in {\mathbb{V}}_{\!h}^0. {\tag{2.4}} \end{equation}$$

We denote with $\left\{x_i\right\}_{i\in {\mathcal{N}}_h}$ the set of all nodes of the triangulation ${\mathcal{T}}_h$. For each $i \in {\mathcal{N}}_h$, $\phi _i$ denotes the canonical piecewise linear basis function corresponding to $x_i$. The star ${\omega _i}$ is the interior relative to $\bar{\Omega }$ of the support of $\phi _i$, and $h_i$ is the maximum size of the elements forming ${\omega _i}$. Finally, ${\Gamma _i}$ will denote the union of the sides touching $x_i$ that are contained in $\Omega$ (see Figure 1).

As a motivation, we consider momentarily the case $g = 0$ and recall the local error indicators $\zeta _i$, introduced in 5, which cannot be computed exactly, but will give us an idea of how to proceed in order to define computable error indicators.

Definition 2.1.

For each star ${\omega _i}$, $i \in {\mathcal{N}}_h$ let $H^1_{\text{loc}}({\omega _i})$ be the space of functions in ${\omega _i}$ that are in $H^1(K)$ for all compact $K\subset \{\phi _i>0\}\cap \bar{\Omega }$. Next we define

$$\begin{equation*} W({\omega _i})= \left\{ v \in H^1_{\text{loc}}({\omega _i}): \int _{\omega _i}v \,\phi _i = 0 \text{ and } \int _{\omega _i}\left|\nabla v\right|^2 \phi _i < \infty \right\}, \end{equation*}$$

if $x_i$ is an interior node, and

$$\begin{equation*} W({\omega _i})= \left\{ v \in H^1_{\text{loc}}({\omega _i}): v = 0 \text{ on } \partial {\omega _i}\cap \partial \Omega \text{ and } \int _{\omega _i}\left|\nabla v\right|^2 \phi _i < \infty \right\} \end{equation*}$$

otherwise. Next we define $\zeta _i \in W({\omega _i})$ to be the solution to

$$\begin{equation} \int _{\omega _i}A \nabla \zeta _i\nabla \phi \, \phi _i = \int _{\omega _i}f \phi \, \phi _i + \int _{\Gamma _i}J \phi \, \phi _i \quad \text{for all } \phi \in W({\omega _i}), {\tag{2.5}} \end{equation}$$

where $J = [A \nabla u_h\cdot \nu ]$ stands for the jump of flux across the sides of the triangulation, which is independent of the orientation of the unit normal $\nu$. We then define the global error estimator $\tilde{\mathcal{ E}}_h = \bigl ( \sum _{i\in {\mathcal{N}}_h} \tilde{{\mathcal{ E}_i}}^2\bigr )^{1/2}$ in terms of the local error indicators

$$\begin{equation} \tilde{{\mathcal{ E}_i}}^2 = \left|\negthinspace \left|\negthinspace \left|{{\zeta _i}} \right|\negthinspace \right|\negthinspace \right|_{\phi _i}^2 := \int _\Omega A \nabla \zeta _i \nabla \zeta _i \, \phi _i. {\tag{2.6}} \end{equation}$$

Remark 2.2.

There is a minor difference between the local spaces used here and those used in 5. For interior nodes, the functions $v \in W({\omega _i})$ in 5 are required to satisfy $\int _{\omega _i}v = 0$ instead of $\int _{\omega _i}v \phi _i = 0$.

Remark 2.3.

The existence of a solution to equation (2.5) is guaranteed by the Lax-Milgram lemma. In fact, it is immediate to see that the bilinear form in (2.5) is coercive and bounded in $W({\omega _i})$ with respect to the norm $\left|\negthinspace \left|\negthinspace \left|{{\cdot }} \right|\negthinspace \right|\negthinspace \right|_{\phi _i}$, whereas the boundedness of the right hand side of this equation is a consequence of the trace theorem and the following proposition, whose proof is postponed to §7.

Proposition 2.4 (Weighted Poincaré inequality).

The space $W({\omega _i})$ endowed with the norm $\left|\negthinspace \left|\negthinspace \left|{{\cdot }} \right|\negthinspace \right|\negthinspace \right|_{\phi _i}$ is continuously embedded in $L^2({\omega _i})$ for any $i \in {\mathcal{N}}_h$. More specifically, there exists a constant $C$, depending only on the minimum angle of the triangulation but otherwise independent of the star being considered such that

$$\begin{equation} \left\| v \right\|_{L^2({\omega _i})} \le C \, h_i \, \left\| \nabla v \right\|_{\phi _i}, {\tag{2.7}} \end{equation}$$

where $\left\| \cdot \right\|_{\phi _i} = \bigl \| \cdot \,\phi _i^{1/2}\bigr \|_{L^2({\omega _i})}$.

Now, if we denote by $e$ the error $u - u_h$, for any $\phi \in H_0^1(\Omega )$, we have

$$\begin{align*} \left(\hskip -1pt\left({e}\, , \,{\phi }\right)\hskip -1pt\right) &= \sum _{i\in {\mathcal{N}}_h} \left(\hskip -1pt\left({u-u_h}\, , \,{\phi _i \, \phi }\right)\hskip -1pt\right) = \sum _{i\in {\mathcal{N}}_h} \left[ \int _{\omega _i}f \phi \, \phi _i + \int _{\Gamma _i}J \phi \, \phi _i \right] \\ &= \sum _{i\in {\mathcal{N}}_h} \left[ \int _{\omega _i}f (\phi - c_i) \phi _i + \int _{\Gamma _i}J (\phi -c_i) \phi _i \right], \end{align*}$$

where $c_i = 0$ for boundary nodes and $\displaystyle c_i = \textstyle {(\int _{\omega _i}\phi _i)^{-1}\int _{\omega _i}\phi \, \phi _i}$ otherwise. The last step follows from the Galerkin orthogonality

$$\begin{equation*} \int _{\omega _i}f \phi _i + \int _{\Gamma _i}J \phi _i = 0, \end{equation*}$$

which holds for every interior node $i$. Since $\phi - c_i \in W({\omega _i})$ we get by the definition of $\zeta _i$,

$$\begin{align*} \left(\hskip -1pt\left({e}\, , \,{\phi }\right)\hskip -1pt\right) &= \sum _{i\in {\mathcal{N}}_h} \int _{\omega _i}A \nabla \zeta _i \nabla \phi \, \phi _i \le \sum _{i\in {\mathcal{N}}_h} \left|\negthinspace \left|\negthinspace \left|{{\zeta _i}} \right|\negthinspace \right|\negthinspace \right|_{\phi _i} \left|\negthinspace \left|\negthinspace \left|{{\phi }} \right|\negthinspace \right|\negthinspace \right|_{\phi _i} \\ &\le \left(\sum _{i\in {\mathcal{N}}_h} \left|\negthinspace \left|\negthinspace \left|{{\zeta _i}} \right|\negthinspace \right|\negthinspace \right|_{\phi _i}^2 \right)^{1/2} \left(\sum _{i\in {\mathcal{N}}_h} \int _\Omega A \nabla \phi \nabla \phi \,\phi _i \right)^{1/2} = \tilde{\mathcal{ E}}_h \left( \int _\Omega A \nabla \phi \nabla \phi \right)^{1/2}. \end{align*}$$

Since $g=0$, taking $\phi = e \in H_0^1(\Omega )$ we obtain $\left|\negthinspace \left|\negthinspace \left|{e} \right|\negthinspace \right|\negthinspace \right|^2 \le \tilde{\mathcal{ E}}_h \left|\negthinspace \left|\negthinspace \left|{e} \right|\negthinspace \right|\negthinspace \right|$, which implies the constant free estimate

$$\begin{equation*} \left|\negthinspace \left|\negthinspace \left|{e} \right|\negthinspace \right|\negthinspace \right| \le \tilde{\mathcal{ E}}_h. \end{equation*}$$

Despite this constant free upper bound, the estimator has an important drawback: it is not computable because it requires the knowledge of exact solutions to local problems. These problems are located on smaller domains, but they are still of an infinite dimensional nature. In order to circumvent such a difficulty, we now define finite dimensional local spaces ${\mathcal{P}}_0^2({\omega _i})$. This leads to computable error indicators ${\mathcal{ E}_i}$, which will be used in the subsequent analysis.

Definition 2.5.

For $i \in {\mathcal{N}}_h$, let ${\mathcal{P}}^2({\omega _i})$ denote the space of piecewise quadratic functions on the star ${\omega _i}$ that vanish on $\partial {\omega _i}$. Let ${\mathcal{P}}_0^2({\omega _i})= {\mathcal{P}}^2({\omega _i})\cap W({\omega _i})$, that is ${\mathcal{P}}_0^2({\omega _i})$ coincides with ${\mathcal{P}}^2({\omega _i})$ for boundary nodes, and for interior nodes ${\mathcal{P}}_0^2({\omega _i})$ is the subspace of functions $\phi \in {\mathcal{P}}^2({\omega _i})$ satisfying $\int _{\omega _i}\phi \phi _i = 0$. For each star ${\omega _i}$, $i \in {\mathcal{N}}_h$, we now define $\eta _i \in {\mathcal{P}}_0^2({\omega _i})$ to be the solution to

$$\begin{equation} \int _{\omega _i}A \nabla \eta _i \nabla \phi \, \phi _i = \int _{\omega _i}f \phi \, \phi _i + \int _{\Gamma _i}J \phi \, \phi _i \quad \text{for all } \phi \in {\mathcal{P}}_0^2({\omega _i}), {\tag{2.8}} \end{equation}$$

where $J$ denotes, as before, the jump of flux across the sides of the triangulation. We then define the global error estimator $\mathcal{ E}_h = \bigl ( \sum _{i\in {\mathcal{N}}_h} {\mathcal{ E}_i}^2\bigr )^{1/2}$ in terms of the local error indicators

$$\begin{equation} {\mathcal{ E}_i}^2 = \left|\negthinspace \left|\negthinspace \left|{{\eta _i}} \right|\negthinspace \right|\negthinspace \right|_{\phi _i}^2 . {\tag{2.9}} \end{equation}$$

Remark 2.6.

The computation of the estimators just defined requires the solution of a small linear system for each star ${\omega _i}$ instead of an infinite dimensional eigenvalue problem as in 5. This issue is further discussed in §8.

Remark 2.7.

Since the functions $\phi \in {\mathcal{P}}_0^2({\omega _i})$ vanish on $\partial {\omega _i}$, satisfaction of $\int _{\omega _i}\phi \phi _i=0$ cannot be achieved by adding a convenient constant and thus does have an effect in the computation of ${\mathcal{ E}_i}$ of (2.9) (compare with Remark 2.2). In contrast with 5, our numerical experiments of §8 show a superior performance and effectivity indices close to 1 for the Laplacian.

3. Equivalence

In this section we prove the first of the two main results of this article: the equivalence (up to oscillation terms) of the error and the new estimator introduced in 2.9. We do this in two steps, we first prove a lower bound for the error without oscillation, and afterwards an upper bound.

Theorem 3.1 (Lower bound).

There exists a constant $C_1>0$, depending on the minimum angle of the triangulation and $A$ such that, for any $i \in {\mathcal{N}}_h$,

$$\begin{equation} C_1 {\mathcal{ E}_i}\le \left|\negthinspace \left|\negthinspace \left|{u-u_h} \right|\negthinspace \right|\negthinspace \right|_{\omega _i}. {\tag{3.1}} \end{equation}$$

Remark 3.2.

This data oscillation free lower bound implies a strong concept of reliability: the relative size of ${\mathcal{ E}_i}$ dictates mesh refinement regardless of fine structure of $f$ and $g$, and thus works even in the pre-asymptotic regime.

Remark 3.3.

Consider the following modified definition of the local indicators:

$$\begin{equation*} \eta _i \in {\mathcal{P}}_0^2({\omega _i}): \quad \int _{\omega _i}A \nabla \eta _i \nabla \phi \, \phi _i = \int _{\Gamma _i}J \phi \, \phi _i \quad \text{for all } \phi \in {\mathcal{P}}_0^2({\omega _i}), \end{equation*}$$

instead of (2.8). Then all the results of this paper remain true except (3.1), because the star oscillation of $f$ defined below in 3.2, namely $h_i \| (f - f_i) \phi _i^{1/2}\|_{\omega _i}$, would also appear on the right hand side.

Proof of Theorem 3.1.

Let $i\in {\mathcal{N}}_h$, and let ${\mathcal{ E}_i}$ and $\eta _i$ be as in Definition 2.5. Then

$$\begin{align*} {\mathcal{ E}_i}^2 &= \left|\negthinspace \left|\negthinspace \left|{{\eta _i}} \right|\negthinspace \right|\negthinspace \right|_{\phi _i}^2 = \int _{\omega _i}A \nabla \eta _i \nabla \eta _i \,\phi _i = \int _{\omega _i}f \eta _i \phi _i + \int _{\Gamma _i}J \eta _i \phi _i \\ &= \int _{\omega _i}A \nabla (u-u_h) \nabla (\eta _i \phi _i) = \int _{\omega _i}A \nabla (u-u_h) \nabla \eta _i \phi _i + \int _{\omega _i}A \nabla (u-u_h) \eta _i \nabla \phi _i. \end{align*}$$

Now, $\int _{\omega _i}A \nabla (u-u_h) \nabla \eta _i \phi _i \le \left|\negthinspace \left|\negthinspace \left|{{u-u_h}} \right|\negthinspace \right|\negthinspace \right|_{\phi _i} \left|\negthinspace \left|\negthinspace \left|{{\eta _i}} \right|\negthinspace \right|\negthinspace \right|_{\phi _i}$, and by Proposition 2.4

$$\begin{equation*} \int _{\omega _i}A \nabla (u-u_h) \eta _i \nabla \phi _i \le C \left|\negthinspace \left|\negthinspace \left|{u-u_h} \right|\negthinspace \right|\negthinspace \right|_{\omega _i}\frac{1}{h_i} \| \eta _i \|_{\omega _i}\le C \left|\negthinspace \left|\negthinspace \left|{u-u_h} \right|\negthinspace \right|\negthinspace \right|_{\omega _i}\left|\negthinspace \left|\negthinspace \left|{{\eta _i}} \right|\negthinspace \right|\negthinspace \right|_{\phi _i}, \end{equation*}$$

which immediately implies the desired result.

■

Before stating the upper bound for the error, we need to define two quantities, related to data oscillation. First we define for ${\mathcal{M}}_h\subset {\mathcal{N}}_h$

$$\begin{equation} \operatorname {osc}(f;{{\mathcal{M}}_h}) := \left( \sum _{\,i \in {\mathcal{M}}_h} h_i^2 \bigl \| (f - f_i) \phi _i^{1/2} \bigr \|_{\Omega }^2 \right)^{1/2}, {\tag{3.2}} \end{equation}$$

where $\displaystyle f_i = \frac{\int _{\omega _i}f \phi _i }{ \int _{\omega _i}\phi _i}$ for interior nodes, and $f_i = 0$ otherwise. Secondly, we define

$$\begin{equation} \operatorname {osc}_\partial (g;{{\mathcal{M}}_h}) := \Biggl ( \sum _{\,i \in {\mathcal{M}}_h,\, x_i \in \partial \Omega } |S_i| \, \bigl \|\partial _\tau (g -g_h) \bigl \|_{L^2(S_i)}^2 \Biggr )^{1/2}, {\tag{3.3}} \end{equation}$$

where $S_i$ denotes the union of the two sides touching $x_i$ that lie on the boundary of $\Omega$ and $|S_i|$ denotes the length of $S_i$. This $\operatorname {osc}_\partial (g;{{\mathcal{N}}_h})$ measures information missing when approximating boundary data by piecewise linear functions (see Lemma 3.4 below). These oscillation terms are generically of higher order than the error. We will prove that they get reduced by a fixed proportion of themselves using an appropriate selection of elements for refinement. This will be a key step in proving convergence of the Adaptive Algorithm C.

The following lemma states the relation between $\operatorname {osc}_\partial (g;{{\mathcal{N}}_h})$ and the information of $g$ missed by $g_h$.

Lemma 3.4.

Let $\tilde{u}$ be the solution of the continuous problem 2.3 with boundary data $g_h$ instead of $g$. Then there exists a constant $C_\partial$ depending on the minimum angle of the triangulation and $A$ such that

$$\begin{equation} \left|\negthinspace \left|\negthinspace \left|{u-\tilde{u}} \right|\negthinspace \right|\negthinspace \right| \le C_\partial \operatorname {osc}_\partial (g;{{\mathcal{N}}_h}). {\tag{3.4}} \end{equation}$$

Remark 3.5.

This result is equivalent to the bound

$$\begin{equation*} \| g-g_h\|_{H^{1/2}(\partial \Omega )} \le C_\partial \operatorname {osc}_\partial (g;{{\mathcal{N}}_h}), \end{equation*}$$

which is already known. For completeness and since the proof is simple, we present it here. The dependence on the minimal angle of ${\mathcal{T}}_h$ is a consequence of the method of proof.

Proof of Lemma 3.4.

By definition, $u-\tilde{u} \in H^1_{g-g_h}(\Omega )$ satisfies $\left(\hskip -1pt\left({u-\tilde{u}}\, , \,{ v}\right)\hskip -1pt\right) = 0$ for all $v \in H_0^1(\Omega )$. Consequently $\left|\negthinspace \left|\negthinspace \left|{u-\tilde{u}} \right|\negthinspace \right|\negthinspace \right| \le \left|\negthinspace \left|\negthinspace \left|{v} \right|\negthinspace \right|\negthinspace \right|$ for all $v \in H^1_{g-g_h}(\Omega )$.

Let us now construct a particular extension $v$ of $g-g_h$ to $\Omega$. If $T \in {\mathcal{T}}_h$ has one side on $\partial \Omega$, then let $F$ be a linear mapping from $T$ onto the reference triangle $\hat{T} = \{ (x,y) : 0\le x\le 1, 0\le y\le x\}$ such that $F(\partial T \cap \partial \Omega ) = \{ 0\le x\le 1, y=0\}$. Let $\hat{g} = (g-g_h)\circ F^{-1} \in H_0^1(0,1)$, and let $\hat{v}$ be the extension of $\hat{g}$ to $\hat{T}$ given by

$$\begin{equation*} \hat{v}(x,y) = \left( 1 - \frac{y}{x}\right) \hat{g}(x), \end{equation*}$$

i.e., for every $x_0 \in [0,1]$, $\hat{v}(x_0,\cdot )$ is the linear interpolant on the segment $\{x=x_0, 0\le y \le x_0\}$ of $\hat{g}(x_0)$ and $0$. Then a straightforward calculation shows that $\int _{\hat{T}} |\nabla \hat{v}|^2 \le C \int _0^1 |\hat{g}'(x)|^2 \,dx$. Defining $v_{|T} = \hat{v}_{|\hat{T}} \circ F$ we obtain

$$\begin{equation*} v_{| \partial T \cap \partial \Omega } = g-g_h, \quad v_{| \partial T \setminus \partial \Omega } = 0, \quad \int _T |\nabla v|^2 \le C h_T \int _{\partial T \cap \partial \Omega } | \partial _\tau (g-g_h) |^2. \end{equation*}$$

If $T$ has two sides in $\partial \Omega$ we divide it by two and apply the above argument to each of the two sub-triangles. If $T$ has no side on $\partial \Omega$ we define $v_{|T}=0$. Consequently, $v \in H^1_{g-g_h}(\Omega )$ and $v_{|T} = 0$ for all $T$ having no side on $\partial \Omega$. Then

$$\begin{equation*} \left|\negthinspace \left|\negthinspace \left|{u-\tilde{u}} \right|\negthinspace \right|\negthinspace \right|^2 \le \left|\negthinspace \left|\negthinspace \left|{v} \right|\negthinspace \right|\negthinspace \right|^2 \le C_a \sum _{\substack{T \in {\mathcal{T}}_h\\ \partial T \cap \partial \Omega \neq \emptyset }} \| \nabla v \|_{T}^2 \le C \sum _{\substack{T \in {\mathcal{T}}_h\\ \partial T \cap \partial \Omega \neq \emptyset }} h_T \| \partial _\tau (g-g_h) \|_{\partial T \cap \partial \Omega }^2, \end{equation*}$$

which implies the claim.

■

We now state and prove the second main result of this section.

Theorem 3.6 (Upper bound).

There exist two positive constants $C_2$ and $C_3$, depending only on the minimum angle of ${\mathcal{T}}_h$ and $A$ such that

$$\begin{equation*} \left|\negthinspace \left|\negthinspace \left|{u - u_h } \right|\negthinspace \right|\negthinspace \right|_\Omega \le C_2 \mathcal{ E}_h + C_3 \operatorname {osc}(f,g;{{\mathcal{N}}_h}), \end{equation*}$$

where $\operatorname {osc}(f,g;{{\mathcal{N}}_h}) := \left( \operatorname {osc}(f;{{\mathcal{N}}_h})^2 + \operatorname {osc}_\partial (g;{{\mathcal{N}}_h})^2 \right)^{1/2}$.

To prove this theorem we will need the following lemma, whose proof is postponed to the end of this section.

Lemma 3.7.

For each node $i \in {\mathcal{N}}_h$ there exists an operator ${\Pi _i}:W({\omega _i})\to {\mathcal{P}}_0^2({\omega _i})$, such that for any $v \in W({\omega _i})$ the following conditions hold:

(i)

$\int _S ( v - {\Pi _i}v ) \phi _i = 0$, for all $S \subset {\Gamma _i}$,

(ii)

$\| \nabla {\Pi _i}v \|_{\phi _i} \le C^* \| \nabla v \|_{\phi _i}$,

where the constant $C^*$ depends only on the minimum angle of ${\mathcal{T}}_h$.

We are now ready to prove the upper bound.

Proof of Theorem 3.6.

Let $\tilde{u}$ be the solution to the continuous problem (1.2) with boundary data $g_h$. By Lemma (3.4), $\left|\negthinspace \left|\negthinspace \left|{\tilde{u} - u} \right|\negthinspace \right|\negthinspace \right| \le C_\partial \operatorname {osc}_\partial (g;{{\mathcal{N}}_h})$. Let us now denote with $\tilde{e}$ the error $\tilde{u}-u_h$. If $\phi \in H_0^1(\Omega )$, then as before

$$\begin{align*} \left(\hskip -1pt\left({\tilde{e}}\, , \,{\phi }\right)\hskip -1pt\right) &= \sum _{i\in {\mathcal{N}}_h} \left[ \int _{\omega _i}f (\phi - c_i) \phi _i + \int _{\Gamma _i}J (\phi -c_i) \phi _i \right], \end{align*}$$

where $c_i = 0$ for boundary nodes and $\displaystyle c_i = \frac{\int _{\omega _i}\phi \, \phi _i}{ \int _{\omega _i}\phi _i}$ otherwise. By Lemma 3.7(i), since $\phi - c_i \in W({\omega _i})$, we have

$$\begin{align*} \left(\hskip -1pt\left({\tilde{e}}\, , \,{\phi }\right)\hskip -1pt\right) &= \sum _{i\in {\mathcal{N}}_h} \left[ \int _{\omega _i}f {\Pi _i}(\phi - c_i) \phi _i + \int _{\Gamma _i}J {\Pi _i}(\phi -c_i) \phi _i \right] \\ &\quad + \sum _{i\in {\mathcal{N}}_h} \int _{\omega _i}f \bigl ( (\phi - c_i) - {\Pi _i}(\phi - c_i)\bigr ) \phi _i \\ &= \sum _{i\in {\mathcal{N}}_h} \int _{\omega _i}A \nabla \eta _i \nabla {\Pi _i}(\phi - c_i) \phi _i + \sum _{i\in {\mathcal{N}}_h} \int _{\omega _i}f \bigl ( (\phi - c_i) - {\Pi _i}(\phi - c_i)\bigr ) \phi _i, \end{align*}$$

where the last equality follows from Definition 2.5. Using Lemma 3.7(ii), together with the Hölder and Cauchy-Schwarz inequalities, we obtain

$$\begin{align*} \sum _{i\in {\mathcal{N}}_h} \int _{\omega _i}A \nabla \eta _i \nabla {\Pi _i}(\phi - c_i) \phi _i &\le \left(\sum _{i\in {\mathcal{N}}_h} \left|\negthinspace \left|\negthinspace \left|{{\eta _i}} \right|\negthinspace \right|\negthinspace \right|_{\phi _i}^2 \right)^{\frac{1}{2}}\left(\sum _{i\in {\mathcal{N}}_h} \left|\negthinspace \left|\negthinspace \left|{{{\Pi _i}(\phi - c_i)}} \right|\negthinspace \right|\negthinspace \right|_{\phi _i}^2 \right)^{\frac{1}{2}}\\ &\le C_{\!A} C^* \left(\sum _{i\in {\mathcal{N}}_h} {\mathcal{ E}_i}^2 \right) ^{\frac{1}{2}}\left(\sum _{i\in {\mathcal{N}}_h} \left|\negthinspace \left|\negthinspace \left|{{\phi }} \right|\negthinspace \right|\negthinspace \right|_{\phi _i} \right)^{\frac{1}{2}}= C_{\!A} C^* \mathcal{ E}_h \left|\negthinspace \left|\negthinspace \left|{\phi } \right|\negthinspace \right|\negthinspace \right| . \end{align*}$$

Observing that both $(\phi - c_i)$ and ${\Pi _i}(\phi - c_i)$ are in $W({\omega _i})$ and thus have weighted mean value zero for interior nodes, we have

$$\begin{equation*} \int _{\omega _i}f \bigl ( (\phi - c_i) -{\Pi _i}(\phi - c_i)\bigr )\phi _i = \int _{\omega _i}(f-f_i) \bigl ( (\phi - c_i) -{\Pi _i}(\phi - c_i)\bigr )\phi _i \end{equation*}$$

for any constant $f_i \in \mathbb{R}$ for all interior nodes. Defining $f_i=0$ for boundary nodes, we have

$$\begin{equation*} \int _{\omega _i}f \bigl ( (\phi - c_i) -{\Pi _i}(\phi - c_i)\bigr )\phi _i \le \left\| f - f_i \right\|_{\phi _i} \left( \left\|\phi - c_i\right\|_{\phi _i} + \left\|{\Pi _i}(\phi - c_i)\right\|_{\phi _i} \right), \end{equation*}$$

for every node. Using again Lemma 3.7(ii) and Proposition 2.4, we obtain

$$\begin{equation*} \sum _{i\in {\mathcal{N}}_h} \int f \bigl ( (\phi - c_i) -{\Pi _i}(\phi - c_i)\bigr )\phi _i \le C \left(\sum _{i\in {\mathcal{N}}_h} \left\| h_i (f - f_i)\right\|_{\phi _i}^2\right)^{1/2} \| \nabla \phi \|_\Omega . \end{equation*}$$

Summarizing, for any $\phi \in H_0^1(\Omega )$, we have

$$\begin{align*} \left(\hskip -1pt\left({e}\, , \,{\phi }\right)\hskip -1pt\right) &\le C_2 \mathcal{ E}_h \left|\negthinspace \left|\negthinspace \left|{\phi } \right|\negthinspace \right|\negthinspace \right| + C_3 \left(\sum _{i\in {\mathcal{N}}_h} \left\| h_i (f - f_i)\right\|_{\phi _i}^2\right)^{1/2} \left|\negthinspace \left|\negthinspace \left|{\phi } \right|\negthinspace \right|\negthinspace \right|, \end{align*}$$

and thus $\left|\negthinspace \left|\negthinspace \left|{\tilde{e}} \right|\negthinspace \right|\negthinspace \right| \le C_2 \mathcal{ E}_h + C_3 \operatorname {osc}(f;{{\mathcal{N}}_h})$. The claim follows from the triangle inequality $\left|\negthinspace \left|\negthinspace \left|{e} \right|\negthinspace \right|\negthinspace \right| \le \left|\negthinspace \left|\negthinspace \left|{\tilde{e}} \right|\negthinspace \right|\negthinspace \right| + \left|\negthinspace \left|\negthinspace \left|{ \tilde{u} - u } \right|\negthinspace \right|\negthinspace \right|$.

■

We now complete the proof of the equivalence result by proving Lemma 3.7.

Proof of Lemma 3.7.

Let us first consider the case of $x_i$ being an interior node. Given $v \in W({\omega _i})$, we construct ${\Pi _i}v \in {\mathcal{P}}_0^2({\omega _i})$ such that (i)–(ii) are satisfied. Let $q_i$, $q_S$ be the canonical piecewise quadratic basis functions corresponding to the center node of ${\omega _i}$, and to the midpoint of $S \subset {\Gamma _i}$, respectively. Let ${\Pi _i}v = \sum _{S} \alpha _S q_S + \beta q_i$, where the sum ranges over all $S \subset {\Gamma _i}$. Then ${\Pi _i}v \in {\mathcal{P}}^2({\omega _i})$ for any choice of $\alpha _S$ and $\beta$. To fulfill (i) we ask

$$\begin{equation*} \alpha _S = \frac{ \int _S v \phi _i - \beta \int _S q_i \phi _i }{ \int _S q_S \phi _i }, \qquad \text{for all } S \subset {\Gamma _i}. \end{equation*}$$

On the one hand, ${\Pi _i}v \in {\mathcal{P}}_0^2({\omega _i})$ if $\int _{\omega _i}{\Pi _i}v \,\phi _i = 0$, and on the other hand

$$\begin{align*} \int _{\omega _i}{\Pi _i}v \, \phi _i &= \int _{\omega _i}\beta q_i \phi _i + \sum _S \alpha _S \int _{\omega _i}q_S \phi _i \\ &= \beta \left( \int _{\omega _i}q_i \phi _i - \sum _S \frac{\int _S q_i \phi _i \int _{\omega _i}q_S \phi _i}{ \int _S q_S \phi _i } \right) + \sum _S \frac{\int _S v \phi _i\int _{\omega _i}q_S \phi _i}{ \int _S q_S \phi _i }. \end{align*}$$

Thus, ${\Pi _i}v \in {\mathcal{P}}_0^2({\omega _i})$ if

$$\begin{equation} \beta \left( \int _{\omega _i}q_i \phi _i - \sum _S \frac{\int _S q_i \phi _i \int _{\omega _i}q_S \phi _i}{ \int _S q_S \phi _i } \right) = -\sum _S \frac{\int _S v \phi _i\int _{\omega _i}q_S \phi _i}{ \int _S q_S \phi _i }. {\tag{3.5}} \end{equation}$$

Denoting with ${\omega _S}$ the support of $q_S$, that is, the union of the two elements sharing $S$, we have

$$\begin{equation*} \int _S q_S \phi _i = \frac{|S|}{3},\quad \int _{\omega _i}q_i \phi _i = \frac{|{\omega _i}|}{30}, \quad \int _{\omega _i}q_S \phi _i = \frac{2|{\omega _S}|}{15}, \quad \text{and}\quad \int _S q_i \phi _i = \frac{|S|}{6}, \end{equation*}$$

which imply that

$$\begin{equation} \int _{\omega _i}q_i \phi _i - \sum _S \frac{\int _S q_i \phi _i \int _{\omega _i}q_S \phi _i}{ \int _S q_S \phi _i } = \frac{|{\omega _i}|}{30} - \frac{\sum _S|{\omega _S}|}{15} = - \frac{|{\omega _i}|}{10}. {\tag{3.6}} \end{equation}$$

Finally, ${\Pi _i}v \in {\mathcal{P}}_0^2({\omega _i})$ satisfies (i) by choosing

$$\begin{equation} \alpha _S = \frac{3}{|S|} \left( \int _S v \phi _i - \beta \frac{|S|}{6}\right), \qquad \beta = \frac{4}{|{\omega _i}|} \sum _S \frac{|{\omega _S}|}{|S|} \int _S v \phi _i . {\tag{3.7}} \end{equation}$$

It remains to prove that choosing $\alpha _S$ and $\beta$ as in (3.7) we obtain (ii). To this end, observe that (3.7), the trace theorem in $L^1$ and the Poincaré inequality for functions with weighted mean value zero given in Proposition 2.4, yield

$$\begin{equation*} \int _S |v\phi _i| \le C h_i \|\nabla v\|_{\phi _i} \quad \Rightarrow \quad |\beta | \le C \left\| \nabla v \right\|_{\phi _i}; \end{equation*}$$

the constant $C$ depends only on the minimum angle of ${\mathcal{T}}_h$. Therefore

$$\begin{equation} \left\| \nabla (\beta q_i) \right\|_{\phi _i} \le C \left\| \nabla v \right\|_{\phi _i} \left\| \nabla q_i \right\|_{\omega _i}\le C \left\| \nabla v \right\|_{\phi _i} . {\tag{3.8}} \end{equation}$$

Analogously, from (3.7) we obtain

$$\begin{equation*} |\alpha _S| \le \frac{3}{|S|} \left\| v \phi _i \right\|_{L^1(S)} + \frac{3}{|S|} C \left\| \nabla v \right\|_{\phi _i} \frac{|S|}{6} \le C \left\| \nabla v \right\|_{\phi _i}, \end{equation*}$$

whence

$$\begin{equation} \left\| \nabla (\alpha _S q_S) \right\|_{\omega _i}\le C \left\| \nabla v \right\|_{\phi _i}. {\tag{3.9}} \end{equation}$$

The claim (ii) then follows from (3.8) and (3.9).

Since $q_i \notin {\mathcal{P}}_0^2({\omega _i})$ for a boundary node $x_i$, the above arguments still apply upon taking $\beta = 0$.

■

4. Convergence

We start this section with some ingredients that will be essential for the convergence of an adaptive algorithm. The proof of Proposition 4.1 is given in §5.

Proposition 4.1 (Error reduction).

Let ${\mathcal{T}}_H$ be a triangulation of $\Omega$ and let $\theta \in (0,1)$. Let $\hat{{\mathcal{N}}_H}$ be a subset of ${\mathcal{N}}_H$ satisfying

$$\begin{equation} \bigg (\sum _{ i \in \hat{{\mathcal{N}}_H}} {\mathcal{ E}_i}^2\bigg )^{1/2} \ge \theta \, \mathcal{ E}_H. {\tag{4.1}} \end{equation}$$

Let ${\mathcal{T}}_h$ be a refinement of ${\mathcal{T}}_H$ such that:

$$\begin{equation} \begin{gathered} \text{for every element $T\in {\mathcal{T}}_H$ having a node in $\hat{{\mathcal{N}}_H}$,}\\ \text{each of the four similar triangles resulting from joining the midside}\\ \text{points of $T$ can be expressed as a union of triangles of ${\mathcal{T}}_h$.} \end{gathered} {\tag{4.2}} \end{equation}$$ Then there exist constants $\alpha \in (0,1)$, $C_{\operatorname {osc}}$, depending on the minimum angle of ${\mathcal{T}}_H$ and ${\mathcal{T}}_h$, $\theta$, $C_a$ and $c_a$ such that

$$\begin{equation*} \left|\negthinspace \left|\negthinspace \left|{u-u_h} \right|\negthinspace \right|\negthinspace \right| \le \alpha \left|\negthinspace \left|\negthinspace \left|{u-u_H} \right|\negthinspace \right|\negthinspace \right| + C_{\operatorname {osc}} \operatorname {osc}(f,g;{{\mathcal{N}}_H}), \end{equation*}$$

where $\operatorname {osc}(f,g;{{\mathcal{N}}_H})$ is as defined in Theorem 3.6.

Remark 4.2.

The marking strategy (4.1) ensures that we choose sufficiently many stars such that their contributions ${\mathcal{ E}_i}$ constitute a fixed proportion of the total error estimator $\mathcal{ E}_H$. This strategy was first introduced by Dörfler 6, and was also used in 11.

Remark 4.3.

Using newest-vertex bisectioning, the crucial property (4.2) can be achieved via three bisection steps as depicted in Figure 2: an element $T\in {\mathcal{T}}_H$ is first bisected twice, thereby ensuring new nodes of ${\mathcal{T}}_h$ at the midside points of $T$; secondly, bisecting once more the two grandchildren with index 1 yields (4.2). This procedure guarantees that ${\mathcal{T}}_h$ and ${\mathcal{T}}_H$ are nested, and that their minimum angle is bounded below by a positive number depending only on the initial mesh. Local regular (red) refinement, combined with the green closure procedure, would also enforce (4.2) but not that ${\mathcal{T}}_h$ and ${\mathcal{T}}_H$ are nested. Therefore, the above bisection procedure is the method of choice, and is used in the numerical experiments of §8.

Since data oscillation is present in the error bound of two consecutive solutions, a convergent algorithm must also ensure a reduction of data oscillation. The following proposition states one possible way of doing this; its proof is given in §6.

Proposition 4.4 (Oscillation reduction).

Let ${\mathcal{T}}_H$ be a triangulation of $\Omega$ and let $\theta _f, \theta _g \in (0,1)$. Assume that $\hat{{\mathcal{N}}_H}\subset {\mathcal{N}}_H$ satisfies

$$\begin{equation} \operatorname {osc}(f;{\hat{{\mathcal{N}}_H}}) \ge \theta _f \operatorname {osc}(f;{{\mathcal{N}}_H}). {\tag{4.3}} \end{equation}$$

If ${\mathcal{T}}_h$ is a refinement of ${\mathcal{T}}_H$ such that all the elements having a node in $\hat{{\mathcal{N}}_H}$ are at least bisected twice, then there exists $\beta _f \in (0,1)$ depending only on the minimum angle of ${\mathcal{T}}_H$ and ${\mathcal{T}}_h$ and $\theta _f$ such that

$$\begin{equation*} \operatorname {osc}(f;{{\mathcal{N}}_h}) \le \beta _f \operatorname {osc}(f;{{\mathcal{N}}_H}). \end{equation*}$$

Analogously, if

$$\begin{equation} \operatorname {osc}_\partial (g;{\hat{{\mathcal{N}}_H}}) \ge \theta _g \operatorname {osc}_\partial (g;{{\mathcal{N}}_H}) {\tag{4.4}} \end{equation}$$

and ${\mathcal{T}}_h$ is as above, there exists $\beta _g \in (0,1)$ depending only on the minimum angle of ${\mathcal{T}}_H$ and ${\mathcal{T}}_h$ and $\theta _g$ such that

$$\begin{equation*} \operatorname {osc}_\partial (g;{{\mathcal{N}}_h}) \le \beta _g \operatorname {osc}_\partial (g;{{\mathcal{N}}_H}). \end{equation*}$$

Inspired in these two results we now propose the adaptive algorithm and prove its convergence.

Convergent Algorithm C.

Given parameters $0<\theta , \theta _f, \theta _g, \nu _f, \nu _g <1$:

(1)

Choose any initial mesh ${\mathcal{T}}_0$ where $A$ is piecewise constant.

(2)

Compute the discrete solution $u_0$ on ${\mathcal{T}}_0$.

(3)

Let $k=0$.

(4)

Compute the local indicators ${\mathcal{ E}_i}$, $\operatorname {osc}(f;{{\mathcal{N}}_k})$, $\operatorname {osc}_\partial (g;{{\mathcal{N}}_k})$.

(5)

Construct a subset $\hat{{\mathcal{N}}_k}$ of ${\mathcal{N}}_k$ such that (4.1) holds.

(6)

If $\operatorname {osc}(f;{{\mathcal{N}}_k}) > \operatorname {osc}(f;{{\mathcal{N}}_0}) \nu _f^{k}$, enlarge $\hat{{\mathcal{N}}_k}$ to satisfy (4.3).

(7)

If $\operatorname {osc}_\partial (g;{{\mathcal{N}}_k}) > \operatorname {osc}_\partial (g;{{\mathcal{N}}_0}) \nu _g^{k}$, enlarge $\hat{{\mathcal{N}}_k}$ to satisfy (4.4).

(8)

Refine ${\mathcal{T}}_k$ into ${\mathcal{T}}_{k+1}$ according to (4.2).

(9)

Compute the discrete solution $u_{k+1}$ on ${\mathcal{T}}_{k+1}$.

(10)

Set $k \leftarrow k + 1$.

(11)

Go to Step 4.

Remark 4.5.

Step 8 must use a refinement procedure that guarantees the minimum angle property. In the experiments of §8 we use the newest-vertex bisection.

Remark 4.6.

In contrast with 11, we no longer impose a data oscillation reduction rate in every step, but rather a decay bounded by exponentials $\nu _f^k$ and $\nu _g^k$. This apparently minor change influences dramatically the role of data oscillation in Algorithm C, as is corroborated experimentally in §8.

Theorem 4.7 (Convergence).

Let $u_k$ be the sequence of finite element solutions produced by Algorithm C. Then, there exist positive constants $C_0$ and $\gamma <1$, depending only on given data and the initial grid, such that

$$\begin{equation} \left|\negthinspace \left|\negthinspace \left|{u - u_k} \right|\negthinspace \right|\negthinspace \right|_\Omega \le C_0 \,\gamma ^k. {\tag{4.5}} \end{equation}$$

Proof.

We claim the existence of $\mu < 1$ and $B_0$, depending only on data and the initial grid, such that

$$\begin{equation} \operatorname {osc}(f,g;{{\mathcal{N}}_k}) \le B_0 \mu ^k, \qquad \forall\nobreakspace k\ge 0. {\tag{4.6}} \end{equation}$$

We postpone the proof of (4.6) and show first how to exploit it. Let $e_k$ be the energy norm of the error in step $k$. Then, by Proposition 4.1,

$$\begin{equation*} {e_{k+1}} \le \alpha {e_k} + C_{\operatorname {osc}} \operatorname {osc}(f,g;{\mathcal{N}}_k)\le \alpha e_k + C_{\operatorname {osc}} B_0 \mu ^{k} \le \delta \left( e_k + C \delta ^{k-1} \right), \end{equation*}$$

where $\delta := \max \{ \alpha , \mu \} < 1$. By induction we obtain

$$\begin{equation*} {e_{k}} \le \delta ^{k} e_0 + k C \delta ^{k-1} , \end{equation*}$$

which implies (4.5) if we choose $\gamma = \sqrt {\delta }$.

To complete the proof it remains to prove (4.6). Regarding $\operatorname {osc}(f;{{\mathcal{N}}_k})$ we claim that

$$\begin{equation*} \operatorname {osc}(f;{{\mathcal{N}}_k}) \le \frac{\operatorname {osc}(f;{{\mathcal{N}}_0})}{\nu _f} \max (\nu _f,\beta _f)^k \qquad \forall\nobreakspace k\ge 0, \end{equation*}$$

and resort to an induction argument. Such a bound holds trivially for $k=0$. We assume it holds for $k$. Then we have either

$$\begin{equation*} \text{(i)}\quad \operatorname {osc}(f;{{\mathcal{N}}_k}) > \operatorname {osc}(f;{{\mathcal{N}}_0}) \nu _f^{k} \qquad \text{or}\qquad \text{(ii)}\quad \operatorname {osc}(f;{{\mathcal{N}}_k}) \le \operatorname {osc}(f;{{\mathcal{N}}_0}) \nu _f^{k}. \end{equation*}$$

In case (i), we see from Step 6 of Algorithm C and Proposition 4.4 that

$$\begin{equation*} \operatorname {osc}(f;{{\mathcal{N}}_{k+1}}) \le \beta _f \operatorname {osc}(f;{{\mathcal{N}}_k}) \le \frac{\operatorname {osc}(f;{{\mathcal{N}}_0})}{\nu _f} \max (\nu _f,\beta _f)^{k+1} . \end{equation*}$$

On the other hand, exploiting that ${\mathcal{T}}_{k+1}$ is a refinement of ${\mathcal{T}}_k$, and thus the oscillation must not increase (see Lemma 6.1 below), we can handle case (ii) as follows:

$$\begin{equation*} \operatorname {osc}(f;{{\mathcal{N}}_{k+1}}) \le \operatorname {osc}(f;{{\mathcal{N}}_k}) \le \operatorname {osc}(f;{{\mathcal{N}}_0}) \nu _f^{k} \le \frac{\operatorname {osc}(f;{{\mathcal{N}}_0})}{\nu _f} \max (\nu _f,\beta _f)^{k+1}. \end{equation*}$$

The same analysis for $\operatorname {osc}_\partial (g;{{\mathcal{N}}_k})$ leads to an analogous result for boundary data oscillation, thereby concluding the proof of (4.6).

■

Remark 4.8.

This result is also true for the usual residual type error estimators. In fact, we will see in §5 that the only assumption for Proposition 4.1 to be true is that the estimators provide a global upper bound of the error and a local lower bound for the difference between two consecutive solutions (e.g., Lemma 5.2).

5. Error reduction

In this section we are going to prove Proposition 4.1. We start with the following result that in the case of $g=g_H\nobreakspace (=g_h)$ reduces to Pythagoras’ theorem.

Lemma 5.1.

Let ${\mathcal{T}}_H$ and ${\mathcal{T}}_h$ be two triangulations satisfying ${\mathbb{V}}_{\!H}\subset {\mathbb{V}}_{\!h}$. Then, there exists a constant $C$, depending on $\Omega$, the minimum angle of ${\mathcal{T}}_H$ and ${\mathcal{T}}_h$, $C_a$ and $c_a$ such that the following relation holds for any $\delta \in (0,1)$

$$\begin{equation} (1-\delta )\left|\negthinspace \left|\negthinspace \left|{u - u_h} \right|\negthinspace \right|\negthinspace \right|_\Omega ^2 \le \left|\negthinspace \left|\negthinspace \left|{u - u_H} \right|\negthinspace \right|\negthinspace \right|_\Omega ^2 - \left|\negthinspace \left|\negthinspace \left|{u_H- u_h} \right|\negthinspace \right|\negthinspace \right|_\Omega ^2 + \frac{C}{\delta } \operatorname {osc}_\partial (g;{{\mathcal{N}}_H})^2. {\tag{5.1}} \end{equation}$$

Proof.

The bilinearity and symmetry of $\left(\hskip -1pt\left({\cdot }\, , \,{\cdot }\right)\hskip -1pt\right)$ imply that

$$\begin{equation*} \left|\negthinspace \left|\negthinspace \left|{u - u_H} \right|\negthinspace \right|\negthinspace \right|^2 = \left|\negthinspace \left|\negthinspace \left|{u - u_h} \right|\negthinspace \right|\negthinspace \right|^2 + \left|\negthinspace \left|\negthinspace \left|{u_h- u_H} \right|\negthinspace \right|\negthinspace \right|^2 + 2 \left(\hskip -1pt\left({u - u_h}\, , \,{u_h- u_H}\right)\hskip -1pt\right). \end{equation*}$$

Let $\bar{g}_h$ (resp. $\bar{g}_H$) denote the extension of $g_h$ (resp. $g_H$) to ${\mathbb{V}}_{\!h}$ defined as 0 at all the interior nodes of ${\mathcal{T}}_h$. Then $\bar{g}_h - \bar{g}_H \in {\mathbb{V}}_{\!h}$ and $u_h - u_H - (\bar{g}_h - \bar{g}_H) \in {\mathbb{V}}_{\!h}^0$, which imply

$$\begin{align*} \left|\negthinspace \left|\negthinspace \left|{u - u_H} \right|\negthinspace \right|\negthinspace \right|^2 &= \left|\negthinspace \left|\negthinspace \left|{u - u_h} \right|\negthinspace \right|\negthinspace \right|^2 + \left|\negthinspace \left|\negthinspace \left|{u_h- u_H} \right|\negthinspace \right|\negthinspace \right|^2 + 2 \left(\hskip -1pt\left({u - u_h}\, , \,{\bar{g}_h - \bar{g}_H}\right)\hskip -1pt\right) \\ &\ge \left|\negthinspace \left|\negthinspace \left|{u - u_h} \right|\negthinspace \right|\negthinspace \right|^2 + \left|\negthinspace \left|\negthinspace \left|{u_h- u_H} \right|\negthinspace \right|\negthinspace \right|^2 - \delta \left|\negthinspace \left|\negthinspace \left|{u - u_h} \right|\negthinspace \right|\negthinspace \right|^2 - \frac{1}{\delta } \left|\negthinspace \left|\negthinspace \left|{\bar{g}_h - \bar{g}_H} \right|\negthinspace \right|\negthinspace \right|_\Omega ^2 \end{align*}$$

for any $\delta >0$. Let ${\mathcal{N}}_h^\partial$ denote the set of indices of nodes $x_i \in \partial \Omega$ of ${\mathcal{N}}_h$ which are not nodes of ${\mathcal{N}}_H$. Since $(g_h - g_H) (x_i) \neq 0$ only for $i \in {\mathcal{N}}_h^\partial$, it is easily seen that

$$\begin{align*} \|\nabla (\bar{g}_h - \bar{g}_H)\|_\Omega ^2 &\le C \sum _{i\in {\mathcal{N}}_h^\partial } | (g_h-g_H) (x_i) |^2 \le C \sum _{i\in {\mathcal{N}}_h^\partial } |S_i| \| \partial _\tau (g_h-g_H)\|_{L^2(S_i)}^2 \\ &\le C (\operatorname {osc}_\partial (g;{{\mathcal{N}}_H})^2 + \operatorname {osc}_\partial (g;{{\mathcal{N}}_h})^2), \end{align*}$$

where $C$ depends only on the minimum angle. In view of $\operatorname {osc}_\partial (g;{{\mathcal{N}}_h}) \le \operatorname {osc}_\partial (g;{{\mathcal{N}}_H})$, which is proved in Lemma 6.4 below, we deduce the claim.

■

The following lemma is crucial for proving Proposition 4.1.

Lemma 5.2.

Let ${\mathcal{T}}_h$ be a refinement of ${\mathcal{T}}_H$ satisfying 4.2. If $i \in \hat{{\mathcal{N}}_H}$, then

$$\begin{equation} {\mathcal{ E}_i}^2 \le C_4 \left|\negthinspace \left|\negthinspace \left|{ u_h- u_H} \right|\negthinspace \right|\negthinspace \right|_{\omega _i}^2 + C_5 \bigl \| h_i (f-f_i) \bigr \|_{\phi _i}^2, {\tag{5.2}} \end{equation}$$

where $C_4$ and $C_5$ depend only on the minimum angle, $C_a$ and $c_a$.

Remark 5.3.

Property (4.2) is essential for error reduction. In fact, if we solve $-\Delta u = 1$ in $\Omega =(0,1)^2$, $u = 0$ on $\partial \Omega$, the solutions obtained in the two grids of Figure 3 coincide: $u_H=u_h=\frac{1}{12}\phi _1$. Hence, no error reduction is obtained in the interior star of ${\mathcal{T}}_H$ even with no data oscillation. For details about the important role of data oscillation, see 11.

To be able to prove Lemma 5.2 we need to introduce a projection operator analogous to that of Lemma 3.7. Since we want to compare two consecutive solutions $u_H$ and $u_h$, we now need the target space of the projection to be contained in the space of piecewise linear functions on the finer mesh. To this end, we first introduce the following definition.

Definition 5.4.

Let ${\mathcal{T}}_h$ be a refinement of ${\mathcal{T}}_H$ satisfying (4.2) and let $i\in \hat{{\mathcal{N}}_H}$. If the center node $x_i$ of the star ${\omega _i}$ is a boundary node, we define ${\mathcal{P}}_0^1({\omega _i})$ as the space of functions in ${\mathbb{V}}_{\!h}^0$ vanishing outside ${\omega _i}$. In case $x_i$ is an interior node, ${\mathcal{P}}_0^1({\omega _i})$ will denote the space of functions $v \in {\mathbb{V}}_{\!h}^0$ vanishing outside ${\omega _i}$ that also satisfy $\int _{\omega _i}v = 0$ (note the absence of the weight $\phi _i$ in the last integral).

Lemma 5.5.

Let $H^1(\omega _i)$ be the space of functions in $H^1({\omega _i})$ with mean value zero provided $x_i$ is an interior node, and let $H^1(\omega _i)= H_0^1({\omega _i})$ for boundary nodes. Then for each node $x_i$, $i \in {\mathcal{N}}_H$, there exists an operator $P_i:H^1(\omega _i)\to {\mathcal{P}}_0^1({\omega _i})$ such that for any $v \in H^1(\omega _i)$, the following conditions hold:

(i)

$\int _S v - P_i v = 0$, for all $S \subset {\Gamma _i}$,

(ii)

$\| \nabla P_i v \|_{{\omega _i}} \le C^* \| \nabla v \|_{{\omega _i}}$,

(iii)

$\| P_i v \|_{\phi _i^{-1}} \le C^* h_i \| \nabla v \|_{{\omega _i}}$,

where the constant $C^*$ depends only on the minimum angle of ${\mathcal{T}}_H$ and ${\mathcal{T}}_h$.

The proof of this lemma is similar to that of Lemma 3.7 and we postpone it to the end of this section.

Proof of Lemma 5.2.

Let ${\mathcal{ E}_i}$ and $\eta _i$ be as in Definition 2.5. Then, by Lemma 5.5(i)

$$\begin{align*} {\mathcal{ E}_i}^2 &= \left|\negthinspace \left|\negthinspace \left|{{\eta _i}} \right|\negthinspace \right|\negthinspace \right|_{\phi _i}^2 = \int _{\omega _i}A \nabla \eta _i \nabla \eta _i \,\phi _i = \int _{\omega _i}f \eta _i \phi _i + \int _{\Gamma _i}J \eta _i \phi _i \\ &= \int _{\omega _i}f P_i(\eta _i \phi _i) + \int _{\Gamma _i}J P_i(\eta _i \phi _i) + \int _{\omega _i}(f-f_i) (\eta _i\phi _i - P_i(\eta _i \phi _i)), \end{align*}$$

where $\displaystyle f_i = \frac{\int _{\omega _i}f \phi _i }{ \int _{\omega _i}\phi _i}$ for interior nodes and $f_i = 0$ otherwise. Since $P_i (\eta _i \phi _i) \in {\mathbb{V}}_{\!h}^0$,

$$\begin{align*} {\mathcal{ E}_i}^2 &= \int _\Omega A \nabla (u_h-u_H) \nabla P_i (\eta _i \, \phi _i) + \int _{\omega _i}(f-f_i) (\eta _i\phi _i - P_i(\eta _i \phi _i)), \\ &\le C_{\!A} \left|\negthinspace \left|\negthinspace \left|{u_h- u_H} \right|\negthinspace \right|\negthinspace \right|_{\omega _i}\left\| \nabla P_i (\eta _i\phi _i)\right\|_{\omega _i}+ \| f - f_i \|_{\phi _i} \left( \| \eta _i \|_{{\omega _i}} + \| P_i (\eta _i\phi _i) \|_{\phi _i^{-1}} \right). \end{align*}$$

On the one hand, by Lemma 5.5(ii) and (iii) we have, respectively,

$$\begin{equation*} \left\| \nabla P_i (\eta _i\phi _i)\right\|_{\omega _i}\le C^* \left\| \nabla (\eta _i\phi _i) \right\|_{\omega _i}\quad \text{and}\quad \left\| P_i(\eta _i\phi _i) \right\|_{\phi _i^{-1}} \le C^* h_i \left\| \nabla (\eta _i\phi _i) \right\|_{\omega _i}. \end{equation*}$$

On the other hand, by Proposition 2.4

$$\begin{equation*} \left\| \nabla (\eta _i\phi _i) \right\|_{\omega _i}\le \left\| \nabla \eta _i \, \phi _i \right\|_{\omega _i}+ \left\| \eta _i\, \nabla \phi _i \right\|_{\omega _i}\le \left\| \nabla \eta _i \right\|_{ \phi _i} + \frac{C}{h_i} \left\| \eta _i \right\|_{\omega _i}\le (1 + C) \left\| \nabla \eta _i \right\|_{ \phi _i}. \end{equation*}$$

Hence,

$$\begin{equation*} {\mathcal{ E}_i}^2 \le \left(C_4 \left|\negthinspace \left|\negthinspace \left|{u_h- u_H} \right|\negthinspace \right|\negthinspace \right|_{\omega _i}+ C_5 h_i \| f - f_i \|_{\phi _i} \right) \left\| \nabla \eta _i \right\|_{\phi _i}, \end{equation*}$$

which immediately implies the claim.

■

Before proving the main result we need the following corollary of Lemma 5.2.

Corollary 5.6.

Suppose the assumptions of Proposition 4.1 hold. Then, we have the following global lower bound for the error reduction

$$\begin{equation*} \left|\negthinspace \left|\negthinspace \left|{u_h- u_H} \right|\negthinspace \right|\negthinspace \right|^2 \ge \frac{\theta ^2}{C_6} \left|\negthinspace \left|\negthinspace \left|{u - u_H} \right|\negthinspace \right|\negthinspace \right|^2 - C_7 \operatorname {osc}(f,g;{{\mathcal{N}}_H})^2, \end{equation*}$$

where the constants $C_6$ and $C_7$ depend only on the minimum angle of ${\mathcal{T}}_H$ and ${\mathcal{T}}_h$, $C_a$ and $c_a$.

Proof.

By Lemma 5.2 and 4.1 we have

$$\begin{align*} \theta ^2 \mathcal{ E}_H^2 &\le \sum _{i \in \hat{{\mathcal{N}}_H}} {\mathcal{ E}_i}^2 \le C_4 \sum _{i \in \hat{{\mathcal{N}}_H}} \left|\negthinspace \left|\negthinspace \left|{u_h- u_H} \right|\negthinspace \right|\negthinspace \right|_{\omega _i}^2 + C_5 \sum _{i \in \hat{{\mathcal{N}}_H}} \bigl \| h_i (f- f_i) \bigr \|_{\phi _i}^2 \\ &\le 3 \, C_4 \left|\negthinspace \left|\negthinspace \left|{u_h- u_H} \right|\negthinspace \right|\negthinspace \right|_\Omega ^2 + C_5 \operatorname {osc}(f;{{\mathcal{N}}_H})^2. \end{align*}$$

Hence,

$$\begin{align*} \left|\negthinspace \left|\negthinspace \left|{u_h- u_H} \right|\negthinspace \right|\negthinspace \right|_\Omega ^2 &\ge \frac{\theta ^2}{3 \, C_4} \mathcal{ E}_H^2 - \frac{C_5}{3 \, C_4} \operatorname {osc}(f;{{\mathcal{N}}_H})^2, \end{align*}$$

and the assertion follows from Theorem 3.6.

■

Proof of Proposition 4.1.

In view of Lemma 5.1 and Corollary 5.6, for any $\delta \in (0,1)$, we have

$$\begin{equation*} \begin{split} (1-\delta ) \left|\negthinspace \left|\negthinspace \left|{u - u_h} \right|\negthinspace \right|\negthinspace \right|_\Omega ^2 &\le \left|\negthinspace \left|\negthinspace \left|{u - u_H} \right|\negthinspace \right|\negthinspace \right|_\Omega ^2 - \left|\negthinspace \left|\negthinspace \left|{u_h- u_H} \right|\negthinspace \right|\negthinspace \right|_\Omega ^2 + \frac{C}{\delta } \operatorname {osc}_\partial (g;{{\mathcal{N}}_H})^2 \\ & \le \left|\negthinspace \left|\negthinspace \left|{u - u_H} \right|\negthinspace \right|\negthinspace \right|_\Omega ^2 \left( 1 - \frac{\theta ^2}{C_6} \right) + C_\delta \operatorname {osc}(f,g;{{\mathcal{N}}_H})^2. \end{split} \end{equation*}$$

Therefore, if $\delta$ is sufficiently small such that $\alpha ^2 := \displaystyle \frac{1 - \frac{\theta ^2}{C_6}}{1-\delta } < 1$, we have

$$\begin{equation*} \left|\negthinspace \left|\negthinspace \left|{u - u_h} \right|\negthinspace \right|\negthinspace \right|_\Omega \le \alpha \left|\negthinspace \left|\negthinspace \left|{u - u_H} \right|\negthinspace \right|\negthinspace \right|_\Omega + C_{\operatorname {osc}} \operatorname {osc}(f,g;{{\mathcal{N}}_H}). \end{equation*}$$■

We now complete the proof of the error reduction result by proving Lemma 5.5.

Proof of Lemma 5.5.

Let us assume that ${\omega _i}$ is an interior star; the case of a boundary star can be treated similarly and is in fact simpler. By Definition 5.4, for every $S\subset {\Gamma _i}$, the set ${\omega _S^h}$ (resp. ${\omega _i^h}$) depicted in Figure 4 (middle) (resp. right) is the support of the piecewise linear funcion $\varphi _S$ (resp. $\varphi _i$) $\in {\mathbb{V}}_{\!h}^0$ that equals 1 at the midpoint of $S$ (resp. $x_i$).

The proof of (i) and (ii) follows that of Lemma 3.7 upon replacing ${\Pi _i}v$, $q_S$, $q_i$ and $\phi _i$ by $P_i v$, $\varphi _S$, $\varphi _i$ and $1$, respectively. In this vein, it is important to mention that the construction of $P_i v$ is possible because the requirement $\int _{\omega _i}P_i v = 0$ corresponds to the equation (see (3.5) in the proof of Lemma 3.7)

$$\begin{equation*} \beta \left( \int _{\omega _i}\varphi _i - \sum _S \frac{\int _S \varphi _i \, \int _{\omega _i}\varphi _S }{\int _S \varphi _S } \right) = - \sum _S \frac{\int _S v \, \int _{\omega _i}\varphi _S }{\int _S \varphi _S }. \end{equation*}$$

Since

$$\begin{equation*} \int _S \varphi _S = \frac{|S|}{2}, \quad \int _{\omega _i}\varphi _i = \frac{|{\omega _i}|}{12},\quad \int _{\omega _i}\varphi _S = \frac{|{\omega _S}|}{4},\quad \int _S \varphi _i = \frac{|S|}{4}, \end{equation*}$$

and ${\omega _S}$ is the union of the two elements sharing $S$, we deduce the solvability condition

$$\begin{equation*} \int _{\omega _i}\varphi _i - \sum _S \frac{\int _S \varphi _i \, \int _{\omega _i}\varphi _S }{\int _S \varphi _S } = \frac{|{\omega _i}|}{12} - \sum _S \frac{ \frac{|S|}{4} \frac{|{\omega _S}|}{4} }{ \frac{|S|}{2} } = \frac{|{\omega _i}|}{12} - \sum _S \frac{|{\omega _S}|}{8} = -\frac{|{\omega _i}|}{6} \neq 0. \end{equation*}$$

Regarding (iii) it is sufficient to prove

$$\begin{equation} \left\| w \right\|_{\phi _i^{-1}} \le C \, h_i \, \left\| \nabla w \right\|_{\omega _i}, \qquad \forall w \in {\mathcal{P}}_0^1({\omega _i}), {\tag{5.3}} \end{equation}$$

which follows from the equivalence of norms on the finite dimensional space ${\mathcal{P}}_0^1({\omega _i})$ and a scaling argument. The claim then follows by taking $w = P_i v$ in 5.3 and using (ii).

■

Remark 5.7.

If we perform only two steps of the newest-vertex bisectioning, we may not meet condition (4.2). Then a pathological situation could occur whenever the center node of the star ${\omega _i}$ is the newest vertex of all the triangles in ${\omega _i}$. In that case, two bisections would lead to a refinement like the one depicted in Figure 5, which would imply a unique choice of ${\omega _S^h}$, ${\omega _i^h}$, $\varphi _S$ and $\varphi _i$. This choice would give

$$\begin{equation*} \int _S \varphi _S = \frac{|S|}{2}, \quad \int _{\omega _i}\varphi _i = \frac{|{\omega _i}|}{12},\quad \int _{\omega _i}\varphi _S = \frac{|{\omega _S}|}{6},\quad \int _S \varphi _i = \frac{|S|}{4}, \end{equation*}$$

whence

$$\begin{equation*} \int _{\omega _i}\varphi _i - \sum _S \frac{\int _S \varphi _i \, \int _{\omega _i}\varphi _S }{\int _S \varphi _S } = \frac{|{\omega _i}|}{6} - \sum _S \frac{ \frac{|S|}{4} \frac{|{\omega _S}|}{6} }{ \frac{|S|}{2} } = \frac{|{\omega _i}|}{6} - \sum _S \frac{|{\omega _S}|}{12} = 0. \end{equation*}$$

Therefore, for $\int _{\omega _i}P_i v$ to be zero we need $\beta$ to satisfy

$$\begin{equation*} \beta \cdot 0 = - \sum _S \frac{ \int _S v \, \frac{|{\omega _S}|}{4} }{ \frac{|S|}{2} }, \end{equation*}$$

but the right hand side will not vanish in general. For this pathological case not to occur, it would be sufficient to refine only one of the star elements as depicted in Figure 2. Since we would not know which one to choose, and to be fair to everyone and thus politically correct, we decided to refine all star triangles that way.

6. Oscillation reduction

As we can see in the proof of Theorem 4.7, to be able to obtain a bound like $\operatorname {osc}(f,g;{{\mathcal{N}}_k})\le B_0 \mu ^k$ we need two different results. The first one is that the oscillation does not increase by refinement, and the second one is that marking elements according to (4.3) and (4.4) implies an oscillation reduction by a fixed factor. We now present four lemmas, which show these two results for both $f$ and $g$.

Lemma 6.1.

Let ${\mathcal{T}}_H$ and ${\mathcal{T}}_h$ be two triangulations satisfying ${\mathbb{V}}_{\!H}\subset {\mathbb{V}}_{\!h}$. Then

$$\begin{equation*} \operatorname {osc}(f;{{\mathcal{N}}_h}) \le \operatorname {osc}(f;{{\mathcal{N}}_H}). \end{equation*}$$

Proof.

Let us indicate with $\{ \phi _j \}_{j\in {\mathcal{N}}_H}$ the canonical basis of ${\mathbb{V}}_{\!H}$, and with $\{ \varphi _i \}_{i\in {\mathcal{N}}_h}$, $\{x_i\}_{i\in {\mathcal{N}}_h}$ the canonical basis of ${\mathbb{V}}_{\!h}$ and the set of nodes of ${\mathcal{T}}_h$, respectively. Then, if we define $\alpha _{ji} := \phi _j(x_i)$, we have that

$$\begin{equation*} \phi _j = \sum _i \alpha _{ji} \, \varphi _i, \qquad \alpha _{ji} \ge 0, \quad \forall i,j, \qquad \text{and}\quad \sum _{j} \alpha _{ji} = 1, \end{equation*}$$

where hereafter the sums on $i$ and $j$ range over ${\mathcal{N}}_h$ and ${\mathcal{N}}_H$, respectively. For the integrals that will appear in the rest of the proof, the domain of integration will be the support of $\varphi _i$ or $\phi _j$ depending on the context, and we will omit it. Since $f_i^h = \int f \varphi _i / \int \varphi _i$ is optimal for interior nodes, we can write

$$\begin{equation*} \int |f - f_i^h |^2 \varphi _i \le \int |f - \sum _j \alpha _{ji}\, f_j^H |^2 \varphi _i. \end{equation*}$$

This inequality is also valid for boundary nodes because $\alpha _{ji} \neq 0$ only for boundary $j$’s, for which also $f_j^H = 0$.

Since $\sum _j \alpha _{ji} = 1$ and $\alpha _{ij}\ge 0$, by the convexity of the function $s \to |s|^2$ in $\mathbb{R}$, we obtain

$$\begin{equation} \begin{aligned} \operatorname {osc}(f;{{\mathcal{N}}_h})^2 &\le \sum _{i} h_i^2 \int | f - \sum _j \alpha _{ji}\, f_j^H |^2 \varphi _i \\ &\le \sum _i h_i^2 \int \sum _j \alpha _{ji} | f - f_j^H |^2 \varphi _i = \sum _j \int |f - f_j^H|^2 \sum _i h_i^2 \alpha _{ji} \, \varphi _i. \end{aligned} {\tag{6.1}} \end{equation}$$

Now, the sum on $i$ is taken over all $i$’s such that $\operatorname {supp}(\varphi _i) \subset \operatorname {supp}(\phi _j)$ since $\alpha _{ji}=0$ otherwise, and consequently for $h_i \le H_j$, thereby giving

$$\begin{align*} \operatorname {osc}(f;{{\mathcal{N}}_h})^2 &\le \sum _j H_j^2 \int |f - f_j^H|^2 \sum _i \alpha _{ji}\,\varphi _i \\ & = \sum _j H_j^2 \int |f - f_j^H|^2 \phi _j = \operatorname {osc}(f;{{\mathcal{N}}_H})^2. \end{align*}$$■

Lemma 6.2 (Reduction of the oscillation of $f$).

Let $0<r<1$ be the reduction factor of element size associated with one refinement step, i.e., $r$ is the smallest number such that $\operatorname {diam}(T') \leq r\,\operatorname {diam}(T)$ for all subelements $T'$ of a refined element $T$. Given $0 < \theta _f < 1$, let $\hat{{\mathcal{N}}_H}\subset {\mathcal{N}}_H$ satisfy

$$\begin{equation} \operatorname {osc}(f;{\hat{{\mathcal{N}}_H}}) \ge \theta _f \operatorname {osc}(f;{{\mathcal{N}}_H}). \tag{6.2} \end{equation}$$

Let ${\mathcal{T}}_h$ be a refinement of ${\mathcal{T}}_H$ such that all the elements having a node in $\hat{{\mathcal{N}}_H}$ are refined. Then, for $\beta _f := (1 - (1-r^2)\theta _f^2)^{1/2} < 1$, we have

$$\begin{equation} \operatorname {osc}(f;{{\mathcal{N}}_h}) \le \beta _f \operatorname {osc}(f;{{\mathcal{N}}_H}). \tag{6.3} \end{equation}$$

Remark 6.3.

Using refinement by bisection, we need at least two bisections to guarantee $r < 1$. We point out that $r$ depends only on the minimal angle of the initial mesh, and that $r\le 1/2$ provided (4.2) is enforced.

Proof of Lemma 6.2.

If $j \in \hat{{\mathcal{N}}_H}$, then $h_i \le r H_j$ for all stars of ${\mathcal{T}}_h$ contained in ${\omega _j}$. Using (6.1) we arrive at

$$\begin{align*} \operatorname {osc}(f;{{\mathcal{N}}_h}) &\le r^2 \sum _{j \in \hat{{\mathcal{N}}_H}} H_j^2 \int |f - f_j^H|^2 \phi _j + \sum _{j \in {\mathcal{N}}_H\setminus \hat{{\mathcal{N}}_H}} H_j^2 \int |f - f_j^H|^2 \phi _j \\ &= (r^2 - 1) \sum _{j \in \hat{{\mathcal{N}}_H}} H_j^2 \int |f - f_j^H|^2 \phi _j + \operatorname {osc}(f;{{\mathcal{N}}_H})^2 \\ &\le [(r^2 - 1)\theta _f^2 + 1] \operatorname {osc}(f;{{\mathcal{N}}_H})^2. \end{align*}$$■

Lemma 6.4.

Let ${\mathcal{T}}_H$ and ${\mathcal{T}}_h$ be two triangulations satisfying ${\mathbb{V}}_{\!H}\subset {\mathbb{V}}_{\!h}$. Then

$$\begin{equation*} \operatorname {osc}_\partial (g;{{\mathcal{N}}_h}) \le \operatorname {osc}_\partial (g;{{\mathcal{N}}_H}). \end{equation*}$$

Proof.

The claim follows from the fact that ${\mathbb{V}}_{\!H}\subset {\mathbb{V}}_{\!h}$, and the piecewise linear interpolant $g_h$ over ${\mathcal{T}}_h$ coincides with the local projection of $g$ into ${\mathbb{V}}_{\!h}$ in the $H^1(\partial \Omega )$ seminorm, whence

$$\begin{equation*} \|\partial _\tau (g-g_h)\|_{L^2(S_i)} \le \|\partial _\tau (g-g_H)\|_{L^2(S_i)} \qquad \forall\nobreakspace i\in {\mathcal{N}}_h, \quad x_i \in \partial \Omega . \end{equation*}$$■

The following lemma is the counterpart of Lemma 6.2 for the oscillation of the boundary data. Its proof is very similar to that of Lemma 6.2 and is thus omitted.

Lemma 6.5 (Reduction of the oscillation of $g$).

Let $0<r <1$ be the reduction factor of boundary sides associated with one refinement step, i.e., $r$ is the smallest number such that $\operatorname {diam}(\gamma ') \leq r\,\operatorname {diam}(\gamma )$ for all the subsides $\gamma '$ of a refined boundary side $\gamma$. Given $0 < \theta _g < 1$, let $\hat{{\mathcal{N}}_H}$ satisfy

$$\begin{equation} \operatorname {osc}_\partial (g;{\hat{{\mathcal{N}}_H}}) \ge \theta _g \operatorname {osc}_\partial (g;{{\mathcal{N}}_H}). \tag{6.4} \end{equation}$$

Let ${\mathcal{T}}_h$ be a refinement of ${\mathcal{T}}_H$ such that all the elements having a node in $\hat{{\mathcal{N}}_H}$ are refined. Then, for $\beta _g := (1 - (1-r^2)\theta _g^2)^{1/2} < 1$, we have

$$\begin{equation} \operatorname {osc}_\partial (g;{{\mathcal{N}}_h}) \le \beta _g \operatorname {osc}_\partial (g;{{\mathcal{N}}_H}). \tag{6.5} \end{equation}$$

Remark 6.6.

Using two bisection steps for the refinement of a boundary element ensures a reduction factor $r \leq 1/2$.

7. Weighted Poincaré Inequality

In this section we present a constructive proof of Proposition 2.4, which yields an explicit dependence of constants on the geometry of stars $\star$ and makes the paper self-contained. For interior stars, (2.7) can be obtained as a consequence of the weighted Poincaré inequality of 10, Theorem 8.8, which actually goes back to 13. Our result can also be derived from the element-oriented Lemmas 5.1-3 of 5 via a contradiction argument, which does not reveal the role of stars’ geometry though.

We split the proof of Proposition 2.4 into two parts. Lemma 7.1 establishes the result for reference boundary stars, and Proposition 2.4 then follows by scaling with a piecewise linear transformation. Lemma 7.2, instead, establishes the assertion for interior stars, and Proposition 2.4 follows by scaling with a linear mapping (a dilation) which preserves the constraint $\int v \phi _i = 0$.

Lemma 7.1.

Let $n \in \mathbb{N}$ be fixed and let $\varphi (x)$, $x \in [-1,1]$, be the polygonal function that satisfies $\varphi (\sqrt 2 \cos (\frac{\pi }{4}+\frac{k\, \pi }{2n})) = \sqrt {2}\sin (\frac{\pi }{4}+\frac{k\, \pi }{2n})$, $k = 0,1,\dots ,n$. Let ${\hat{\omega }}$ be the reference star with $n$ elements defined by

$$\begin{equation*} {\hat{\omega }}:= \left\{(x,y): |x| <y<\varphi (x), -1<x<1 \right\} \end{equation*}$$

(see Figure 6), and let $\phi$ be the continuous piecewise linear function on ${\hat{\omega }}$ that is equal to $1$ at the origin and zero on the graph of $\varphi$ for $-1 \le x \le 1$. Then for any smooth function $v$ with $v=0$ on $\{y=|x|\}$ we have

$$\begin{equation*} \int _{\hat{\omega }}|v|^2 \le C \int _{\hat{\omega }}|\nabla v|^2 \phi ^\alpha , \end{equation*}$$

for any $0 \le \alpha < 2$, where $C$ depends only on $\alpha$.

Proof.

Let $(x,y) \in {\hat{\omega }}$, and extend $v$ and $\phi$ by zero for $y<|x|$ and $y>\varphi (x)$, respectively, to write

$$\begin{equation*} v(x,y)^2 = \left(\int _0^y \partial _y v(x,t) \,dt\right)^2 \le \int _0^{\sqrt 2} |\nabla v(x,t) |^2 \phi (x,t)^\alpha \, dt \, \int _{|x|}^y \frac{1}{\phi (x,t)^\alpha }\, dt, \end{equation*}$$

which implies

$$\begin{equation*} \int _{\hat{\omega }}v(x,y)^2 \le \int _{-1}^{1} \left\{ \int _0^{\sqrt 2} |\nabla v(x,t) |^2 \phi (x,t)^\alpha \, dt \int _{|x|}^{\varphi (x)} \int _{|x|}^y \frac{1}{\phi (x,t)^\alpha }\, dt \,dy \right\} dx. \end{equation*}$$

Let us assume momentarily that for all $(x,t)\in {\hat{\omega }}$

$$\begin{equation} \phi (x,t) \ge \frac{1}{2} \left(1 - \frac{t}{\varphi (x)}\right) = \frac{1}{2} \frac{\varphi (x)-t}{\varphi (x)}. {\tag{7.1}} \end{equation}$$

Then, the assertion of the lemma follows from

$$\begin{align*} \int _{|x|}^{\varphi (x)} \int _{|x|}^y \frac{1}{\phi (x,t)^\alpha }\, dt \,dy &\le 2^\alpha \int _{|x|}^{\varphi (x)} \int _{|x|}^y \left(\frac{\varphi (x)}{\varphi (x)-t}\right)^\alpha \, dt \,dy \\ &= 2^\alpha \frac{\varphi (x)^\alpha (\varphi (x)-|x|)^{2-\alpha }}{2-\alpha } \le \frac{2^\alpha }{2-\alpha }, \end{align*}$$

which holds for $0 \le \alpha < 2$.

It only remains to prove (7.1). Observe first that for each fixed $-1<x<1$, the function $\phi (x,t)$ is a concave function of $t$ on $|x| \le t \le \varphi (x)$ (the domain of $\phi (x,\cdot )$). For $t$ at the endpoints of this interval we have

$$\begin{equation*} \text{(i)} \quad \phi (x,|x|) = 1-|x|, \qquad \qquad \text{(ii)} \quad \phi (x,\varphi (x)) = 0. \end{equation*}$$

Let us now define $h_x(t) = \frac{(1-|x|)(\varphi (x)-t)}{\varphi (x)-|x|}$, i.e., $h_x(t)$ is the linear function (in $t$) that coincides with $\phi (x,t)$ at $t=|x|$ and $t=\varphi (x)$, the endpoints of the interval of interest (see Figure 7).

Then, by (i)–(ii) and the concavity of $\phi (x,t)$, we have that

$$\begin{equation*} \phi (x,t) \ge h_x(t),\qquad \forall t \in [|x|,\varphi (x)]. \end{equation*}$$

Now, since $|x| < \varphi (x) < 2-|x|$ on $(-1,1)$ (see Figure 7 (left)), we have

$$\begin{equation*} 0 < \varphi (x) - |x| \le 2-2|x| = 2 (1-|x|), \qquad \forall x \in (-1,1). \end{equation*}$$

Consequently,

$$\begin{equation*} \phi (x,t) \ge h_x(t) \ge \frac{1}{2} \frac{(1-|x|)(\varphi (x)-t)}{1-|x|} = \frac{\varphi (x)}{2}\left(1 - \frac{t}{\varphi (x)}\right) \ge \frac{1}{2}\left(1 - \frac{t}{\varphi (x)}\right), \end{equation*}$$

where in the last inequality we used $\varphi (x) \ge 1$. This completes the proof.

■

To establish the result for interior stars we first observe that any interior star $\star$ is star-shaped with respect to a ball $B_i$ of radius comparable with $h_i$. That is, there exists a ball $B_i \subset \star$ such that

(1)

for all $x\in \star$ and $y \in B_i$, the segment joining $x$ and $y$ is contained in $\star$;

(2)

the radius $r_i$ of $B_i$ satisfies $r_i \ge C h_i$ for some constant $C$ solely depending on the minimum angle of the mesh.

To prove these properties, we consider an interior star $\star$, as depicted in Figure 8, and a triangle $T$ of $\star$. We realize that if $r<b=h\sin \alpha$, then the ball centered at $x_i$ with radius $r$ is contained in the halfspace generated by the boundary side of $T$ and $x_i$. Now let $\rho _i$ and $\alpha _i$ denote the minimum element side $h$ and angle $\alpha$ in $\star$, which in turn satisfy $\rho _i\ge Ch_i$ and $\alpha _i\ge C$ due to mesh shape regularity. It then suffices to take $r_i=\frac{1}{2} \rho _i\sin \alpha _i$.

By dilation, the star $\star$ can be mapped into a star $\hat{\omega }_i$ of diameter 1 which is star-shaped with respect to a ball $\hat{B}_i$ whose radius is bounded below by a positive constant independent of $i$ and above by 1. Lemma 7.2 states the assertion of Proposition 2.4 for a generic star ${\hat{\omega }}$ of diameter 1. Since the constraint $\int v \phi _i = 0$ is preserved by a linear map, a scaling argument completes the proof.

Lemma 7.2.

Let ${\hat{\omega }}$ be a star of diameter $1$ which is star-shaped with respect to a ball $\hat{B}$. Let $\phi$ be the piecewise linear function that is equal to one at the center node and zero outside ${\hat{\omega }}$. Then for any smooth function $v$ with $\int v \phi =0$ we have

$$\begin{equation*} \int _{\hat{\omega }}|v|^2 \le C \int _{\hat{\omega }}|\nabla v|^2 \phi ^\alpha \end{equation*}$$

for any $0 \le \alpha < 2$, where $C$ depends only on $\alpha$ and the radius $\hat{r}$ of $\hat{B}$.

Proof.

Let $B$ be the ball with the same center as $\hat{B}$ but with radius $\hat{r}/2$. Let $\rho$ be a nonnegative function in $C_0^\infty (B)$ with $\int \rho =1$. Then, defining $\bar{v} = \int v \rho$ and by a simple change of variables we have that

$$\begin{align*} v(y)-\bar{v} &= \int _{\hat{\omega }}(v(y)-v(z)) \rho (z) \, dz = \int _{\hat{\omega }}\int _0^1 (y-z) \cdot \nabla v(y+s(z-y))\rho (z)\,ds\,dz \\ &= \int _{\hat{\omega }}\int _0^1 \frac{1}{s^{3}} (y-x) \cdot \nabla v(x) \rho \left(y + \frac{x-y}{s}\right) \,dx\,ds = \int _{\hat{\omega }}G(x,y) \nabla v(x) \, dx, \end{align*}$$

where $G(x,y) = \int _0^1\frac{1}{s^{3}} (y-x)\rho (y + \frac{x-y}{s})\,ds$.

On the one hand, if $y \in {\hat{\omega }}$ and if $\frac{|x-y|}{s}> 2$, then the integrand in the definition of $G(x,y)$ is 0 because in that case $\rho (y + \frac{x-y}{s})=0$. Then,

$$\begin{equation*} |G(x,y)| \le \left\|\rho \right\|_{L^\infty } \int _{\frac{|x-y|}{2}}^1 \frac{|x-y|}{s^{3}} \le C |x-y|^{-1}. \end{equation*}$$

On the other hand, if $x$ is outside the convex hull of $B \cup \{y\}$ we have $G(x,y)$ $=0$. Instead, if $x$ is inside that hull, since the radius of $B$ is $\hat{r}/2$, it follows that $|x-y| \le C \phi (x)$ (with $C$ depending only on $\hat{r}$). Therefore,

$$\begin{align*} |v(y) - \bar{v}| &\le \int _{\hat{\omega }}|G(x,y)| | \nabla v(x)| \, dx = \int _{\hat{\omega }}\frac{|G(x,y)|}{\phi (x)^{\alpha /2}} \, \phi (x)^{\alpha /2} | \nabla v(x)|\, dx \\ &\le C \int _{\hat{\omega }}|x-y|^{-1-\alpha /2} \, \phi (x)^{\alpha /2} | \nabla v(x)|\, dx = K * \left(\phi ^{\alpha /2}| \nabla v|\right) (y) \end{align*}$$

where $K(x) = C |x|^{-1-\alpha /2} \chi _{\{|x|\le 3\}}$ is a function in $L^1(\mathbb{R}^2)$ for any $\alpha <2$. Consequently, by Young’s inequality

$$\begin{equation*} \left\|v - \bar{v}\right\|_{L^2} \le \left\|K\right\|_{L^1} \left\| \phi ^{\alpha /2}| \nabla v|\right\|_{L^2}. \end{equation*}$$

Let us assume first that $\int _\star v = 0$. Then $\left\|v\right\|_{L^2}\le \left\|v - \bar{v}\right\|_{L^2}$ and the result is proved in this case. If $\int _\star v \phi = 0$ instead, by defining $\tilde{v} = \int _{\hat{\omega }}v/|{\hat{\omega }}|$ we have

$$\begin{align*} \tilde{v} &= \frac{\int _{\hat{\omega }}v}{|{\hat{\omega }}|} = \frac{\int _{\hat{\omega }}v}{|{\hat{\omega }}|} - \frac{\int _{\hat{\omega }}v\phi }{|{\hat{\omega }}|} = \frac{1}{|{\hat{\omega }}|} \int _{\hat{\omega }}v (1 - \phi ) \\ &\le \frac{1}{|{\hat{\omega }}|} \left\|v\right\|_{L^2} \left( \int _{\hat{\omega }}(1-\phi )^2\right)^{1/2} =: \frac{\nu }{|{\hat{\omega }}|^{1/2}} \left\|v\right\|_{L^2} , \end{align*}$$

where $\nu < 1$. Then $\left\| \tilde{v}\right\|_{L^2} \le \nu \left\| v\right\|_{L^2}$, which by the triangle inequality implies

$$\begin{equation*} \left\|v\right\|_{L^2} \le \left\|v - \tilde{v}\right\|_{L^2} + \left\| \tilde{v}\right\|_{L^2} \le C \left\| \nabla v\right\|_{\phi ^\alpha } + \nu \left\|v\right\|_{L^2}, \end{equation*}$$

and the claim follows.

■

Remark 7.3.

The assertion of Lemma 7.2, as well as its proof, is valid for any domain ${\hat{\omega }}\subset \mathbb{R}^d$ ($d\ge 2$) which is star-shaped with respect to a ball $B$, and $\phi (x) = \text{dist}(x,\partial {\hat{\omega }})$. The proof of Lemma 7.2 is inspired by an argument of Durán and Muschietti 8 for finding an explicit right inverse of the divergence operator.

8. Numerical experiments

The purpose of this section is to illustrate with several examples the good performance of the estimators introduced in §2, as well as the quasi-optimality of the meshes generated by Algorithm C.

8.1. Implementation

For all the computations presented in this section, the flexible adaptive finite element toolbox ALBERT 1516 was used. The discrete system (2.4) is assembled by the standard assembling tools of ALBERT, and the resulting linear system is solved by a conjugate gradient method using hierarchical basis preconditioning.

In order to compute the solutions $\eta _i$ to the discrete local problems, we loop around each center node, collect data about this star and assemble the small linear system (2.8) which is solved by Gaussian elimination. Since for interior nodes a basis of ${\mathcal{P}}_0^2({\omega _i})$ is not directly at hand, we assemble the system in ${\mathcal{P}}^2({\omega _i})$ and impose the constraint $\int _{\omega _i} \eta _i \phi _i = 0$; recall that ${\mathcal{P}}^2({\omega _i})$ stands for the space of piecewise quadratic polynomials vanishing on $\partial {\omega _i}$. All integrals involving only discrete functions are computed exactly, whereas those also involving data functions are computed element-wise by a quadrature formula which is exact for polynomials of degree 7.

In Algorithm C we have to mark nodes in such a way that (4.1) holds. This can be easily achieved by a slight modification of the marking algorithms proposed in 611: Let $0<\kappa <1$ be a given parameter, where $\kappa \approx 0.1$: We use the same marking procedure for marking nodes due to data oscillation. Here we use the set $\hat{{\mathcal{N}}_H}$ of already marked nodes as an initial guess for $\hat{{\mathcal{N}}_H}$ and set $\mathtt{sum} \coloneq \operatorname {osc}(f;{\hat{{\mathcal{N}}_H}})^2$, or $\mathtt{sum} \coloneq \operatorname {osc}_\partial (g;{\hat{{\mathcal{N}}_H}})^2$. In our experiments we use the parameters (see Convergent Algorithm C):

$$\begin{equation*} \theta = \theta _f = \theta _g = 0.3, \qquad \nu _f = \nu _g = 0.7. \end{equation*}$$

In the examples that follow, the new estimator $\mathcal{ E}$, without any scaling constants, is compared with the Bank-Weiser estimator $\mathcal{ E}_{BW}$ 4 and the residual type estimator $\mathcal{ E}_R$ 18. For getting enhanced results from the Bank-Weiser estimator, the jumps of the fluxes have to be equilibrated 12417. This equilibration procedure additionally needs solutions of small algebraic equations on stars, and thus extra computing time. For the comparison we only use the standard Bank-Weiser estimator computed element-wise with quadratic bubble functions for all interior edges. The residual-type error estimator is scaled by a factor of 0.25 in all problems.

The true error is computed by using the following identity for the solution $u$ of (1.2) and any $\varphi \in H^1(\Omega )$:

$$\begin{align*} \left|\negthinspace \left|\negthinspace \left|{u - \varphi } \right|\negthinspace \right|\negthinspace \right|^2 &= \int _\Omega f (u - 2\,\varphi ) + \int _{\partial \Omega } A \nabla u \cdot \nu \, (u - 2\,\varphi ) + \left|\negthinspace \left|\negthinspace \left|{\varphi } \right|\negthinspace \right|\negthinspace \right|. \end{align*}$$

When dealing with singular solutions, this formula has the following advantages for computing an approximation to the true error $\left|\negthinspace \left|\negthinspace \left|{u - u_h} \right|\negthinspace \right|\negthinspace \right|$ by quadrature: In the interior the gradient of $u$ is replaced by $u$, which leads to a better approximation. On the boundary, the solution $u$ is smooth (even zero) where $\nabla u$ is singular.

8.2. Example: Interior layer

The objective in this example is to illustrate the behaviour of Algorithm C when studying a problem with a rough right hand side and strongly varying boundary values. For this, we consider the exact solution

$$\begin{equation*} u(x) = \operatorname {arctan}(60\,(|x|^2 -1)) \end{equation*}$$

for $A = I$, $\Omega$ as in Figure 9 and boundary data $g$ and right hand side $f$ defined accordingly. This solution exhibits a strong interior layer and the right hand side $f$, as well as $g$, is rather rough. Due to a bad resolution of data on coarse grids, the discrete solution exhibits strong oscillations (top pictures in Figure 9) which vanish after improving data resolution (bottom pictures in Figure 9).

Table 1 presents the degrees of freedom (DOFs), error and the effectivity indices for the new estimator ($\mathcal{ E}^k$), the Bank-Weiser estimator ($\mathcal{ E}_{BW}^k$) and the residual type estimator ($\mathcal{ E}_R^k$) for the $k$-th adaptive cycle. All estimators suffer from approximating rough data $f$ and $g$ on coarse grids, since they ignore data oscillation. It is remarkable that after resolving data, the effectivity index of the new estimator $\mathcal{ E}^k$ is close to $1$; recall that the Definition 2.5 does not involve any constant. As can be seen in Figure 9, the grids obtained after some iterations of Algorithm C are strongly graded and do not exhibit any uniformity typical of superconvergence effects.

Although data oscillation is rather big, there were no additional nodes marked for reducing the data oscillation error. This means that in all adaptive cycles $k$ where nodes had to be marked due to oscillation of $f$, $\operatorname {osc}(f;{\hat{{\mathcal{N}}_k}})$ already satisfies (4.3) where $\hat{{\mathcal{N}}_k}$ denotes the set of nodes marked due to (4.1). The same holds for oscillation of boundary data. The average reduction rates of $\left|\negthinspace \left|\negthinspace \left|{u-u_{k}} \right|\negthinspace \right|\negthinspace \right|$, $\operatorname {osc}(f;{{\mathcal{N}}_k})$, and $\operatorname {osc}_\partial (g;{{\mathcal{N}}_k})$ are equal to 0.68, 0.57, and 0.53, respectively. Although data is rather rough, data oscillation is decreasing faster than the error (see Figure 10). This behaviour sheds light on why adaptive methods seem to converge in practice even without taking data oscillation into account.

Quite revealing is Figure 11. It shows the asymptotic relation $\left|\negthinspace \left|\negthinspace \left|{u-u_{k}} \right|\negthinspace \right|\negthinspace \right| = C$ DOFs$(k)^{-1/2}$ typical of quasi-optimal meshes in 2d and thus illustrates the quasi-optimal numerical complexity of Algorithm C. In the $\log$-$\log$ plot the optimal decay of $\left|\negthinspace \left|\negthinspace \left|{u-u_{k}} \right|\negthinspace \right|\negthinspace \right|$ is a straight line with slope $-1/2$, which is also plotted in Figure 11 for reference.

Finally, in Table 2 we give a comparison of CPU times on an SGI O2 for the four finest grids. The times shown for $\mathcal{ E}^k$ include, unlike those for $\mathcal{ E}_{BW}^k$ and $\mathcal{ E}_{R}^k$, the computation of data oscillation. The computation of $\mathcal{ E}^k$ requires less than half the time needed for solving the discrete system. However, it is about four and five times more costly than the Bank-Weiser and residual type estimator, respectively. This accounts for the fact that each triangle belongs to three stars and for the additional cost of computing data oscillation. Since the effectivity index of $\mathcal{ E}^k$ is close to $1$, this additional effort pays off in practice.

8.3. Example: Crack problem

For analyzing a problem with a singularity of the type $r^{1/2}$, we consider the domain $\Omega =\{|x|+|y|<1\}\backslash \{0\le x\le 1, y=0\}$ with a crack and the exact solution $u$ in polar coordinates $(r,\theta )$:

$$\begin{equation*} u(r,\theta )=r^{{\textstyle \frac{1}{2}}} \sin {\textstyle \frac{\theta }{2}} - {\textstyle \frac{1}{2}} r^2\sin ^2\theta . \end{equation*}$$

We solve (1.2) with $A=I$ and $f=1$, and nonvanishing boundary values $g$ on $\partial \Omega$. Table 3 demonstrates the strong reliability of the new estimator. The effectivity index of the new estimator is $\approx$ 1 for all grids, whereas those of the Bank-Weiser estimator and the residual type estimator are completely different from the previous example. The good performance of the new estimator even on coarse grids accounts for the fact that data oscillation is small in this example. As before, no elements are marked due to data oscillation. The average reduction rates for $\left|\negthinspace \left|\negthinspace \left|{u-u_{k}} \right|\negthinspace \right|\negthinspace \right|$, $\operatorname {osc}(f;{{\mathcal{N}}_k})$, and $\operatorname {osc}_\partial (g;{{\mathcal{N}}_k})$ are 0.72, 0.61, and 0.41, respectively (see Figure 12). Finally, the quasi-optimality of the meshes produced by Algorithm C is shown in Figure 13, revealing the asymptotic performance of $\left|\negthinspace \left|\negthinspace \left|{u-u_{k}} \right|\negthinspace \right|\negthinspace \right| = C\,$DOFs$(k)^{-1/2}$.

8.4. Example: Discontinuous coefficients

We invoke the formulas derived by Kellogg 9 to construct an exact solution of an elliptic problem with piecewise constant coefficients and vanishing right hand side $f$. We use these formulas in the particular case $\Omega =(-1,1)^2$, $A = a_1 I$ in the first and third quadrants, and $A = a_2 I$ in the second and fourth quadrants with $R= a_1/a_2 \cong 161.4476387975881$. For a representation of the exact solution $u$ and boundary values $g$, see 11. The solution behaves like $r^{0.1}$ at the origin and thus is barely in $H^1(\Omega )$.

Due to the large ratio of $a_1$ and $a_2$ the estimator $\mathcal{ E}$ is underestimating the error by 30% (compare Table 4). Altogether, this is not that surprising since the equivalence constants are sensitive to the ratio of smallest and largest eigenvalues of $A$ within a star. We stress that $A$ is not monotone around the origin but rather presents a checkerboard pattern that leads to the worst possible singularity.

Algorithm C produces a convergent sequence of discrete solutions with an average reduction rate of 0.91 and 0.84 for $\left|\negthinspace \left|\negthinspace \left|{u-u_{k}} \right|\negthinspace \right|\negthinspace \right|$ and $\operatorname {osc}_\partial (g;{{\mathcal{N}}_k})$. Note that $\operatorname {osc}(f;{{\mathcal{N}}_k}) = 0$. Usually, only some nodes at the origin are marked for refinement, resulting in a highly graded grid. Starting with a uniform mesh size of 1 on the macro triangulation, the minimal mesh size after 20 iterations is $\approx$ $6.75e-07$ with 1669 DOFs in total. Figure 14 demonstrates that the grids and associated numerical complexity are quasi-optimal: $\left|\negthinspace \left|\negthinspace \left|{u-u_{k}} \right|\negthinspace \right|\negthinspace \right|= C\,$DOFs$(k)^{-1/2}$ is valid asymptotically (the performance of an optimal method is again indicated by the additional straight line).

Figure 1.

Star ${\omega _i}$ and sides forming ${\Gamma _i}$ (indicated by thick lines) for an interior node $x_i$ and for a boundary one.

Figure 2.

Refinement of a triangle by the newest-vertex bisection leading to (4.2). Dashed lines indicate the refinement edges (opposite the newest vertex).

Figure 3.

Refinement by bisecting all triangles twice.

Figure 4.

Refinement of a star by bisection ensuring (4.2) (left); ${\omega _S^h}= \operatorname {supp}\varphi _S$ (middle); ${\omega _i^h}= \operatorname {supp}\varphi _i$ (right).

Figure 5.

Refinement of a star by two bisections whenever the center node is the newest vertex of all the triangles (left); ${\omega _S^h}= \operatorname {supp}\varphi _S$ (middle); ${\omega _i^h}= \operatorname {supp}\varphi _i$ (right).

Figure 6.

Reference stars for boundary nodes with $n$ elements for $n=1,2,3,4$. The nodes lie on the circle of radius $\sqrt 2$.

Figure 7.

A generic boundary reference star, and the functions $\phi (x,t)$, $\frac{1}{2}\left(1-t/\varphi (x)\right)$, and $h_x(t)$ for a fixed $x$.

Figure 8.

All interior stars $\star$ are star-shaped with respect to a ball centered at $x_i$ and of radius $r_i\ge Ch_i$.

Figure 9.

Graphs of the discrete solutions and meshes for adaptive cycles 1 (top left), 4 (top right), 7 (bottom left), and 10 (bottom right) of the interior layer example. The diameter of the domain is $\approx$ 2.5 and the height of the layer in the bottom right picture is $\approx$ 3. The scaling (0.17) of the height of the graphs is the same in all pictures.

Table 1.

Error and effectivity indices for the interior layer example.

$k$	DOFs	$\left|\negthinspace \left|\negthinspace \left|{u-u_{k}} \right|\negthinspace \right|\negthinspace \right|$	${\mathcal{ E}^k}/{\left|\negthinspace \left|\negthinspace \left|{u-u_{k}} \right|\negthinspace \right|\negthinspace \right|}$	${\mathcal{ E}_{BW}^k}/{\left|\negthinspace \left|\negthinspace \left|{u-u_{k}} \right|\negthinspace \right|\negthinspace \right|}$	${\mathcal{ E}_R^k}/{\left|\negthinspace \left|\negthinspace \left|{u-u_{k}} \right|\negthinspace \right|\negthinspace \right|}$
0	18	$2.294e+02$	0.62	0.66	1.45
3	269	$8.478e+01$	0.52	0.86	1.25
6	1225	$3.809e+01$	0.46	0.86	0.80
9	6377	$7.946e+00$	0.85	1.76	1.12
12	43190	$2.256e+00$	1.02	1.91	1.27
13	93799	$1.495e+00$	1.03	1.91	1.28
14	206422	$9.965e-01$	1.04	1.93	1.28
15	452991	$6.707e-01$	1.04	1.91	1.28

Table 2.

CPU times on an SGI O2 for assembling and solving the discrete system (A&S), and for computing the estimators $\mathcal{ E}^k$, $\mathcal{ E}_{BW}^k$, $\mathcal{ E}_{R}^k$. The times for $\mathcal{ E}^k$ include, unlike those for $\mathcal{ E}_{BW}^k$ and $\mathcal{ E}_{R}^k$, the computation of data oscillation.

$k$	DOFs	time A&S	time $\mathcal{ E}^k$	time $\mathcal{ E}_{BW}^k$	time $\mathcal{ E}_R^k$
12	43190	22.90	11.15	2.99	2.46
13	93799	57.27	24.78	6.43	5.33
14	206422	143.97	55.96	14.30	11.67
15	452991	343.57	125.41	31.76	25.86

Figure 10.

Error and data oscillation for the interior layer example.

Figure 11.

Quasi-optimality of Algorithm C. The optimal decay is indicated by the dashed line with slope $-$1/2.

Table 3.

Error and effectivity indices for the Crack problem.

$k$	DOFs	$\left|\negthinspace \left|\negthinspace \left|{u-u_{k}} \right|\negthinspace \right|\negthinspace \right|$	${\mathcal{ E}^k}/{\left|\negthinspace \left|\negthinspace \left|{u-u_{k}} \right|\negthinspace \right|\negthinspace \right|}$	${\mathcal{ E}_{BW}^k}/{\left|\negthinspace \left|\negthinspace \left|{u-u_{k}} \right|\negthinspace \right|\negthinspace \right|}$	${\mathcal{ E}_R^k}/{\left|\negthinspace \left|\negthinspace \left|{u-u_{k}} \right|\negthinspace \right|\negthinspace \right|}$
0	6	$7.839e-01$	0.97	0.45	0.80
2	40	$4.079e-01$	0.99	0.58	0.77
4	124	$2.418e-01$	1.00	0.76	0.89
6	300	$1.341e-01$	1.01	0.93	0.97
8	961	$7.403e-02$	1.04	0.81	0.97
10	7845	$3.024e-02$	1.01	0.76	0.90
12	31123	$1.561e-02$	1.02	0.75	0.91

Figure 12.

Error and data oscillation for the Crack problem.

Figure 13.

Quasi-optimality of Algorithm C for the Crack problem. The optimal decay is indicated by the dashed line with slope $-$1/2.

Table 4.

Error and effectivity indices for the discontinuous coefficient example.

$k$	DOFs	$\left|\negthinspace \left|\negthinspace \left|{u-u_{k}} \right|\negthinspace \right|\negthinspace \right|$	${\mathcal{ E}^k}/{\left|\negthinspace \left|\negthinspace \left|{u-u_{k}} \right|\negthinspace \right|\negthinspace \right|}$	${\mathcal{ E}_{BW}^k}/{\left|\negthinspace \left|\negthinspace \left|{u-u_{k}} \right|\negthinspace \right|\negthinspace \right|}$	${\mathcal{ E}_R^k}/{\left|\negthinspace \left|\negthinspace \left|{u-u_{k}} \right|\negthinspace \right|\negthinspace \right|}$
0	13	$1.051e+00$	0.81	1.11	6.55
5	133	$5.210e-01$	0.68	0.95	5.63
10	317	$3.306e-01$	0.68	0.90	5.89
15	677	$2.299e-01$	0.71	0.89	6.26
20	1669	$1.531e-01$	0.67	0.83	5.55

Figure 14.

Quasi-optimality of Algorithm C for the discontinuous coefficient example. The optimal decay is indicated by the dashed line with slope $-$1/2.