Directed acyclic graphs as conceptual and analytical tools in applied and theoretical epidemiology: advances, setbacks and future possibilities

1.   Introduction

In this paper, we consider the following diffusion equation on $\Omega\subset \mathbb{R}^2$,

To approximate (1), taking advantage of the adaptive mesh refinement (AMR) to save valuable computational resources, the adaptive finite element method on quadtree mesh is among the most popular ones in the engineering and scientific computing community [20]. Compared with simplicial meshes, quadtree meshes provide preferable performance in the aspects of the accuracy and robustness. There are lots of mature software packages (e.g., [1,2]) on quadtree meshes. To guide the AMR, one possible way is through the a posteriori error estimation to construct computable quantities to indicate the location that the mesh needs to be refined/coarsened, thus to balance the spacial distribution of the error which improves the accuracy per computing power. Residual-based and recovery-based error estimators are among the most popular ones used. In terms of accuracy, the recovery-based error estimator shows more appealing attributes [28,3].

More recently, newer developments on flux recovery have been studied by many researchers on constructing a post-processed flux in a structure-preserving approximation space. Using (1) as an example, given that the data $f\in L^2(\Omega)$, the flux $-\alpha \nabla u$ is in $\mathit{\boldsymbol{H}}(\mathrm{div}): = \{\mathit{\boldsymbol{v}}\in \mathit{\boldsymbol{L}}^2(\Omega): \nabla \cdot \mathit{\boldsymbol{v}}\in L^2(\Omega)\}$, which has less continuity constraint than the ones in [28,3] which are vertex-patch based with the recovered flux being $H^1(\Omega)$-conforming. The $\mathit{\boldsymbol{H}}(\mathrm{div})$-flux recovery shows more robustness than vertex-patch based ones (e.g., [11,10]).

However, these $\mathit{\boldsymbol{H}}(\mathrm{div})$-flux recovery techniques work mainly on conforming meshes. For nonconforming discretizations on nonmatching grids, some simple treatment of hanging nodes exists by recovering the flux on a conforming mother mesh [22]. To our best knowledge, there is no literature about the local $\mathit{\boldsymbol{H}}(\mathrm{div})$-flux recovery on a multilevel irregular quadtree meshes. One major difficulty is that it is impossible to recover a robust computable polynomial flux to satisfy the $\mathit{\boldsymbol{H}}(\mathrm{div})$-continuity constraint, that is, the flux is continuous in the normal direction on edges with hanging nodes.

More recently, a new class of methods called the virtual element methods (VEM) were introduced in [4,8], which can be viewed as a polytopal generalization of the tensorial/simplicial finite element. Since then, lots of applications of VEM have been studied by many researchers. A usual VEM workflow splits the consistency (approximation) and the stability of the method as well as the finite dimensional approximation space into two parts. It allows flexible constructions of spaces to preserve the structure of the continuous problems such as higher order continuities, exact divergence-free spaces, and many others. The VEM functions are represented by merely the degrees of freedom (DoF) functionals, not the pointwise values. In computation, if an optimal order discontinuous approximation can be computed elementwisely, then adding an appropriate parameter-free stabilization suffices to guarantee the convergence under common assumptions on the geometry of the mesh.

The adoption of the polytopal element brings many distinctive advantages, for example, treating rectangular element with hanging nodes as polygons allows a simple construction of $\mathit{\boldsymbol{H}}(\mathrm{div})$-conforming finite dimensional approximation space on meshes with multilevel irregularities. We shall follow this approach to perform the flux recovery for a conforming $\mathcal{Q}_k$ discretization of problem (1). Recently, arbitrary level of irregular quadtree meshes have been studied in [21,26,15]. Analyses of the residual-based error estimator on 1-irregular (balanced) quadtree mesh can be found, e.g., in [14]. In the virtual element context, Zienkiewicz-Zhu (ZZ)-type recovery techniques are studied for linear elasticity in [18], and for diffusion problems in [24]. In [18,24], the recovered flux is in $\mathit{\boldsymbol{H}}^1$ and associated with nodal DoFs, thus cannot yield a robust estimate when the diffusion coefficient has a sharp contrast [11,10]. The first equilibrated flux recovery in $\mathit{\boldsymbol{H}}(\mathrm{div})$ for virtual element methods is studied in [19]. While [19] recovers a flux by solving a mixed problem globally, we opt for a cheap and simple weighted averaging locally.

The major ingredient in our study is an $\mathit{\boldsymbol{H}}(\mathrm{div})$-conforming virtual element space modified from the ones used in [8,5] (Section 2.2). Afterwards, an $\mathit{\boldsymbol{H}}(\mathrm{div})$-conforming flux is recovered by a robust weighted averaging of the numerical flux, in which some unique properties of the tensor-product type element $\mathcal{Q}_k$ are exploited (Section 3). The a posteriori error estimator is constructed based on the projected flux elementwisely. The efficiency of the local error indicator is then proved by bounding it above by the residual-based error indicator (Section 4.1). The reliability of the recovery-based error estimator is then shown under certain assumptions (Section 4.2). These estimates are verified numerically by some common AMR benchmark problems implemented in a publicly available finite element software library $i$FEM [16] (Section 5).

2.   Preliminaries

2.1.   Discretization and notations

If $\Omega$ is not a rectangle, $u$ is extended by 0 to an $\widetilde{\Omega}$ that is rectangular, therefore without loss of generality, we assume $\Omega$ is partitioned into a shape-regular $\mathcal{T} = \{K\}$ with rectangular elements, and $\alpha : = \alpha_K$ is assumed to be a piecewise, positive constant with respect to $\mathcal{T}$. The weak form of problem (1) is then discretized in a tensor-product finite element space as follows,

in which the standard notation is opted. $(\cdot,\cdot)_D$ denotes the inner product on $L^2(D)$, and $\left\Vert{\cdot}\right\Vert_D: = \sqrt{(\cdot,\cdot)_D}$, with the subscript omitted when $D = \Omega$. The discretization space is

and on $K = [a,b]\times [c,d]$

where $\mathbb{P}_k(D)$ stands for the degree no more than $k$ polynomial defined on $D$. Henceforth, we shall simply denote $\mathcal{Q}_k(\mathcal{T}) = : \mathcal{Q}_k$ when no ambiguity arises.

On $K$, the sets of 4 vertices, as well as 4 edges of the same generation with $K$, are denoted by $\mathcal{N}_K$ and $\mathcal{E}_K$, respectively. The sets of nodes and edges in $\mathcal{T}$ are denoted by $\mathcal{N}: = \bigcup_{K\in \mathcal{T}} \mathcal{N}_K$ and $\mathcal{E}: = \bigcup_{K\in \mathcal{T}} \mathcal{E}_K$. A node $\mathit{\boldsymbol{z}} \in \mathcal{N}$ is called a hanging node if it is on $\partial K$ but is not counted as a vertex of $K\in \mathcal{T}$, and we denote the set of hanging nodes as $\mathcal{N}_H$

Otherwise the node $\mathit{\boldsymbol{z}} \in \mathcal{N}$ is a regular node. If an edge $e\in \mathcal{E}$ contains at most $l$ hanging nodes, the partition $\mathcal{T}$, as well as the element these hanging nodes lie on, is called $l$-irregular.

For each edge $e\in \mathcal{E}$, a unit normal vector $\mathit{\boldsymbol{n}}_e$ is fixed by specifying its direction pointing rightward for vertical edges, and upward for horizontal edges. If an exterior normal of an element on this edge shares the same orientation with $\mathit{\boldsymbol{n}}_e$, then this element is denoted by $K_-$, otherwise it is denoted by $K_+$, i.e., $\mathit{\boldsymbol{n}}_e$ is pointing from $K_-$ to $K_+$. The intersection of the closures of $K_+, K_-$ is always an edge $e\in \mathcal{E}$. However, we note that by the definition in (3) it is possible that $e\in\mathcal{E}_{K_+}$ but not in $\mathcal{E}_{K_-}$ or vice versa, if there exists a hanging node on $e$ (see e.g., Figure 1). For any function or distribution $v$ well-defined on the two elements, define $\lbrack\lbrack {{v}} \rbrack\rbrack_{e} = v^- - v^+$ on an edge $e\not\subset\partial \Omega$, in which $v^-$ and $v^+$ are defined in the limiting sense $v^{\pm} = \lim_{\epsilon\to 0^{\pm}}v(\mathit{\boldsymbol{x}}+\epsilon \mathit{\boldsymbol{n}}_e)$ for $\mathit{\boldsymbol{x}}\in e$. If $e$ is a boundary edge, the function $v$ is extended by zero outside the domain to compute $\lbrack\lbrack {{v}} \rbrack\rbrack_{e}$. Furthermore, the following notation denotes a weighted average of $v$ on edge $e$ for a weight $\gamma\in [0,1]$,

2.2.   Virtual element spaces

In this subsection, the quadtree mesh $\mathcal{T}$ of interest is embedded into a polygonal mesh $\mathcal{T} \hookrightarrow \mathcal{T}_{\mathrm{poly}} = \{K_{\mathrm{poly}}\}$. On any given quadrilateral element $K$, for example we consider a $v_{\mathcal{T}} \in \mathbb{Q}_1(K)$, it has 4 degrees of freedom associated with 4 nodes $\{z\}$. Its numerical flux $-\alpha \nabla v_{\mathcal{T}} \cdot \mathit{\boldsymbol{n}}$ is well-defined on the 4 edges $\{e\}$ locally on $K$, such that on each edge it is a polynomial defined on the whole edge, regardless of the number of hanging nodes on that edge. Using Figure 1 as an example, on the upper right element $K$, $\nabla {v_{\cal T}}{|_K} \cdot {\boldsymbol{n}}|\mathop \rightharpoonup \limits_{{z_2}{z_6}} \in {\mathbb{P}_1}(\mathop \rightharpoonup \limits_{{z_2}{z_6}} )$ is a linear function in $y$-variable.

For the embedded element $K_{\mathrm{poly}} \in \mathcal{T}_{\mathrm{poly}}$, which geometrically coincides with $K$, it includes all the hanging nodes, while the set of edges are formed accordingly as the edges of the cyclic graph of the vertices. We shall denote the set of all edges on $\mathcal{T}_{\mathrm{poly}}$ as $\mathcal{E}_{\mathrm{poly}}$. Using Figure 1 as example, it is possible to define a flux on $K$ with piecewise linear normal component on $\mathop \rightharpoonup \limits_{{z_2}{z_6}}$ which now consists of three edges on $\partial K_{\mathrm{poly}}$.

Subsequently, $K_{\mathrm{poly}}\in \mathcal{T}_{\mathrm{poly}}$ shall be denoted by simply $K\in \mathcal{T}_{\mathrm{poly}}$ in the context of flux recovery, and the notion $e\subset \partial K$ denotes an edge on the boundary of $K$, which takes into account of the edges formed with one end point or both end points as the hanging nodes.

On $\mathcal{T}_{\mathrm{poly}}$, we consider the following Brezzi-Douglas-Marini-type virtual element modification inspired by the ones used in [8,5]. The local space on a $K\in \mathcal{T}_{\mathrm{poly}}$ is defined as for $k\geq 1$

An $\mathit{\boldsymbol{H}}(\mathrm{div})$-conforming global space for recovering the flux is then

Next we turn to define the degrees of freedom (DoFs) of this space. To this end, we define the set of scaled monomials $\mathbb{P}_{k}(e)$ on an edge $e$. $e$ is parametrized by $[0,h_e]\ni s\mapsto \mathit{\boldsymbol{a}} + s\mathit{\boldsymbol{t}}_e$, where $\mathit{\boldsymbol{a}}$ is the starting point of $e$, and $\mathit{\boldsymbol{t}}_e$ is the unit tangential vector of $e$. The basis set for $\mathbb{P}_{k}(e)$ is chosen as:

where $m_e = h_e/2$ representing the midpoint when using this parametrization. Similar to the edge case, $\mathbb{P}_{k}({K})$'s basis set is chosen as follows (see e.g., [4]):

The degrees of freedom (DoFs) are then set as follows for a $\mathit{\boldsymbol{\tau}}\in \mathcal{V}_{k}$:

Remark 1. We note that in our construction, the degrees of freedom to determine the curl of a VEM function originally in [8] are replaced by a curl-free constraint thanks to the flexibility to virtual element. The reason why we opt for this subspace is that the true flux $-\alpha \nabla u$ is locally curl-free since we have assumed that $\alpha$ is a piecewise constant. The unisolvency of the set of DoFs (8) including the curl-part can be found in [8]. While for the modified space (4), a simplified argument is in the proof of Lemma 7.3.

3.   Flux recovery

As the data $f\in L^2(\Omega)$, the true flux $\mathit{\boldsymbol{\sigma}} = -\alpha\nabla u\in \mathit{\boldsymbol{H}}(\mathrm{div})$. Consequently, we shall seek a postprocessed flux $\mathit{\boldsymbol{\sigma}}_{\mathcal{T}}$ in $\mathcal{V}_{k}\subset \mathit{\boldsymbol{H}}(\mathrm{div})$ by specifying the DoFs in (8). Throughout this section, whenever considering an element $K\in \mathcal{T}$, we treat it a polygon as $K\in \mathcal{T}_{\mathrm{poly}}$.

3.1.   Virtual element-based flux recovery

Consider $-\alpha_K \nabla u_{\mathcal{T}}$ which is the numerical flux on $K$. We note that $-\alpha_K \nabla u_{\mathcal{T}}|_K \in \mathbb{P}_{k-1,k}(K) \times\mathbb{P}_{k,k-1}(K)$. The normal flux on each edge $e\in \mathcal{E}_{\mathrm{poly}}$ is in $\mathbb{P}_{k}(e)$ as $n_e = (\pm 1, 0)$ and $x = \mathrm{const}$ on vertical edges, $n_e = (0, \pm 1)$ and $y = \mathrm{const}$ on horizontal edges. Therefore, the edge-based DoFs can be computed by a simple averaging thanks to the matching polynomial degrees of the numerical flux to the functions in $\mathcal{V}_k$.

On each $e = \partial K_+ \cap \partial K_-$, define

where

First for both $k = 1$ and $k\geq 2$ cases, we set the normal component of the recovered flux is set as

In the lowest order case $k = 1$, $\nabla \cdot \mathit{\boldsymbol{\sigma}}_{\mathcal{T}}$ is a constant on $K$ by (4), thus the construction (11) alone, which consists the edge DoFs $(\mathfrak{e})$ in (8), can determine the divergence $\nabla \cdot \mathit{\boldsymbol{\sigma}}_{\mathcal{T}}$ in $K$ as follows

If $k\geq 2$, after the normal component (11) is set, furthermore on each $K$, denote $\Pi_{k-1}$ stands for the $L^2$-projection to $\mathbb{P}_{k-1}(K)$, and we let

The reason to add $c_K$ is that we have set the normal components of the recovered flux first without relying on the divergence information. While in general $\nabla \cdot \mathit{\boldsymbol{\sigma}}_{\mathcal{T}} \neq \Pi_{k-1} f$ as otherwise the divergence theorem will be rendered invalid in (12). As a result, an element-wise constant $c_K$ is added to ensure the compatibility of $\mathit{\boldsymbol{\sigma}}_{\mathcal{T}}$ locally on each $K$. It is straightforward to verify that $c_K$ has the following form, and later we shall show that $c_K$ does not affect the efficiency as well as the reliability of the error estimates.

Consequently for $k\geq 2$, the set $(\mathfrak{i})$ of DoFs can be set as: $\forall q\in \mathbb{P}_{k-1}(K)$

3.2.   Locally projected flux

To the end of constructing a computable local error indicator, inspired by the VEM formulation [8], the recovered flux is projected to a space with a much simpler structure. A local oblique projection ${\Pi}: \mathit{\boldsymbol{L}}^2(K) \to \nabla \mathbb{P}_k(K), \; \mathit{\boldsymbol{\tau}}\mapsto {\Pi} \mathit{\boldsymbol{\tau}}$ is defined as follows:

Next we are gonna show that this projection operator can be straightforward computed for vector fields in $\mathcal{V}_k(K)$.

3.2.1.   $k=1$

When $k = 1$, we can compute the right hand side of (16) as follows:

By definition of the space (4) when $k = 1$, $\nabla \cdot \mathit{\boldsymbol{\tau}}$ is a constant on $K$ and can be determined by edge DoFs $(\mathfrak{e})$ in (8) similar to (12). Moreover, $p|_e\in \mathbb{P}_1(e)$, thus the boundary term can be evaluated using DoFs $(\mathfrak{e})$ in (8).

3.2.2.   $k\geq 2$

When $k\geq 2$, the right hand side of (16) can be evaluated following a similar procedure as (17), if we exploit the fact that $\nabla \cdot \mathit{\boldsymbol{\tau}}\in \mathbb{P}_{k-1}(K)$, we have

which can be evaluated using both DoF sets $(\mathfrak{e})$ and $(\mathfrak{i})$.

4.   A posteriori error estimation

Given the recovered flux \(\mathit{\boldsymbol{\sigma}}_{\mathcal{T}}\) in Section 3, the recovery-based local error indicator $\eta_{\mathrm{flux},K}$ and the element residual $\eta_{\mathrm{res},K}$ as follows:

then

A computable $\widehat{\eta}_{\mathrm{flux},K}$ is defined as:

with the oblique projection ${\Pi}$ defined in (16). The stabilization part $\widehat{\eta}_{\mathrm{stab},K}$ is

Here $|{\cdot}|_{S,K}: = \big(S_K(\cdot, \cdot)\big)^{1/2}$ is seminorm induced by the following stabilization

where $\Lambda$ is the index set for the monomial basis of $\mathbb{P}_{k-1}(K)/\mathbb{R}$ with cardinality $k(k+1)/2 - 1$, i.e., the second term in (23) is dropped in the $k = 1$ case. We note that this is a slightly modified version of the standard stabilization for an $\mathit{\boldsymbol{H}}(\mathrm{div})$-function in [8] as we have replaced the edge DoFs by an integral. In Section 7.1 it is shown that the integral-based stabilization still yields the crucial norm equivalence result.

The computable error estimator $\widehat{\eta}$ is then

4.1.   Efficiency

In this section, we shall prove the proposed recovery-based estimator $\widehat{\eta}_{K}$ is efficient by bounding it above by the residual-based error estimator. In the process of adaptive mesh refinement, only the computable $\widehat{\eta}_{K}$ is used as the local error indicator to guide a marking strategy of choice.

Theorem 4.1. Let $u_{\mathcal{T}}$ be the solution to problem (2), and $\widehat{\eta}_{\mathrm{flux},K}$ be the error indicator in (24). On $K\in \mathcal{T}_{\mathrm{poly}}$, $\widehat{\eta}_{\mathrm{flux},K}$ can be locally bounded by the residual-based ones:

where

In the edge jump term, $K_e$ is the element on the opposite side of $K$ with respect to an edge $e\subset \partial K$. The constant depends on $k$ and the number of edges on $\partial K$.

Proof. Let $\alpha^{-1}_K{\Pi}(\mathit{\boldsymbol{\sigma}}_{\mathcal{T}} + \alpha_K \nabla u_{\mathcal{T}}) = :\nabla p$ on $K$, then $p\in \mathbb{P}_k(K)/\mathbb{R}$ and we have

By (11), without loss of generality we assume $K = K_-$ (the local orientation of $e$ agrees with the global one, i.e., $\mathit{\boldsymbol{n}}_{\partial K}\big|_{e} = \mathit{\boldsymbol{n}}_e$), and $K_e = K_+$ which is the element opposite to $K$ with respect to $e$, and $\gamma_e : = {\alpha_{K_e}^{1/2}}/({\alpha_{K_e}^{1/2} + \alpha_{K}^{1/2}})$, we have on edge $e\subset\partial K$

The boundary term in (26) can be then rewritten as

By a trace inequality on an edge of a polygon (Lemma 7.1), and the Poincaré inequality for $p\in \mathbb{P}_k(K)/\mathbb{R}$, we have,

As a result,

For the bulk term on $K$'s in (26), when $k = 1$, by (12), the representation in (28), and the Poincaré inequality for $p\in \mathbb{P}_k(K)/\mathbb{R}$ again with $h_K\simeq |K|^{1/2}$, we have

When $k\geq 2$, by (13),

The first two terms can be handled by combining the weights $\alpha^{-1/2}$ and $h_K$ from $\left\Vert{p}\right\Vert_K\leq h_K \left\Vert{\nabla p}\right\Vert_K$. For $c_K$, it can be estimated straightforwardly as follows

The two terms on $K$ can be treated the same way with the first two terms in (29) while the edge terms are handled similarly as in the $k = 1$ case. As a result, we have shown

and the theorem follows.

Theorem 4.2. Under the same setting with Theorem 4.1, let $\widehat{\eta}_{\mathrm{stab},K}$ as the estimator in (22), we have

The constant depends on $k$ and the number of edges on $\partial K$.

Proof. This theorem follows directly from the norm equivalence Lemma 7.3:

while evaluating the DoFs $(\mathfrak{e})$ and $(\mathfrak{i})$ using (11) and (15) reverts us back to the proof of Theorem 4.1.

Theorem 4.3. Under the same setting with Theorem 4.1, on any $K\in \mathcal{T}_{\mathrm{poly}}$ with $\omega_K$ defined as the collection of elements in $\mathcal{T}$ which share at least 1 vertex with $K$

with a constant independent of $\alpha$, but dependent on $k$ and the maximum number of edges in $K\in \mathcal{T}_{\mathrm{poly}}$.

Proof. This is a direct consequence of Theorem 4.1 and 4.2 and the fact that the residual-based error indicator is efficient by a common bubble function argument.

4.2.   Reliability

In this section, we shall prove that the computable error estimator $\widehat{\eta}$ is reliable under two common assumptions in the a posteriori error estimation literature. For the convenience of the reader, we rephrase them here using a "layman" description, for more detailed and technical definition please refer to the literature cited.

Assumption 1 ($\mathcal{T}$ is $l$-irregular [14]). Any given $\mathcal{T}$ is always refined from a mesh with no hanging nodes by a quadsecting red-refinement. For any two neighboring elements in $\mathcal{T}$, the difference in their refinement levels is $\leq l$ for a uniformly bounded constant $l$, i.e., for any edge $e\in \mathcal{E}$, it has at most $l$ hanging nodes.

By Assumption 1, we denote the father $1$-irregular mesh of $\mathcal{T}$ as $\mathcal{T}_1$. On $\mathcal{T}_1$, a subset of all nodes is denoted by $\mathcal{N}_{1}$, which includes the regular nodes $\mathcal{N}_R$ on $\mathcal{T}_1$, as well as $\mathcal{N}_E$ as the set of end points of edges with a hanging node as the midpoint. By [14,Theorem 2.1], there exists a set of bilinear nodal bases $\{\phi_z\}$ associated with $z\in \mathcal{N}_{1}$, such that $\{\phi_z\}$ form a partition of unity and can be used to construct a Clément-type quasi-interpolation. Furthermore, the following assumption assures that the Clément-type quasi-interpolant is robust with respect to the coefficient distribution on a vertex patch, when taking nodal DoFs as a weighted average.

Assumption 2 (Quasi-monotonicity of $\alpha$ [6]). On $\mathcal{T}$, let $\phi_z$ be the bilinear nodal basis associated with $z\in \mathcal{N}_{1}$, with $\omega_{z} : = \operatorname{supp} \phi_z$. For every element $K \subset \omega_{z}, K\in \mathcal{T}$, there exists a simply connected element path leading to $\omega_{m(z)}$, which is a Lipschitz domain containing the elements where the piecewise constant coefficient $\alpha$ achieves the maximum (or minimum) on $\omega_{z}$.

Denote

We note that if $\alpha$ is a constant on $\omega_z$, $(1, \left(v-\pi_{z} v\right) \phi_{z})_{\omega_z} = 0$. A quasi-interpolation $\mathcal{I}: L^2(\Omega) \to \mathcal{Q}_1(\mathcal{T}_1)$ can be defined as

Lemma 4.4 (Estimates for $\pi_z$ and $\mathcal{I}$). Under Assumption 1 and 2, the following estimates hold for any $v\in H^1(\omega_K)$

and for $\mathit{\boldsymbol{z}}\in \mathcal{N}_1$

in which $h_z : = \max_{K\subset\omega_z} h_K$, and here $\omega_K$ denotes the union of elements in $\mathcal{T}_1$ sharing at least a node (hanging or regular) with $K$.

Proof. The estimate for $\pi_z$ follows from [6,Lemma 2.8]. For $\mathcal{I}$, its error estimates and stability only rely on the partition of unity property of the nodal basis set $\{\phi_z\}$ (see e.g., [27]), therefore the proof follows the same argument with the ones used on triangulations in [6,Lemma 2.8].

Denotes the subset of nodes $\{\mathit{\boldsymbol{z}}\}\subset \mathcal{N}_1$ (i) on the boundary as $\mathcal{N}_{\partial \Omega}$ and (ii) with the coefficient $\alpha$ on patch $\omega_z$ as $\mathcal{N}_I$. For the lowest order case, we need the following oscillation term for $f$

with $f_z : = {\int_{\omega_{z}} v \phi_z }/{\int_{\omega_{z}} \phi_z }$.

Theorem 4.5. Let $u_{\mathcal{T}}$ be the solution to problem (2), and $\widehat{\eta}$ be the computable error estimator in (24), under Assumption 2 and 1, we have for $k = 1$

For $k\geq 2$,

where the constant depends on $l$ and $k$.

Proof. Let $\varepsilon : = u - u_{\mathcal{T}} \in H^1_0(\Omega)$, and $\mathcal{I}\varepsilon\in \mathcal{Q}_1(\mathcal{T}_1) \subset \mathcal{Q}_1(\mathcal{T})$ be the quasi-interpolant in (34) of $\varepsilon$, then by the Galerkin orthogonality, $\alpha\nabla u + \mathit{\boldsymbol{\sigma}}_{\mathcal{T}} \in \mathit{\boldsymbol{H}}(\mathrm{div})$, the Cauchy-Schwarz inequality, and the interpolation estimates (35), we have for $k\geq 2$,

Applying the norm equivalence of $\eta$ to $\widehat{\eta}$ by Lemma 7.3, as well as the fact that the number of elements in $\omega_K$ is uniformly bounded by Assumption 1, yields the desired estimate.

When $k = 1$, the residual term on $K$ can be further split thanks to $\Delta \mathbb{Q}_1(K) = \{0\}$. First we notice that by the fact that $\{\phi_z\}$ form a partition of unity,

in which a patch-wise constant $f_z$ (weighted average of $f$) can be further inserted by the definition of $\pi_z$ (33) if $\alpha$ is a constant on $\omega_z$. Therefore, by the assumption of $\alpha_K$ being a piecewise constant, splitting (40), we have

Applied an inverse inequality in Lemma 7.2 on $\big\Vert \nabla\cdot (\mathit{\boldsymbol{\sigma}}_{\mathcal{T}} + \alpha_K\nabla u_{\mathcal{T}}) \big\Vert_{K}$ and the projection estimate for $\pi_z$ (36), the rest follows the same argument with the one used in the $k\geq 2$ case.

5.   Numerical examples

The numerics is prepared using the bilinear element for common AMR benchmark problems. The codes for this paper are publicly available on https://github.com/lyc102/ifem implemented using $i$FEM [16]. The linear algebraic system on an $l$-irregular quadtree is implemented following the conforming prolongation approach [15] by $\mathbf{P}^{\top}\mathbf{A} \mathbf{P}\mathbf{u} = \mathbf{P}^{\top}\mathbf{f}$, where $\mathbf{A}$ is the locally assembled stiffness matrix for all nodes in $\mathcal{N}$, $\mathbf{u}$ and $\mathbf{f}$ are the solution vector associated with $\mathcal{N}_R$ and load vector associated with $\mathcal{N}$, respectively. $\mathbf{P} = (\mathbf{I}, \mathbf{W})^{\top}: \mathbb{R}^{\dim \mathcal{N}_R}\to \mathbb{R}^{\dim \mathcal{N}}$ is a prolongation operator mapping conforming $H^1$-bilinear finite element function defined on regular nodes to all nodes, the weight matrix $\mathbf{W}$ is assembled locally by a recursive $k$NN query in $\mathcal{N}_H$, while the polygonal mesh data structure embedding is automatically built during constructing $\mathbf{P}$. For details we refer the readers to https://github.com/lyc102/ifem/tree/master/research/polyFEM.

The adaptive finite element (AFEM) iterative procedure is following the standard

The linear system is solved by MATLAB $\texttt{mldivide}$. In $\texttt{MARK}$, the Dorfler $L^2$-marking is used with the local error indicator $\widehat{\eta}_{K}$ in that the minimum subset $\mathcal{M} \subset \mathcal{T}$ is chosen such that

Throughout all examples, we fix $\theta = 0.3$. $\mathcal{T}$ is refined by a red-refinement by quadsecting the marked element afterwards. For comparison, we compute the standard residual-based local indicator for $K\in \mathcal{T}_{\mathrm{poly}}$

Let $\eta_{\;{\rm{Residual}}\;}^2 = \sum_{K\in \mathcal{T}} \eta_{\;{\rm{Residual}}\;,K}^2$. The residual-based estimator $\eta_{\;{\rm{Residual}}\;}$ is merely computed for comparison purpose and not used in marking. The AFEM procedure stops when the relative error reaches a threshold. The effectivity indices for different estimators are compared

i.e., the closer to 1 the effectivity index is, the more accurate this estimator is to measure the error of interest. We use an order $5$ Gaussian quadrature to compute $\Vert{\alpha^{1/2} \nabla (u - u_{\mathcal{T}}) } \Vert$ elementwisely. The orders of convergence for various $\eta$'s and $\Vert{\alpha^{1/2} \nabla (u - u_{\mathcal{T}}) } \Vert$ are computed, for which $r_{\eta}$ and $r_{\;{\rm{err}}\;}$ are defined as the slope for the linear fitting of $\ln \eta_n$ and $\ln \Vert\alpha^{1/2}\nabla(u - u_{\mathcal{T},n}) \Vert$ in the asymptotic regime,

where the subscript $n$ stands for the number of iteration in the AFEM cycles, $N_n: = \# ( \mathcal{N}_R\backslash \mathcal{N}_{\partial \Omega})$. $r_{\eta}$ and $r_{\;{\rm{err}}\;}$ are considered optimal when being close to $1/2$.

5.1.   L-shaped domain

In this example, a standard AMR benchmark on the L-shaped domain is tested. The true solution $u = r^{2/3}\sin(2\theta/3)$ in polar coordinates on $\Omega = (-1,1) \times(-1,1) \backslash[0,1) \times(-1,0]$. The AFEM procedure stops if the relative error has reached $0.01$. The adaptively refined mesh can be found in Figure 2a. While both estimators show optimal rate of convergence in Figure 2b, the effectivity index for $\eta_{\mathrm{Residual}}$ is $4.52$, and is $2.24$ for $\widehat{\eta}$.

5.2.   A circular wave front

The solution $u = \tan^{-1}(\alpha(r-r_0))$ is defined on $\Omega = (0,1)^2$ with $r: = \sqrt{(x+0.05)^2+(y+0.05)^2}$, $\alpha = 100$, and $r_0 = 0.7$. The true solution shows a sharp transition layer (Figure 3a). The result of the convergence can be found in Figure 3b. In this example, the AFEM procedure stops if the relative error has reached $0.05$. Additionally, we note that by allowing $l$-irregular ($l\geq 2$), the AMR procedure shows to be more efficient toward capturing the singularity of the solution. A simple comparison can be found in Figure 4. The effectivity indices for $\eta_{\mathrm{Residual}}$ and $\widehat{\eta}$ are $5.49$ and $2.08$, respectively.

5.3.   Kellogg benchmark

This example is a common benchmark test problem introduced in [9], see also [17,12]) for elliptic interface problems. The true solution $u = r^{\gamma}\mu(\theta)$ is harmonic in four quadrants, and $\mu(\theta)$ takes different values within four quadrants:

While $\alpha = R$ in the first and third quadrants, and $\alpha = 1$ in the second and fourth quadrants, and the true flux $\alpha \nabla u$ is glued together using $\mathit{\boldsymbol{H}}(\mathrm{div})$-continuity conditions. We choose the following set of coefficients for $u$

By this choice, this function is very singular near the origin as the maximum regularity it has is $H^{1+\gamma}_{loc}(\Omega\backslash\{\mathit{\boldsymbol{0}}\})$. Through an integration by parts, it can be computed accurately that $\|\alpha^{1/2}\nabla u\| \approx 0.56501154$. For detailed formula and more possible choices of the parameters above, we refer the reader to [17].

The AFEM procedure for this problem stops when the relative error reaches $0.05$, and the resulting mesh and finite element approximation during the refinement can be found in Figure 5, and the AFEM procedure shows optimal rate of convergence in Figure 6. The effectivity index for $\eta_{\mathrm{Residual}}$ is $2.95$, and $1.33$ for $\widehat{\eta}$.

6.   Conclusion

A postprocessed flux with the minimum $\mathit{\boldsymbol{H}}(\mathrm{div})$ continuity requirement is constructed for tensor-product type finite element. The implementation can be easily ported to finite element on quadtree to make use the vast existing finite element libraries in the engineering community. Theoretically, the local error indicator is efficient, and the global estimator is shown to be reliable under the assumptions that (i) the mesh has bounded irregularities, and (ii) the diffusion coefficient is a quasi-monotone piecewise constant. Numerically, we have observed that both the local error indicator and the global estimator are efficient and reliable (in the asymptotic regime), respectively. Moreover, the recovery-based estimator is more accurate than the residual-based one.

However, we do acknowledge that the technical tool involving interpolation is essentially limited to $1$-irregular meshes in reliability. A simple weighted averaging has restrictions and is hard to generalize to $hp$-finite elements, or discretization on curved edges/isoparametric elements. Nevertheless, we have shown that the flexibility of the virtual element framework allows further modification of the space in which we perform the flux recovery to cater the needs.

Acknowledgments

The author is grateful for the constructive advice from the anonymous reviewers. This work was supported in part by the National Science Foundation under grants DMS-1913080 and DMS-2136075, and no additional revenues are related to this work.

7.   Appendix

7.1.   Inverse estimates and the norm equivalence of a virtual element function

Unlike the identity matrix stabilization commonly used in most of the VEM literature, for $\mathit{\boldsymbol{\tau}}\in \mathcal{V}_k(K)$, we opt for a mass matrix/DoF hybrid stabilizer approach. Let $\big\Vert{ \alpha^{-1/2}\mathit{\boldsymbol{\tau}}}\big\Vert_{h,K}^2 : = (\!(\mathit{\boldsymbol{\tau}}, {\mathit{\boldsymbol{\tau}}})\!)_{K}$ and

where $S_{K}(\cdot,\cdot)$ is defined in (23).

To show the inverse inequality and the norm equivalence used in the reliability bound, on each element, we need to introduce some geometric measures. Consider a polygonal element $K$ and an edge $e\subset \partial K$, let the height $l_e$ which measures how far from this edge $e$ one can advance to an interior subset of $K$, and denote $T_e\subset K$ as a right triangle with height $l_e$ and base as edge $e$.

Proposition 1. Under Assumption 1, $\mathcal{T}_{poly}$ satisfies (1) The number of edges in every $K\in \mathcal{T}_{poly}$ is uniformly bounded above. (2) For any edge $e$ on every $K$, $l_e/h_e$ is uniformly bounded below.

Lemma 7.1 (Trace inequality on small edges [13]). If Proposition 1 holds, for $v \in H^1(K)$ and $K\in \mathcal{T}_{\mathrm{poly}}$ we have

Proof. The proof follows essentially equation (3.9) in [13,Lemma 3.3] as a standard scaled trace inequality on $e$ toward $T_e$ reads

Lemma 7.2 (Inverse inequalities). Under Assumption 1, we have the following inverse estimates for $\mathit{\boldsymbol{\tau}} \in \mathcal{V}_k(K)$ (4) on any $K\in \mathcal{T}_{poly}$ with constants depending on $k$ and the number of edges in $K$:

Proof. The first inequality in (43) can be shown using a bubble function trick. Choose $b_K$ be a bubble function of $T_{e'}$ where $e'$ is the longest edge on $\partial K$. Denote $p : = \nabla \cdot \mathit{\boldsymbol{\tau}} \in \mathbb{P}_{k-1}(K)$, we have

and then $\left\Vert{ \nabla (p b_K)}\right\Vert$ can be estimated as follows

Consequently, the first inequality in (43) follows above by the standard inverse estimate for polynomials $\left\Vert{\nabla p}\right\Vert_K\lesssim h_K^{-1} \left\Vert{p}\right\Vert_K$, and the properties of the bubble function $\left\Vert{b_K}\right\Vert_{\infty, K} = O(1)$, and $\left\Vert{\nabla b_K}\right\Vert_{ \infty, K} = O(h_K^{-1})$.

To prove the second inequality in (43), by integration by parts we have

Expand $\nabla\cdot \mathit{\boldsymbol{\tau}} = p$ in the monomial basis $p(\mathit{\boldsymbol{x}}) = \sum_{\alpha\in \Lambda} p_{\alpha} m_{\alpha}(\mathit{\boldsymbol{x}})$, and denote the mass matrix $\mathbf{M}: = \big( (m_{\alpha}, m_{\gamma})_{K} \big)_{\alpha\gamma}$, $\mathbf{p}: = (p_{\alpha})_{\alpha\in \Lambda}$, it is straightforward to see that

since $\int_K (x-x_K)^l (y-y_K)^m\,\mathrm{d} x\mathrm{d} y\geq 0$ for the off-diagonal entries of $\mathbf{M}$ due to $K$ being geometrically a rectangle (with additional vertices). As a result, applying the trace inequality in Lemma 7.1 on (44) yields

As a result, the second inequality in (43) is proved when apply an inverse inequality for $\left\Vert{\nabla p}\right\Vert_K$ and estimate (45).

Remark 2. While the proof in Lemma 7.2 relies on $K$ being a rectangle, the result holds for a much broader class of polygons by changing the basis of $\mathbb{P}_{k-1}(K)$ from the simple scaled monomials to quasi-orthogonal ones in [25,7] and apply the isotropic polygon scaling result in [13].

Lemma 7.3 (Norm equivalence). Under Assumption 1, let ${\Pi}$ be the oblique projection defined in (16), then the following relations holds for $\mathit{\boldsymbol{\tau}} \in \mathcal{V}_k(K)$ (4) on any $K\in \mathcal{T}_{poly}$:

where both $\gamma_*$ and $\gamma^*$ depends on $k$ and the number of edges in $K$.

Proof. First we consider the lower bound, by triangle inequality,

Since ${\Pi}\mathit{\boldsymbol{\tau}} \in \mathcal{V}^{k}(K)$, it suffices to establish the following to prove the lower bound in (46)

To this end, we consider the weak solution to the following auxiliary boundary value problem on $K$:

By a standard Helmholtz decomposition result (e.g. Proposition 3.1, Chapter 1[23]), we have $\mathit{\boldsymbol{\tau}} -\nabla \psi = \nabla^{\perp} \phi$. Moreover, since on $\partial K$, $0 = \nabla^{\perp} \phi \cdot \mathit{\boldsymbol{n}} = \nabla \phi \cdot \mathit{\boldsymbol{t}} = \partial \phi/\partial s$, we can further choose $\phi\in H^1_0(K)$. As a result, by the assumption that $\nabla\times \mathit{\boldsymbol{\tau}} = 0$ for $\mathit{\boldsymbol{\tau}}$ in the modified virtual element space (4), we can verify that

Consequently, we proved essentially the unisolvency of the modified VEM space (4) and $\mathit{\boldsymbol{\tau}} = \nabla \psi$. We further note that $\psi$ in (48) can be chosen in $H^1(K)/\mathbb{R}$ and thus

Proposition 1 allows us to apply an isotropic trace inequality on an edge of a polygon (Lemma 7.1), combining with the Poincaré inequality for $H^1(K)/\mathbb{R}$, we have, on every $e\subset \partial K$,

Furthermore applying the inverse estimate in Lemma 7.2 on the bulk term above, we have

which proves the validity of (47), thus yield the lower bound.

To prove the upper bound, by $\big\Vert{{\Pi}\mathit{\boldsymbol{\tau}}}\big\Vert_{K}\leq \Vert{\mathit{\boldsymbol{\tau}}}\Vert_{K}$, it suffices to establish the reversed direction of (47) on a single edge $e$ and for a single monomial basis $m_{\alpha}\in \mathbb{P}_{k-1}(K)$:

To prove the first inequality, by Proposition 1 again, consider the edge bubble function $b_e$ such that $\operatorname{supp} b_e = T_e$. We can let $b_e = 0$ on $e'\subset\partial K$ for $e'\neq e$. It is easy to verify that:

Denote $q_e: = \mathit{\boldsymbol{\tau}}\cdot\mathit{\boldsymbol{n}}_e$, and extend it to $\mathop K\limits^ \circ$ by a constant extension in the normal direction rectangular strip $R_e \subset K$ with respect to $e$ (notice $\operatorname{supp} b_e \subset R_e$), we have

Now by the fact that $\left\Vert{q_e}\right\Vert_{T_e} \lesssim h_e^{1/2} \left\Vert{q_e}\right\Vert_e$, the scaling of the edge bubble function in (51), and the first inverse estimate of $\left\Vert{\nabla\cdot\mathit{\boldsymbol{\tau}}}\right\Vert_K\lesssim h_K^{-1}\left\Vert{\mathit{\boldsymbol{\tau}}}\right\Vert_K$ in Lemma 7.2 yields the first part of (50).

The second inequality in (50) can be estimated straightforward by the scaling of the monomials (7)

Hence, (46) is proved.