Higher-order adaptive virtual element methods with contraction properties

1.   Introduction

Adaptive Finite Element Methods (AFEM) for self-adjoint coercive problems written in the form

iterate the sequence

to produce better and better approximations of $u$. Their practical efficiency is corroborated by sound theoretical results of convergence, complexity, and optimality, which in various cases (such as, e.g., conforming $h$-versions) completely explain the behaviour of the adaptive algorithms ^{[11,13,14,15,18]}.

The standard AFEM realization preserves the conformity of the initial mesh, at the expense of performing a completion step in ${\sf REFINE}$: In addition to elements marked for refinement due to their contribution to the global error estimator, other elements are refined. Without this step, one would obtain nonconforming meshes, containing elements with hanging nodes.

In the new perspective opened by the introduction of Virtual Element Methods (VEM) ^[3,4], elements with hanging nodes can be viewed as polygons with aligned edges, carrying virtual (i.e., non-accessible) functions together with standard polynomial functions. The potential advantage is that all activated degrees of freedom are motivated by error reduction, not just by geometric reasons. On the other hand, in this transformation of an adaptive FEM into an adaptive VEM, one loses the availability of a general convergence theory, which so far is lacking (although results on a posteriori error estimates ^[8,12] have been obtained, together with efficient practical recipes for refining polytopal meshes ^[2,9,10]).

Such a shift in perspective inspired the recent papers ^[5,6], devoted to the analysis of an adaptive VEM generated by the successive newest-vertex bisections of triangular elements without applying completion, in the lowest-order case (polynomial degree $k = 1$). Despite the simple geometric setup, the investigation faced some VEM-specific obstacles in the analysis, giving answers that could prove useful in the study of more general adaptive VEM discretizations. For instance, a VEM solution $u_\mathcal{T} \in \mathbb{V}_\mathcal{T} \subset \mathbb{V}$, defined by the Galerkin projection

satisfies an a posteriori error bound of the type

where $\eta_\mathcal{T}(u_\mathcal{T})$ is a residual-type error estimator, $S_ \mathcal{T}(u_\mathcal{T}, u_\mathcal{T})$ is the stabilization term that makes the discrete bilinear form $\mathcal{B}_ \mathcal{T}(u_\mathcal{T}, v_\mathcal{T})$ coercive in $\mathbb{V}$, and for simplicity we assume piecewise constant data on the mesh $\mathcal{T}$. Unfortunately, the term $S_ \mathcal{T}(u_\mathcal{T}, u_\mathcal{T})$ need not reduce under a mesh refinement, as $\eta_\mathcal{T}^2(u_\mathcal{T})$ does: This makes the convergence analysis problematic. However, one of the key results obtained in ^[5] states that $S_ \mathcal{T}(u_\mathcal{T}, u_\mathcal{T})$ is dominated by $\eta_\mathcal{T}^2(u_\mathcal{T})$, i.e.,

provided an assumption of admissibility of the non-conforming meshes generated by successive refinements is fulfilled; such a restriction, which appears to have little practical impact, amounts to requiring the uniform boundedness of the global index of all hanging node, a useful concept introduced in ^[5] to hierarchically organize the set of hanging nodes. Once the a posteriori error bound is reduced to

the convergence analysis becomes feasible, and a contraction property is proven to hold for a linear combination of the (squared) energy norm of the error and the (squared) residual estimator.

The purpose of this paper is to extend the results in ^[5] to the case of VEMs of order $k\ge2$ built on triangular meshes. The problem at hand is again a variable-coefficient, second-order self-adjoint elliptic equation with Dirichlet boundary conditions. The geometric concept of hanging node (a vertex for some elements, contained inside an edge of some other elements) is replaced by a functional one, referring to the degrees of freedom associated with the node; once the meaning of hanging node is clarified, the definition of global index of a node, and its role in the analysis, is similar to the one given in ^[5].

A significant difference with respect to the content of that paper concerns the control of the stabilization term, which does not involve only the residual estimator, but a new term, called the virtual inconsistency estimator and denoted by $\Psi_ \mathcal{T}(u_\mathcal{T})$. It measures the projection error, upon local spaces of polynomials, of certain expressions depending on the operator coefficients and the discrete solution; it vanishes when $k = 1$ or when the coefficients are constant. The new stabilization bound, which we derive under an admissibility assumption of the mesh, takes the form

which leads to the a posteriori, stabilization-free error control

Correspondingly, we obtain the convergence of the adaptive VEM of order $k$ by proving a contraction result for a linear combination of (squared) energy norm of the error, (squared) residual estimator, and (squared) virtual inconsistency estimator.

Similarly to ^[5], we assume here that the data $\mathcal{D}$ of our boundary-value problem are piecewise polynomials of degrees related to $k-1$, on the initial mesh $\mathcal{T}_0$ and consequently on each mesh $\mathcal{T}$ derived by newest-vertex bisection. This is not a restriction, since we propose to insert the adaptive VEM procedure just described, which we now consider as a module ${\sf GALERKIN}$, into an outer loop ${\sf AVEM}$ of the form

$[\mathcal{T}, u_\mathcal{T}] = {\tt AVEM}(\mathcal{T}_0, \epsilon_0, \omega, \text{tol})$

$j = 0$

while $\epsilon_j > \frac{1}{2} \mathrm{tol}$ do

  $[\hat{\mathcal{T}}_j, \hat{\mathcal{D}}_j] = {\sf DATA}(\mathcal{T}_j, \mathcal{D}, \omega \epsilon_j)$

  $[{\mathcal{T}}_{j+1}, \mathcal{D}_{j+1}] = {\sf GALERKIN}(\hat{\mathcal{T}}_j, \hat{\mathcal{D}}_j, \epsilon_j)$

  $\epsilon_{j+1} \gets\frac{1}{2}\epsilon_j$

  $j \gets j +1$

end while

return

where the module ${\sf DATA}$ produces, via greedy-type iterations, a piecewise polynomial approximation of the input data with prescribed accuracy, defined on a suitable refinement of the input partition. Manifestly, the target accuracy is matched after a finite number of calls to ${\sf DATA}$ and ${\sf GALERKIN}$. Properties of complexity and quasi-optimality of this two-loop algorithm are investigated in ^[6] in the linear case $k = 1$. We plan to do the same for the case $k \geq 2$ in a forthcoming paper.

The outline of this paper is as follows. In Sections 2 and 3, we introduce the model boundary-value problem, and its discretization by an enhanced version of the VEM (^[1]). In Section 4 we define the global index of a node, and we formulate the admissibility assumption on the mesh. Two essential properties for bounding the stabilization term are established in Section 5. The a posteriori error estimators are defined in Section 6, whereas stabilization-free a posteriori error estimates are proven in Section 7. In Section 8, we investigate how the a posteriori error estimators are reduced under mesh refinement. These properties are needed to justify the refinement strategy in our adaptive module ${\sf GALERKIN}$, which is described in Section 9. In Section 9, we discuss the proof of convergence of the loop ${\sf GALERKIN}$. The paper ends with some numerical experiments, reported in Section 11.

2.   VEM spaces of order $k\ge2$

We consider the following Dirichlet boundary value problem in a polygonal domain $\Omega$,

where $A \in (L^\infty(\Omega))^{2\times 2}$ is symmetric and uniformly positive definite in $\Omega$, $c \in L^{\infty}(\Omega)$ and non-negative in $\Omega$, $f \in L^2 (\Omega)$. Data will be denoted by $\mathcal{D} = (A, c, f)$. The variational formulation of this problem is written as

where $(\cdot, \cdot)$ is the scalar product in $L^2(\Omega)$ and $\mathcal{B}(u, v) : = a(u, v) + m(u, v)$ is the bilinear form associated with Problem (2.1), i.e.,

We denote the energy norm as ${\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert {\cdot} \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert} = \sqrt{\mathcal{B}(\cdot, \cdot)}$, which satisfies

for suitable $0 < c_\mathcal{B} \le c^\mathcal{B}$.

Remark 2.1. For the sake of simplicity, in (2.1) we consider the Poisson problem with vanishing Dirichlet conditions on the whole boundary domain. The extension to generic Dirichlet and/or Neumann/Robin boundary conditions does not pose conceptual difficulties. In the numerical examples, we actually provide experiments with more general Dirichlet boundary conditions.

In order to find a discrete approximation of the solution of Problem (2.2), we firstly introduce a fixed initial partition $\mathcal{T}_0$ on the domain $\overline{\Omega}$ made of triangular elements $E$. We will denote by $\mathcal{T}$ any refinement of $\mathcal{T}_0$ obtained by a finite number of newest-vertex element bisections. We underline that we are not requiring $\mathcal{T}$ to be a conforming mesh, since hanging nodes may arise in the refinement. The classification of nodes, which will play a crucial role in the proofs presented in this paper, is postponed in Section 4.

According to the Virtual Element theory ^[3], an element $E$ of the triangulation can be viewed as a polygon with more than three edges, if some hanging nodes are sitting on its boundary. We can then denote by $\mathcal{E}_E$ the set of edges $e$ of element $E$ and $\mathcal{E}: = \bigcup_{E\in\mathcal{T}}\mathcal{E}_E$. We finally define the diameter of an element $E$ as $h_E = |E|^{1/2}$ and $h = \max_{E\in \mathcal{T}}\{h_E\}$.

We introduce the functional spaces needed to apply the VEM. We start by defining the space of functions on the boundary of $E$, $\mathbb{V}_{\partial E, k}$, which is constituted by the functions that are continuous on the boundary of $E$ and that, when restricted to any edge of $\partial E$, are polynomials of degree $k > 0$, i.e.,

Then, we define the "enhanced" VEM space in $E$, as done in ^[1], such that

where $\mathbb{P}_k(E)\setminus \mathbb{P}_{k-2}(E)$ is the space spanned by the monomials of degree equal to $k$ and $k-1$, and $\Pi^\nabla_{E} \; :\; H^1(E)\rightarrow \mathbb{P}_k(E)$ is the projector defined by

We remark that $\mathbb{V}_{E, k}$ contains the polynomial space of degree $k$ on $E$ and its dimension is

where $n^E_e$ is the number of edges of $E$. We notice that in the case $k > 1$ a function $v$ in $\mathbb{V}_{E, k}$ is uniquely defined by

● the set of the values at the vertices of $E$;

● the set of the values at the $k-1$ equally-spaced internal points on each edge of $\partial E$;

● the set of the moments $\frac{1}{|E|}\int_{E} v(\mathit{\boldsymbol{x}}) m(\mathit{\boldsymbol{x}}) d\mathit{\boldsymbol{x}} \; \forall m \in \mathcal{M}_{k-2}(E)$,

where the set $\mathcal{M}_{p}(E)$, $p \geq 0$, is defined as

We will denote by $\mathit{\boldsymbol{\mu}}_p(E, v) = \left(\frac{1}{|E|}\int_{E} v(\mathit{\boldsymbol{x}}) m(\mathit{\boldsymbol{x}}) d\mathit{\boldsymbol{x}} : m \in \mathcal{M}_{p}(E) \setminus \mathcal{M}_{p-1}(E) \right)$ the vector of the moments of $v$ of order $p$. By $|\mathit{\boldsymbol{\mu}}_p(E, v) |$ we will denote the $l^2$-norm of this vector.

We can now introduce the global discrete space as

On $\mathcal{T}$ we need also to give the definition of the space of piecewise polynomial functions on $\mathcal{T}$

and its subspace

which plays a crucial role in the forthcoming analysis.

We now introduce a series of projectors that will be used in the rest of the paper. For any $E\in \mathcal{T}$, we denote by $\Pi^0_{p, E}: L^2(E)\rightarrow \mathbb{P}_p(E)$ the $L^2(E)$-orthogonal projector onto the space of polynomial of degree $p$ on $E$. Thanks to the choice of the enhanced space $\mathbb{V}_{E, k}$ (2.4), we remark that $\Pi^0_{k, E} v$ and $\Pi^0_{k-1, E}\nabla v$ can be computed for any function $v \in \mathbb{V}_{E, k}$, see ^[1] for the details. To simplify the notation, in the following we will drop the symbol $E$ from $\Pi^0_{k, E}$ when no confusion arises. The global $L^2$-orthogonal projector is denoted by $\Pi^0_{p, \mathcal{T}}: L^2(\Omega)\rightarrow \mathbb{W}^p_\mathcal{T}$.

We can also define the Lagrange interpolation operator $\mathcal{I}_E:\mathbb{V}_{E, k}\rightarrow \mathbb{P}_k(E)$ on $E$, which builds a polynomial of degree $k$ using the $3 k$ degrees of freedom on the boundary of $E$ and the moments of order $\leq k-3$, since

Moreover, we will denote by $\mathcal{I}_\mathcal{T}:\mathbb{V}_\mathcal{T}\rightarrow \mathbb{W}^k_{\mathcal{T}}$ the Lagrange interpolation operator that restricts to $\mathcal{I}_E$ on each $E\in\mathcal{T}$.

3.   Discretization with data of degree $k-1$

In the rest of this paper, we assume that data $\mathcal{D} = (A, c, f)$ are piecewise polynomials of degree $k-1$ on the initial partition $\mathcal{T}_0$, hence on each partition $\mathcal{T}$ obtained by newest-vertex refinement. Their values on each element of the triangulation will be denoted by

We here define the bilinear forms that we need for the Galerkin discretization problem, starting from $a_E, m_E : \mathbb{V}_{E, k} \times \mathbb{V}_{E, k} \rightarrow \mathbb{R}$, such that

We also introduce the symmetric bilinear form $s_E: \mathbb{V}_E \times \mathbb{V}_E \rightarrow \mathbb{R}$ as

where $\{\mathit{\boldsymbol{x}}_i\}^{\overline{\mathcal{N}}_E}_{i = 1}$ indicates the set of the degrees of freedom on the boundary of $E$. Indeed, we remark that in this case the stabilization term can be built without using the internal degrees of freedom, as shown in ^[7]. We assume for $s_E$ the existence of two positive constant $c_s$ and $C_s$ independent on $E$, such that

We define the local stabilizing form as

and the global stabilization form

From (3.1), we get

where $\left| {{\, \cdot \, }} \right|$ denotes the broken $H^1$-seminorm over $\mathcal{T}$. Thus, we can now define the bilinear form $\mathcal{B}_{\mathcal{T}}(\cdot, \cdot)$, $\mathcal{B}_{\mathcal{T}}: \mathbb{V}_{\mathcal{T}} \times \mathbb{V}_{\mathcal{T}}\rightarrow \mathbb{R}$, as

with $\gamma$ independent of $\mathcal{T}$ satisfying $\gamma \ge \gamma_0$ for some fixed $\gamma_0 > 0$. For the loading term we introduce $\mathcal{F}_\mathcal{T}: \mathbb{V}_\mathcal{T} \rightarrow \mathbb{R}$ as

since $f_E$ has been already approximated with a polynomial of degree $k-1$. Note that the equality in (3.3) remains true if $f_E$ is an approximation of $f$ of degree $k$ on $E$.

We have now defined all the forms that appear in the discrete formulation of the Problem (2.2). It reads as

The bilinear form $\mathcal{B}_\mathcal{T}$ is continuous and coercive, hence, there exists a unique and stable solution of the Problem (3.4). Furthermore, the following result extends Lemma 2.6 in ^[5].

Lemma 3.1 (Gakerkin quasi-orthogonality). For any $v\in \mathbb{V}_\mathcal{T}$ and $w \in \mathbb{V}^0_\mathcal{T}$, it holds

Consequently,

where $u$ is the solution of (2.2) and $u_\mathcal{T}$ the solution of (3.4).

4.   The index of a node

A crucial concept, firstly introduced in ^[5] for the case $k = 1$, is the global index of a node: It will be used in the proofs of Section 5. In order to extend its definition to the case $k > 1$, we preliminarily introduce some useful definitions.

Let

be the reference element and denote by $\hat{R}_{\hat{E}, k}$ the $k$-lattice built on $\hat{E}$, i.e.,

Considering the affine function $F_E: \hat{E}\rightarrow E$ mapping the reference element onto an element $E\in \mathcal{T}$, we define the physical lattice on $E$ by

and the set of proper nodes of $E$ as the points of the physical lattice sitting on the boundary of $E$, i.e.,

Observe that we implicitly assume that $k\ge2$ is sufficiently small so that interpolation on equally spaced nodes is numerically stable.

Next, we denote by $\mathcal{H}_E$ the set of hanging nodes of $E$, i.e., the set of points $\mathit{\boldsymbol{x}} \in \partial E$ that are not proper nodes of $E$, but that are proper nodes of some other element $E'$, i.e.,

Finally, let $\mathcal{N}_E: = \mathcal{P}_E \cup \mathcal{H}_E$ be the set of all nodes sitting on $E$.

At the global level, $\mathcal{N} : = \bigcup_{E \in\mathcal{T}} \mathcal{N}_E$ will be the set of all nodes of the triangulation $\mathcal{T}$, which we split into the set $\mathcal{P}: = \left\{\mathit{\boldsymbol{x}}\in \mathcal{N}: \mathit{\boldsymbol{x}}\in \mathcal{P}_E \ \forall E \text{ containing }\mathit{\boldsymbol{x}}\right\}$ of the proper nodes of$\mathcal{T}$, and the set $\mathcal{H}: = \mathcal{N}\setminus \mathcal{P}$ of the hanging nodes of$\mathcal{T}$.

Next, let us clarify what happens when a hanging node is created. Let $S$ be an element edge that is being refined, i.e., split into two contiguous edges $S^-$ and $S^+$. Before the refinement, $S$ contains $k+1$ equally-spaced nodes $\mathit{\boldsymbol{\xi}}_n$, $n = 1, \dots k+1$: The endpoints and the $k-1$ internal ones. After the refinement, $S$ contains $2k+1$ nodes, precisely $k+1$ equally-spaced nodes on each sub-edge $S^\pm$, with the midpoint in common; see Figure 1. The spacing of the 'old' nodes on $S$ was $\frac{|S|}k$ (where $|S|$ denotes the length of $S$), whereas the spacing of the 'new' nodes is $\frac{|S|}{2k}$. Consequently, $k+1$ of these nodes coincide with those initially on $S$, and the new nodes introduced in the refinement are only $k$. We will denote these latter by $\mathit{\boldsymbol{\zeta}}_i$, $i = 1, \dots, k$.

This suggests the following definition.

Definition 4.1 (closest neighbors of a node). With the previous notation, if $\mathit{\boldsymbol{x}} : = \mathit{\boldsymbol{\zeta}}_i$ is created as the midpoint of the segment $[\mathit{\boldsymbol{x}}', \mathit{\boldsymbol{x}}'']: = [\mathit{\boldsymbol{\xi}}_{n_i}, \mathit{\boldsymbol{\xi}}_{n_i+1}]$ for some $n_i$, we define the set $\mathcal{B}(\mathit{\boldsymbol{x}}): = \{ \mathit{\boldsymbol{x}}', \mathit{\boldsymbol{x}}''\}$.

We are ready to give the announced definition of global index of a node of the triangulation $\mathcal{T}$.

Definition 4.2. (global index of a node). Given a node $\mathit{\boldsymbol{x}}\in\mathcal{N}$, we define its global index $\lambda$ recursively as follows:

● If $\mathit{\boldsymbol{x}}$ is a proper node, then $\lambda(\mathit{\boldsymbol{x}}): = 0$;

● If $\mathit{\boldsymbol{x}}$ is a hanging node, with $\mathit{\boldsymbol{x}}', \mathit{\boldsymbol{x}}'' \in \mathcal{B}(\mathit{\boldsymbol{x}})$, then set

Figure 2 shows the evolution of the global index after three refinements in the cases $k = 2$ (a) and $k = 3$ (b). We remark that, for instance, the midpoint of the horizontal edge is a proper node in case (a), and a hanging node in case (b).

The largest global index in $\mathcal{T}$ will be denoted by $\Lambda_\mathcal{T}: = \max_{\mathit{\boldsymbol{x}}\in \mathcal{N}}\{\lambda(\mathit{\boldsymbol{x}})\}$. In this paper, as in ^[5], we will consider sequences of successively refined triangulations $\{ \mathcal{T} \}$ whose global index does not blow up.

Assumption 4.3. There exists a constant $\Lambda > 0$ such that, for any triangulation $\mathcal{T}$ generated by successive refinements of $\mathcal{T}_0$, it holds

Any such triangulation will be called $\Lambda$-admissible.

5.   Two key properties

In this section we discuss the validity of some results for the degree $k > 1$ that will be used in the rest of the paper. We will highlight in particular the differences from the case $k = 1$.

Proposition 5.1 (scaled Poincaré inequality in $\mathbb{V}_\mathcal{T}$). There exists a constant $C_P > 0$, independent of $\mathcal{T}$, such that

Proof. Let $E\in \mathcal{T}$ be an element of the triangulation. If $E$ is an element of the original partition $\mathcal{T}_0$, all its vertices are proper nodes. Otherwise, $E$ has been generated after some refinements by splitting an element $\widetilde{E}$ into two elements, $E$ and $E'$. Let $L$ be the common edge shared by $E$ and $E'$. If $L$ is not further refined, then all the nodes on $L$ are proper because they are shared by $E$ and $E'$. If $L$ is refined and $k$ is even, then the midpoint of $L$ is a proper node.

So, let us consider the case $k$ odd and let us assume that $L$ is refined $M\ge1$ times. We focus in particular on the internal node $\bar{\mathit{\boldsymbol{x}}}$ of $L$ is at distance $\frac{\left| {{L}} \right|}{k}$ from one of the endpoints, Figure 3 shows the case $k = 3$. This point belongs to one of the $M+1$ intervals in which $L$ is refined, having width $\left| {{L}} \right|/2^s$, for some $1\le s \le M$. We remark that $s$ depends on how $L$ has been refined (in the case of uniform refinements of $L$, one has $2^s = M +1$). We localize the chosen node $\bar{\mathit{\boldsymbol{x}}}$ in $L$ by defining an $m \ge 0$ such that

or, equivalently,

The interval going from $\frac{\left| {{L}} \right|\; m}{2^s}$ to $\frac{\left| {{L}} \right|(m +1)}{2^s}$ is an edge for a smaller element $E'$, thus it contains $k-1$ internal nodes. Since they are equi-spaced, their positions are at

By taking $n = 2^s - m\; k$, which is compatible with conditions (5.2), we conclude that one of the internal nodes of $E'$ coincides with $\bar{\mathit{\boldsymbol{x}}}$.

This guarantees that $E$ has at least one proper node $\mathit{\boldsymbol{x}}$ on its boundary. By hypothesis $v(\mathit{\boldsymbol{x}}) = 0$, and so we can apply the classical Poincaré inequality,

that concludes the proof. □

Remark 5.2. The previous proof exploits the fact that when $k > 1$, each element of the triangulation contains at least a proper node. This differs from the case $k = 1$ in which the edges do not contain internal nodes, and then elements with all hanging nodes as vertices are admissible. As a further difference from the case $k = 1$, we highlight that in Proposition 5.1 the constant $C_P$ does not depend on the constant $\Lambda$, whose existence has been introduced in Assumption 4.3.

The next result we are going to establish is a hierarchical representation of the interpolation error $v - \mathcal{I}_E v$ on the boundary $\partial E$ of an element $E \in \mathcal{T}$. Assume that $v \in {\mathbb{V}_{E, k}}$, and let $L$ be a side of the triangle $E$; for simplicity, in the sequel the restriction of $v$ to $L$, which is a piecewise polynomial of degree $k$, will be still denoted by $v$. The subsequent bisections of $L$ which generate the nodes in $\mathcal{N}_E \cap L$ allow us to write the difference $(v - \mathcal{I}_E v)_{|L}$ telescopically as

here, ${\cal I}_0 = { \mathcal{I}_E}_{|L}$, ${\cal I}_{J_L}$ is the identity operator, whereas ${\cal I}_jv$ for $1 \leq j \leq J_L-1$ is the piecewise polynomial of degree $k$ which interpolates $v$ on the partition of $L$ of level $j$, namely the partition formed by sub-edges of length $\leq \frac{|L|}{2^j}$.

In order to understand the structure of the detail $({\cal I}_j - {\cal I}_{j-1})v$, assume that $S$ is a sub-edge of $L$ of length $= \frac{|L|}{2^{j-1}}$, which is split into two sub-edges $S^\pm$ of length $= \frac{|L|}{2^{j}}$ (see Figure 1). On $S$ we have two interpolation operators, namely

and

which coincides with the interpolation operator $\mathcal{I}^-:C^0(S^-)\rightarrow \mathbb{P}_k(S^-)$ when restricted to $S^-$ and with the analogous operator $\mathcal{I}^+$ when restricted to $S^+$. With the notation introduced just before Definition 4.1, we can quantify the discrepancy between the two interpolation operators by defining the $k$ basis functions

See Figure 4 for a graphical representation of these functions in the cases $k = 1$ (a), $k = 2$ (b), $k = 3$ (c).

Hence, the difference between the two interpolation operators on $S$ can be written as

where $d$ is defined as

The values of $\mathcal{I}v$ at the $k$ nodes $\mathit{\boldsymbol{\zeta}}_i$ are a linear combination of the values of $\mathcal{I}v$ at the $k+1$ nodes $\mathit{\boldsymbol{\zeta}}_n$, where $\mathcal{I}v$ coincides with $v$. Thus, there exist coefficients $\alpha_{i, n}$ such that

The explicit values of these coefficients in the case $k = 2$ for the two new nodes $\mathit{\boldsymbol{\zeta}}_1$ and $\mathit{\boldsymbol{\zeta}}_2$ are given in these expressions:

where $\mathit{\boldsymbol{\xi}}_i \le \mathit{\boldsymbol{\zeta}}_i \le \mathit{\boldsymbol{\xi}}_{i+1}$, $i = 1, 2$. Similarly, in the case $k = 3$, we get

where again $\mathit{\boldsymbol{\xi}}_i \le \mathit{\boldsymbol{\zeta}}_i \le \mathit{\boldsymbol{\xi}}_{i+1}$, $i = 1, 2, 3$. Figure 5 shows both cases. We notice that the coefficients $\alpha_{i, n}$ depend only on the relative positions of the nodes on $S$, not on the level $j$ of refinement.

Summarizing, at the level $j$ of refinement of the edge $L$, we get

where $\mathcal{H}_{L, j}$ is the set of hanging nodes on $L$ created at the level $j$ of refinement, whereas

Summing-up over the levels and recalling (5.3), we obtain

where $\mathcal{H}_{L} = \mathcal{H}_{E} \cap L$, whence

We now introduce the subspace of ${\mathbb{V}_{E, k}}$

which contains $v- \mathcal{I}_E v$ by definition of $\mathcal{I}_E$. On $X_E$, we have two norms, namely the seminorm $|w |_{1, E}$ (which is a norm on $X_E$ due to the vanishing of $w$ at the three vertices of $E$) and the norm

Note that, due to Assumption 4.3, the dimension of $X_E$ is uniformly bounded by a constant depending on $\Lambda$; furthermore, the number of possible patterns of hanging nodes on $\partial E$, which determines the details $d(w, \mathit{\boldsymbol{x}})$, is also bounded in terms of $\Lambda$. As a consequence, the two norms are equivalent, with equivalence constants depending on $\Lambda$. Therefore,

Since $v - \mathcal{I}_E v \in X_E$ and $d (v - \mathcal{I}_E v, \mathit{\boldsymbol{x}}) = d(v, \mathit{\boldsymbol{x}})$ for any $\mathit{\boldsymbol{x}} \in \mathcal{H}_E$, we obtain

Summing-up over all the elements of the triangulation, we arrive at the following result.

Lemma 5.3 (global interpolation error vs hierarchical errors). There exists a constant $C_D > 0$ depending on $\Lambda$ but independent of the triangulation $\mathcal{T}$ such that

Next, we introduce the interpolation operator

where $\mathbb{V}_\mathcal{T}^0$ is defined in (2.8), by the following conditions:

● $(\mathcal{I}_\mathcal{T}^0 v)(\mathit{\boldsymbol{x}}) = v(\mathit{\boldsymbol{x}})$ for all $\mathit{\boldsymbol{x}} \in {\cal P}$,

● $\mathit{\boldsymbol{\mu}}_p(E, \mathcal{I}_\mathcal{T}^0 v) = \mathit{\boldsymbol{\mu}}_p(E, v)$ for all $0 \leq p \leq k-3$ and for all $E \in \mathcal{T}$.

These conditions uniquely identify $\mathcal{I}_\mathcal{T}^0 v$. Indeed, if $\mathit{\boldsymbol{x}} \in {\cal H}$ is generated by a refinement of level $j$ of an edge $L$ (say, $\mathit{\boldsymbol{x}} = \mathit{\boldsymbol{\zeta}}_i$ with the notation introduced before Definition 4.1), then $(\mathcal{I}_\mathcal{T}^0 v)(\mathit{\boldsymbol{x}})$ can be expressed in terms of the values of $\mathcal{I}_\mathcal{T}^0 v$ at the $k+1$ nodes (say, $\mathit{\boldsymbol{\xi}}_n$) created at the previous levels of refinement of $L$, using the same coefficients as in formula (5.5), i.e.,

and so on recursively.

The following result provides a representation of the error $\mathcal{I}_\mathcal{T} v - \mathcal{I}_\mathcal{T}^0 v$.

Lemma 5.4. It holds

where $\delta(v, \mathit{\boldsymbol{x}}): = v(\mathit{\boldsymbol{x}}) - (\mathcal{I}_\mathcal{T}^0 v)(\mathit{\boldsymbol{x}})$.

Proof. Consider an element $E \in \mathcal{T}$. Recall that by construction it holds $\mathit{\boldsymbol{\mu}}_p(E, \mathcal{I}_E v) = \mathit{\boldsymbol{\mu}}_p(E, v) = \mathit{\boldsymbol{\mu}}_p(E, \mathcal{I}_\mathcal{T}^0 v)$, whence $\mathit{\boldsymbol{\mu}}_p(\mathcal{I}_E v - \mathcal{I}_\mathcal{T}^0 v, E) = \mathit{\boldsymbol{0}}$ for all $0 \leq p \leq k-3$. Consequently,

If $\mathit{\boldsymbol{x}}\in\mathcal{P}_E$, $(\mathcal{I}_E v)(\mathit{\boldsymbol{x}}) = v(\mathit{\boldsymbol{x}})$, hence

Summing on all the elements of the partition, we get

since if $\mathit{\boldsymbol{x}}\in\mathcal{P}$, $(\mathcal{I}^0_ \mathcal{T} v)(\mathit{\boldsymbol{x}}) = v(\mathit{\boldsymbol{x}})$. This concludes the proof. □

Concatenating Lemmas 5.3 and 5.4, we can prove the second key property of this section.

Proposition 5.5 (comparison between interpolation operators). Let $\mathcal{T}$be $\Lambda$-admissible. Then, there exists a constant $C_I > 0$, depending on $\Lambda$, but independent of $\mathcal{T}$, such that

Proof. Given a function $v \in \mathbb{V}_ {\mathcal{T}}$, by the triangle inequality

so it is enough to bound the last norm on the right-hand side. To this end, considering the vectors

and recalling the two Lemmas, the proof can be concluded if we show that

Given $\mathit{\boldsymbol{x}} \in {\cal H}$, assume that it is generated by a refinement of level $j$ of an edge $L$ (say, $\mathit{\boldsymbol{x}} = \mathit{\boldsymbol{\zeta}}_i$ with the notation introduced before Definition 4.1). Writing $v^*: = \mathcal{I}^0_\mathcal{T}v$ for short, and exploiting formulas (5.4) and (5.5), we get

Thus, we can build a matrix $\mathit{\boldsymbol{W}}: l^2(\mathcal{H})\rightarrow l^2(\mathcal{H})$ such that $\mathit{\boldsymbol{\delta}} = \mathit{\boldsymbol{W}}\mathit{\boldsymbol{d}}$, and we just need to prove that

We now organize the hanging nodes with respect to the global index $\lambda \in [1, \Lambda_\mathcal{T}]$. Calling $\mathcal{H}_\lambda = \{\mathit{\boldsymbol{x}} \in \mathcal{H}: \lambda(\mathit{\boldsymbol{x}}) = \lambda\}$, and $\mathcal{H} = \bigcup_{1\le \lambda\le \Lambda_\mathcal{T}}\mathcal{H}_\lambda$, the matrix $\mathit{\boldsymbol{W}}$ can be factorized in lower triangular matrices $\mathit{\boldsymbol{W}}_\lambda$, that change the nodes of level $\lambda$, leaving the others unchanged. In particular,

where $\mathit{\boldsymbol{W}}_{1}$ is just the identity matrix $\mathit{\boldsymbol{I}}$, whereas each other matrix $\mathit{\boldsymbol{W}}_\lambda$ differs from the identity only in the rows of block $\lambda$. In each of these rows, all entries are zero, but the entries $\alpha_{i, n}$ in the off-diagonal part and $1$ on the diagonal. In order to estimate $\mathit{\boldsymbol{W}}_\lambda$, we use the Hölder inequality $||\mathit{\boldsymbol{W}}_\lambda||^2_2\le ||\mathit{\boldsymbol{W}}_\lambda||_1||\mathit{\boldsymbol{W}}_\lambda||_\infty$. From the construction of $\mathit{\boldsymbol{W}}_\lambda$ have that

where in the last inequality it has been used the fact that a hanging node of global index $< \lambda$ may appear at most 5 times on the right-hand side of (5.9), since at most five edges meet at a node [5, Proposition 3.2]. These bring us to the following bound

and the proof is concluded. □

6.   A posteriori error estimator

With the aim of discussing the a posteriori error analysis, and following ^[12], we define the a posteriori error estimators, starting from the internal residual over an element $E$, i.e.,

for any $v\in \mathbb{V}_{E, k}$. We highlight that in the case $k = 1$, with piecewise constant data, the diffusion term in the residual vanishes. Furthermore, we define the jump residual over $e$, where $e$ is an edge shared by two elements $E_1$ and $E_2$ of the partition $\mathcal{T}$, as

where $\mathit{\boldsymbol{n_i}}$ denotes the unit normal vector to $e$ pointing outward with respect to $E_i$; we set $j_\mathcal{T}(e; v) = 0$ of $e \in \partial \Omega$. Then, let the local residual estimator associated with $E$ be

and the global residual estimator as the sum of the local residuals

In contrast to what has been done for the case $k = 1$, we also need to introduce the virtual inconsistency terms, defined by

as well as their sum

7.   A posteriori error estimates

In this section we present one of the main results of this paper, a stabilization-free a posteriori error bound. In this view, we firstly start by introducing the classical Clément operator upon the space $\mathbb{V}^0_ \mathcal{T}$, $\tilde{\mathcal{I}}^0_\mathcal{T}: \mathbb{V}\rightarrow \mathbb{V}^0_ \mathcal{T}$; it is defined at the proper nodes on the skeleton of $\mathcal{T}$ as the average of the target function on the support of the associated basis functions, whereas the internal moments (if any) coincide with those of the target function.

The scaled Poincaré inequality (Proposition 5.1) and Proposition 5.5 guarantee the validity of the error estimate for $\tilde{\mathcal{I}}^0_\mathcal{T}$. Given these propositions, its proof does not involve the polynomial degree $k$, hence, it does not change with respect to the one presented in ^[5].

Lemma 7.1 (Clément interpolation estimate). $\forall v \in \mathbb{V}$, it holds

where the hidden constant depends on $\Lambda$ but not on $\mathcal{T}$.

We can now prove the following results, which is similar to Theorem 13 in ^[12], but with a slightly modified proof.

Proposition 7.2 (upper bound). There exists a constant $C_\mathit{\text{apost}} > 0$, independent of $u$, $\mathcal{T}$, $u_\mathcal{T}$ and $\gamma$, such that

Proof. For any $v \in\mathbb{V}$, using the definition of Problem (2.2), we have that

where $v_\mathcal{T}: = \tilde{\mathcal{I}}^0_\mathcal{T} v \in \mathbb{V}_\mathcal{T}^0$. The first term can be written as

The addend $I_1$ can be expressed as

which can be bounded by using Lemma 7.1,

On the other hand, noting that

and applying again Lemma 7.1 and the stability of the Clément operator in the $H^1$ norm, the addend $I_2$ can be bounded as follows:

Looking now at the term $II$, we have by Lemma 3.1

Finally, by taking $v: = u- u_\mathcal{T} \in \mathbb{V}$, we get

which, using the coercivity of $\mathcal{B}$, concludes the proof. □

We now report a bound for the local residual estimator, proved in [12, Theorem 16].

Proposition 7.3 (local lower bound). There holds

where $w_E: = \{E':|\partial E \cap \partial E'| \neq 0\}$. The hidden constant is independent of $\gamma$, $h$, $u$ and $u_\mathcal{T}$.

Summing on all the elements of the partition, we get the following corollary.

Corollary 7.4 (global lower bound). There exists a constant $c_\mathit{\text{apost}} > 0$, independent of $u$, $\mathcal{T}$, $u_\mathcal{T}$ and $\gamma$, such that

In the following proposition we present a bound of the stabilization term. We remark that in the case $k = 1$ the inconsistency term does not appear.

Proposition 7.5 (bound of the stabilization term). There exists a constant $C_B > 0$ independent of $\mathcal{T}$, $u_\mathcal{T}$ and $\gamma$, such that

Proof. From the Definition (3.2) of the form $\mathcal{B}_\mathcal{T}$ and from (3.4), $\forall w \in \mathbb{V}^0_\mathcal{T}$ it holds

Defining $e_\mathcal{T}: = u_\mathcal{T}-w$, we get

We notice that

and

By substituting (7.5) and (7.6) into (7.4), it results

With the same strategy used in ^[5], for any $\delta > 0$, we get

where

Posing now $w = \mathcal{I}^0_\mathcal{T}u_\mathcal{T}$ and applying Proposition 5.1, we get

whereas Proposition 5.5 yields

so we obtain

for a suitable constant $C_B > 0$. □

Combining Propositions 7.2 and 7.5, we arrive at the following key result.

Corollary 7.6 (stabilization-free a posteriori error upper bound). It holds

where $C_{U_1} = C_\mathit{\text{apost}} \left(\frac{C_B}{\gamma^2} +1\right)$ and $C_{U_2} = C_\mathit{\text{apost}}\frac{C_B}{\gamma^2}$.

Remark 7.7. Note that the chosen stabilization affects the value of the constant $c_{apost}$, which in principle may depend on the polynomial degree and the geometry of the mesh. However, this dependence is under control; indeed, (i) we are not proposing a $p$-method, so the polynomial degree $k$ is fixed, (ii) the refinement procedure is obtained by newest-vertex bisection, which guarantees shape regularity on the refined elements, (iii) Assumption 4.3 enforces an upper bound on the number of hanging nodes on each edge.

8.   The effect of a mesh refinement

In view of the convergence analysis of the adaptive algorithm ${\sf GALERKIN}$, in this section we analyse the effect of refining the partition $\mathcal{T}$ by applying one or more newest-vertex bisections to some of its elements. Specifically, in Sect. 8.1 we prove that the residual estimator (6.2) is reduced by a fixed fraction (up to an addend proportional to the stabilization term) when the element $E$ is split into two elements by one bisection. We prove a similar result for the inconsistency term estimator (6.4), provided a suitable number of bisections is applied to $E$. Next, in Sect. 8.2 we establish a quasi-orthogonality property in the energy norm between the solutions on two partitions, one being a refinement of the other.

8.1.   Reduction of estimators under refinement

Let us consider an element $E$ in $\mathcal{T}$ which is bisected into elements $E_1$ and $E_2$; the refined partition containing these two elements will be denoted by ${\mathcal{T}_*}$. Given $v \in \mathbb{V}_\mathcal{T}$, we notice that $v$ is known on $\partial E$, and in particular at the new vertex of $E_1$ and $E_2$ produced by the bisection. Denoting by $e = E_1\cap E_2$ the new edge, we associate a function $v_*\in\mathbb{V}_{\mathcal{T}_*}$ to $v$ such that $v_*|_{\partial E} = v|_{\partial E}$, $v_*|_{e}\in \mathbb{P}_1(e)$, and $\mathit{\boldsymbol{\mu}}_p(E_i, v_*) = \mathit{\boldsymbol{\mu}}_p(E, v)$ for all $0 \leq p \leq k-2$ and for $i = 1, 2$. In the following we will write $v$ instead of $v_*$ when no confusion arises.

8.1.1.   The residual estimator

Let $\eta_\mathcal{T}(E; v, \mathcal{D})$ be defined in (6.2) and $\eta_{\mathcal{T}_*}(E; v, \mathcal{D})$ be the sum of the local residual estimators on the two newly formed elements, defined as follows:

where we recall that $h_ {E_i} = \frac{1}{\sqrt{2}}h_E$, $i = 1, 2$ We notice that, since $\mathcal{D}$ does not change under refinement, the functions $f_{E_i} = f_E|_{E_i}$, $c_{E_i} = c_E|_{E_i}$ and $A_{E_i} = A_E|_{E_i}$ will be denoted again by $f_E$, $c_E$ and $A_E$, respectively.

Lemma 8.1 (local residual estimator reduction). There exist constants $\mu_r \in (0, 1)$ and $c_{er, 1} > 0$ such that for any $v \in \mathbb{V}_\mathcal{T}$

where $S_{\mathcal{T}(E)}(v, v): = \sum_{E' \in \mathcal{T}(E)}S_{E'}(v, v)$ with $\mathcal{T}(E): = \{E'\in \mathcal{T}: \mathcal{E}_E \cap \mathcal{E}_{E'} \not = \emptyset\}$.

Proof. Recalling the Definition (6.1), we have the following residuals

Writing

we get, for any $\epsilon > 0$,

The second term can be bounded by using the inverse inequality and the minimality of $\Pi_{k-1, E_i}^0$ as follows:

while, for the last term, using the Poincaré inequality we have

Finally, taking an appropriate value of $\epsilon$ and setting $\mu: = \frac{1+ \epsilon}{2}\in (0, 1)$ (for instance, if $\epsilon = \frac{1}{2}$, $\mu = \frac{3}{4}$) we get

where $C > 0$ is a constant.

For the jump condition, we will essentially use the proof given in [5, Lemma 5.2]. In particular, we write $j_{\mathcal{T}_*}(e; v) = j_\mathcal{T}(e; v)+\left(j_{\mathcal{T}_*}(e; v)- j_\mathcal{T}(e, v)\right)$ and for any $\epsilon > 0$

with $T_1: = \sum^2_{i = 1}\sum_{e\in\mathcal{E}_{E_i}}h_{E_i}\left\| {{j_\mathcal{T}(e; v)}} \right\|^2$ and $T_2 : = \sum^2_{i = 1}\sum_{e\in\mathcal{E}_{E_i}}h_{E_i}\left\| {{j_{\mathcal{T}_*}(e; v)-j_\mathcal{T}(e; v)}} \right\|^2$. On the new edge we notice that $j_\mathcal{T}(e; v) = 0$, then,

We now define $\mathcal{T}_*(E_i): = \{E'\in \mathcal{T}_*: \mathcal{E}_{E_i}\cap \mathcal{E}_{E'}\neq \emptyset\}$; for any edge $e\in \mathcal{E}_{E_i}$, we denote by $E_{i, e}\in \mathcal{T}_*(E_i)$ the element such that $e = \partial E_i \cap \partial E_{i, e}$. Then,

where $\hat{E}_{i, e}$ indicates the parent of $E_{i, e}$. Using the trace inequality we have

Using now the minimality property of $\Pi^0_{k-1, E'}$ and $\Pi^0_{k-1, \hat{E}'}$, we easily get as above

which, for a sufficiently small $\epsilon$, concludes the proof. □

From this Lemma and the Lipschitz continuity of the residual estimator with respect to the argument $v$ (whose proof is independent of the used polynomial degree, so we refer to [5, Lemma 5.3]), we immediately deduce the following result.

Proposition 8.2 (residual estimator reduction on refined elements). There exist constants ${ \mu}_r \in (0, 1)$, ${c_{er, 1}} > 0$ and ${ c_{er, 2}} > 0$ independent of $\mathcal{T}$ such that for any $v \in \mathbb{V}_\mathcal{T}$ and $w \in \mathbb{V}_{\mathcal{T}_*}$, and any element $E \in \mathcal{T}$ which is split into two children $E_1, E_2 \in {\mathcal{T}_*}$, one has

8.1.2.   The virtual inconsistency estimator

Given $v \in \mathbb{V}_\mathcal{T}$ and $E \in \mathcal{T}$, consider the two virtual inconsistency terms $\Psi_{ \mathcal{T}, A}(E, v, \mathcal{D})$ and $\Psi_{ \mathcal{T}, c}(E, v, \mathcal{D})$ introduced in (6.3). When $E$ is bisected into $E_1$ and $E_2$, the term $\Psi_{ \mathcal{T}, c}(E, v, \mathcal{D})$ is reduced by a factor $\mu_c < 1$ up to an addend proportional to the stabilization term, i.e., there exists $c_{vi, c} > 0$ such that

This stems from the presence of the factor $h_E$ in front of the norm $||\left (I-\Pi^0_k)(c_E \Pi^0_k v \right)||_{0, E}$, with an argument similar to the one used in the proof of Lemma 8.1.

Due to the lack of the factor $h_E$, a reduction result similar to (8.2) does not hold for $\Psi_{ \mathcal{T}, A}(E, v, \mathcal{D})$. Indeed, since $A_E \Pi_{k-1, E}^0 \nabla v \in \mathbb{P}_{2k-2}(E)$, one may ask whether a constant $\mu < 1$ esists such that

Unfortunately, the answer is no, as it can be seen numerically, working on the reference element $\hat{E}$ by affinity and identifying $\mu^2$ as the largest eigenvalue of a generalized eigenvalue problem. However, the same numerics indicates that if $\hat{E}$ is split into $2^m$ triangles of equal area by $m$ successive levels of uniform bisections, then $\mu^2$ becomes $< 1$ for $m$ large enough, as seen in Table 1.

This is indeed predicted by the following result.

Lemma 8.3. Let $E \in \mathcal{T}$. For any polynomial degree $k \geq 1$ there exists a minimal $m \in \mathbb{N}$ and a constant $\mu = \mu_m < 1$ independent of $E$ such that, if $E$ is partitioned into $2^m$ elements $E_i$ of equal area by $m$ levels of uniform newest vertex bisection, it holds

Proof. Since by construction $h_{E_i} = 2^{-m/2} h_E$, classical approximation results give

for some constant $C_k$ depending on $k$. Replacing $q$ by $q- \Pi_{k-1, E}^0 q$ leaves the left-hand side unchanged, whereas on the right-hand side an inverse inequality yields

One concludes taking as $m$ the smallest integer such that $\mu_m^2 : = C_k C_{\text{inv}, k}2^{-m} < 1$. □

Based on these results, let ${\mathcal{T}_*}^{\! m}$ be a refinement of $\mathcal{T}$ in which the element $E$ has undergone $m$ levels of uniform refinements by newest vertex bisection, and has been replaced by $2^m$ subelements $E_i$. Given $v \in \mathbb{V}_\mathcal{T}$, let us set

Lemma 8.4. There exist constants $\rho_A < 1$ and $c_{vi, A} > 0$ such that for any $v \in \mathbb{V}_{E, k}$

Proof. Write

sum over $i$, and conclude using (8.4) and the usual arguments based on the minimality of the $L^2$-orthogonal projections. □

Let us set

with

Applying a bound similar to (8.2) to the successive level of refinements, we arrive at the following result.

Lemma 8.5. There exist constants $\mu_{vi} < 1$ and $c_{vi, 1} > 0$ such that for any $v \in \mathbb{V}_{E, k}$

Combining this estimate with the Lipschitz continuity property of the virtual inconsistency estimator, we obtain the following result.

Proposition 8.6 (virtual inconsistency estimator reduction on refined elements). There exist constants ${ \mu}_{vi} \in (0, 1)$, ${c_{vi, 1}} > 0$ and ${ c_{vi, 2}} > 0$ independent of $\mathcal{T}$ such that for any $v \in \mathbb{V}_\mathcal{T}$ and $w \in \mathbb{V}_{ {\mathcal{T}_*}^{\!\! m}}$, and any element $E \in \mathcal{T}$ which is split into $2^m$ children $E_i \in \mathbb{V}_{ {\mathcal{T}_*}^{\!\! m}}$, one has

8.2.   Quasi-orthogonality property

Let $u_{\mathcal{T}_*} \in \mathbb{V}_{\mathcal{T}_*}$ be the solution of Problem (3.4) on the refined mesh $\mathcal{T}_*$. Hereafter we establish relations between the two energy errors ${\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert {u - u_\mathcal{T}} \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert}$ and ${\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert {u - u_{\mathcal{T}_*}} \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert}$. The first result follows from Proposition 5.5 and Lemma 3.1; the proof is independent of the used polynomial degree, so we refer to [5, Proposition 5.7].

Proposition 8.7 (comparison of the energy error under refinement). For any $\delta \in (0, 1]$ there exists a constant $C_E > 0$ independent of $\mathcal{T}$ and $\delta$ such that

Next result extends Corollary 5.8 in ^[5].

Proposition 8.8 (quasi-orthogonality of energy errors without stabilization). Given any $\delta \in \left(0, \frac{1}{4}\right]$, there exists $\gamma_{\delta} > 0$ such that for any $\gamma > \gamma_\delta$, it holds

Proof. Let $e: = {\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert {u -u_\mathcal{T}} \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert}$, $e_*: = {\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert {u -u_{\mathcal{T}_*}} \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert}$, $S: = S_\mathcal{T}(u_\mathcal{T}, u_\mathcal{T})$, $S_*: = S_{\mathcal{T}_*}(u_{\mathcal{T}_*}, u_{\mathcal{T}_*})$, $\eta: = \eta_\mathcal{T}(u_\mathcal{T}, \mathcal{D})$, $\Psi: = \Psi_\mathcal{T}(u_\mathcal{T}, \mathcal{D})$, $\Psi_*: = \Psi_{\mathcal{T}_*}(u_{\mathcal{T}_*}, \mathcal{D})$ and $E: = {\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert {u_\mathcal{T} - u_{\mathcal{T}_*}} \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert}$. From Corollary 7.4 and (2.3), we get $\eta^2 \le \frac{S}{c_{apost}} + \frac{e^2}{c_{apost}\; c_\mathcal{B}},$ while, from Proposition 7.5, $S\le \frac{C_B}{\gamma^2}\left(\eta^2 + \Psi^2\right).$ Combining them, we have

Doing the same on $\mathcal{T}_*$ and defining

provided $\gamma^2\ge1$, we get

Employing Proposition 8.7, we obtain

If we define $D: = C_E\left(1 + \frac{1}{\delta} \right)\, \overline{C}$,

By choosing $\gamma$ such that

we get

which concludes the proof by observing that $\frac{1 + 2 \delta}{1-\delta}\le 1 + 4 \delta$ and $\frac{\delta}{1-\delta}\le 2\delta$, when $\delta \le \frac{1}{4}$. □

9.   The module ${\sf GALERKIN}$

Let us consider a $\Lambda$-admissible input mesh $\mathcal{T}_0$, a set of approximated data $\mathcal{D}$ which consist of piecewise polynomials of degree $k-1$ on $\mathcal{T}_0$, and a tolerance $\epsilon > 0$. The call

produces a $\Lambda$-admissible refined mesh $\mathcal{T}$ and the Galerkin approximation $u_\mathcal{T}\in \mathbb{V}_\mathcal{T}$, such as

where $u$ is the solution of Problem (2.2) and $C_G = \sqrt{c^\mathcal{B} \max\left\{C_{U_1}, C_{U_2}\right\}}$, with $c^\mathcal{B}$ is defined in (2.3) and $C_{U_1}, C_{U_2}$ in Corollary 7.6. We obtain it by iterating the sequence

At each step, a $\Lambda-$admissible mesh $\mathcal{T}_j$ and the associated solution $u_j$ of the discrete Problem (3.4) are produced. The process stops when the condition $\eta^2_{\mathcal{T}_j}(u_j, \mathcal{D}) + \Psi^2_{ \mathcal{T}_j}(u_j, \mathcal{D}) \le \epsilon^2$ is reached.

In particular, the modules are defined as follows:

● $[u_\mathcal{T}] = {\sf SOLVE}(\mathcal{T}, \mathcal{D})$ produces the solution of Problem (3.4) with data $\mathcal{D}$;

● $[\{\eta_\mathcal{T}(\cdot; u_\mathcal{T}, \mathcal{D})\}, \{\Psi_\mathcal{T}(\cdot; u_\mathcal{T}, \mathcal{D})\}] = {\sf ESTIMATE}(\mathcal{T}, u_\mathcal{T})$ computes the local estimators on $\mathcal{T}$;

● $[\mathcal{M}] = {\sf MARK}(\mathcal{T}, \{\eta_{\mathcal{T}}(\cdot; u_\mathcal{T}, \mathcal{D}) \}, \{\Psi_\mathcal{T}(\cdot; u_\mathcal{T}, \mathcal{D})\}], \theta)$ implements the Dörfler criterion ^[15] and finds an almost minimal set $\mathcal{M}$ of elements in $\mathcal{T}$ such that

for a given parameter $\theta \in (0, 1)$;

● $[\mathcal{T}_*] = {\sf REFINE}(\mathcal{T}, \mathcal{M}, \Lambda)$ returns a $\Lambda$-admissible refined mesh obtained from $\mathcal{T}$ by suitable newest-vertex bisections of the elements in $\mathcal{M}$, and possibly of other elements to fullfil the $\Lambda$-admissibility condition.

It is worth adding some details about the procedure ${\sf REFINE}$. Let $E \in \mathcal{M}$ be an element marked for refinement. For simplicity, hereafter let us set $\eta: = \eta_\mathcal{T}(E; u_\mathcal{T}, \mathcal{D})$ and $\Psi: = \Psi_\mathcal{T}(E; u_\mathcal{T}, \mathcal{D})$. The refinement of $E$ is performed as follows:

● If $\eta \geq \Psi$, then $E$ is bisected once;

● If $\eta < \Psi$, then $E$ is bisected $m$-times, where $m$ has been introduced in Section 8.1.2 (see Lemma 8.3).

Denote by ${\cal P}(E)$ the partition of $E$ so obtained, and set $\eta_*^2 : = \sum_{E' \in {\cal P}(E)} \eta_{{\cal P}(E)}^2(E'; u_ \mathcal{T}, \mathcal{D})$ and $\Psi_*^2 : = \sum_{E' \in {\cal P}(E)} \Psi_{{\cal P}(E)}^2(E'; u_ \mathcal{T}, \mathcal{D}).$ Then, recalling Lemmas 8.1 and 8.5, one gets when $\eta \geq\Psi$

Indeed, $\Psi$ can be written as $\Psi = \lambda \eta$ for a certain $\lambda\in[0, 1]$ and

In the case $\eta < \Psi$,

In all cases, it holds

which shows that in each marked element the sum of the two estimators is reduced under refinement, up to the stabilization term. Note that for values $k = 2$ or $3$ of the polynomial degree of practical use, two bisections ($m = 2$) are enough when $\eta < \Psi$.

This refinement may create non-admissible hanging nodes, i.e., hanging nodes with global index larger than $\Lambda$. To remove them and guarantee $\Lambda$-admissibility of ${\mathcal{T}_*}$, further refinements should be applied. For the realization of this technical part, we refer to Section 11.1 in ^[6].

The following section proves the convergence of the ${\sf GALERKIN}$ algorithm.

10.   Convergence property of ${\sf GALERKIN}$

Proposition 10.1 (global estimators reduction). Let $u_\mathcal{T} \in \mathbb{V}_\mathcal{T}$ be the solution of the discrete variational Problem (3.4). There exist constants $\rho \in (0, 1)$ and $C_{ger, 1}, C_{ger, 2} > 0$ independent of $\mathcal{T}$ such that, if $\mathcal{T}_*$ is the refinement of $\mathcal{T}$ obtained by applying ${\sf REFINE}$, one has for any $w \in \mathbb{V}_{\mathcal{T}_*}$

Proof. One can reach the conclusion e.g., as in [5, proof of Proposition 5.5], using the bound (9.2) in each element $E$ marked for refinement. □

Theorem 10.2 (contraction property of ${\sf GALERKIN}$). Let $\mathcal{M}\subset \mathcal{T}$ be the set of the marked elements relative to the solution $u_\mathcal{T} \in \mathbb{V}_\mathcal{T}$ of the discrete variational Problem (3.4). If $\mathcal{T}_*$ is the refinement of $\mathcal{T}$ obtained by applying ${\sf REFINE}$, then for $\gamma$ sufficiently large there exist $\alpha \in (0, 1)$ and $\beta > 0$, $\zeta > 0$ such that

Proof. To simplify notation, we set again $e = {\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert {u -u_\mathcal{T}} \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert}$, $e_* = {\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert {u -u_{\mathcal{T}_*}} \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert}$, $S = S_\mathcal{T}(u_\mathcal{T}, u_\mathcal{T})$, $S_* = S_{\mathcal{T}_*}(u_{\mathcal{T}_*}, u_{\mathcal{T}_*})$, $\eta = \eta_\mathcal{T}(u_\mathcal{T}, \mathcal{D})$, $\eta = \eta_{\mathcal{T}_*}(u_{\mathcal{T}_*}, \mathcal{D})$, $\Psi = \Psi_\mathcal{T}(u_\mathcal{T}, \mathcal{D})$, $\Psi_* = \Psi_{\mathcal{T}_*}(u_{\mathcal{T}_*}, \mathcal{D})$ and $E = {\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert {u_\mathcal{T} - u_{\mathcal{T}_*}} \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert}$. From Proposition 8.8,

whereas using Proposition 10.1 and Proposition 7.5, we get

Combining them, we get

which suggests choosing $\beta$ such that

Next, we write

and we invoke Corollary 7.6 to write

which gives

We now choose $\gamma$ and $\delta$ such that

which holds true if

(recall that we already assumed $\gamma^2\geq 1$). Similarly, we choose $\gamma$ and $\delta$ such that

which holds true if $\gamma$ satisfies the first condition in (10.3), whereas $\delta$ satisfies

This proves the result, if we define $\zeta : = \beta - 2\delta$, with $\beta$ defined by (10.2) and $\delta < \frac{\beta}2$, and

The conditions on $\gamma$ and $\delta$ which lead to the desired estimate are given in (8.6), (10.3) and (10.4). □

11.   Numerical experiments

The aim of this numerical test is to confirm the convergence of our ${\sf GALERKIN}$ algorithm. We consider a classical test with an $L-$shaped domain $\Omega = (-1, 1)^2 \setminus (-1, 0)^2$ and the reaction-diffusion problem (2.1), with polynomial coefficients of order one for the case $k = 2$, i.e.,

and polynomials of order two for the case $k = 3$,

The forcing term $f$ and the Dirichlet boundary conditions are chosen so that the solution of the problem results

where $r$ and $\beta$ are the polar coordinates. It is possible to prove that there exists a $p$ with $p \in \left(\frac{2}{k+1}, \frac{2}{k+ 2/3}\right)$ such that $u_\text{ex}\in W_p^k(\Omega)$ when $p\ge 1$, and $u_\text{ex}\in L_p^k(\Omega)$ when $p\in (0, 1)$, where $W_p^k(\Omega)$ and $L_p^k(\Omega)$ indicate respectively Sobolev and Lipschitz spaces. Then, according to the theory of approximation classes ^[16,17], we expect the maximal rate of convergence, i.e., ${\tt Ndofs}^\land{-k/2}$, where ${\tt Ndofs}$ is the number of the degrees of freedom. We apply the adaptive algorithm as described in Section 9 and for the marking strategy (9.1) we consider $\theta = 0.5$. In order to compute the VEM error, we consider the computable quantity:

In Figure 6, we represent the evolution of ${\tt H ^\land 1-error}$ and the estimator terms $\eta_\mathcal{T}(u_ \mathcal{T}, \mathcal{D})$ and $\Psi_ \mathcal{T}(u_ \mathcal{T}, \mathcal{D})$, which confirms the results of Corollary 7.6. Furthermore, we notice that after a transient phase, the error and the estimator terms decays reach asymptotically the theoretical optimal rate ${\tt Ndofs ^\land{-1.0}}$ (for the case $k = 2$) and ${\tt Ndofs ^\land{-1.5}}$ (for the case $k = 3$). In Figure 7, we depict the meshes after 20, 35 and 50 loops of the adaptive algorithm in the case $k = 2$. We highlight the presence of hanging nodes in the different meshes loops.

12.   Conclusions

In this paper, we presented an adaptive VEM of order $k\ge2$ on nonconforming triangular meshes. In the analysis, the space $\mathbb{V}^0_ \mathcal{T}$ of continuous, piecewise polynomials functions of degree $k$ on the triangulation $\mathcal{T}$ plays a fundamental role. Indeed, it is contained in the global VEM space, $\mathbb{V}^0_ \mathcal{T} \subseteq \mathbb{V}_ \mathcal{T}$, and guarantees a quasi-orthogonality property for any refinement ${\mathcal{T}_*}$ of $\mathcal{T}$, since $\mathbb{V}^0_ \mathcal{T} \subseteq \mathbb{V}^0_ {\mathcal{T}_*}$. By pivoting on this space, we proved an a posteriori error estimate which does not contain the stabilization term appearing in the VEM discrete formulation. Consequently, we established the convergence of the adaptive VEM algorithm, by a contraction argument.

Extensions of our work include:

● The complexity and optimality analysis of the two step algorithm AVEM mentioned in the Introduction to account for non-polynomial data;

● The study of a variant of the adaptive algorithm in which the polynomial degree $k$ may take large values, in the spirit of a $p$-version;

● The treatment of more general polygonal meshes which, as remarked in ^[5], seems non-trivial. The main difficulties lay in the choice of a suitable refinement strategy in replacement of the newest-vertex bisection used here, and in the lack of a conforming space $\mathbb{V}^0_ \mathcal{T}$ for general polygonal meshes.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The authors performed this research in the framework of the Italian MIUR Award "Dipartimenti di Eccellenza 2018-2022" granted to the Department of Mathematical Sciences, Politecnico di Torino (CUP: E11G18000350001). CC was partially supported by the Italian MIUR through the PRIN grant 201752HKH8; DF thanks the INdAM-GNCS project "Metodi numerici per lo studio di strutture geometriche parametriche complesse" (CUP: E53C22001930001). The authors are members of the Italian INdAM-GNCS research group.

Conflict of interest

The authors declare no conflicts of interest.