# A virtual element-based flux recovery on quadtree

In this paper, we introduce a simple local flux recovery for 𝒬_k finite element of a scalar coefficient diffusion equation on quadtree meshes, with no restriction on the irregularities of hanging nodes. The construction requires no specific ad hoc tweaking for hanging nodes on l-irregular (l≥ 2) meshes thanks to the adoption of some novel ideas borrowed from virtual element families. The rectangular elements with hanging nodes are treated as polygons as in the flux recovery context. An efficient a posteriori error estimator is then constructed based on the recovered flux projected to a space with simpler structure, and its reliability is proved under common assumptions, both of which are further verified in numerics.

## Authors

• 7 publications
• ### Residual-based a posteriori error estimation for immersed finite element methods

In this paper we introduce and analyze the residual-based a posteriori e...
10/17/2019 ∙ by Cuiyu He, et al. ∙ 0

• ### Lower a posteriori error estimates on anisotropic meshes

Lower a posteriori error bounds obtained using the standard bubble funct...
06/13/2019 ∙ by Natalia Kopteva, et al. ∙ 0

• ### Superconvergent flux recovery of the Rannacher-Turek nonconforming element

This work presents superconvergence estimates of the Rannacher-Turek ele...
10/23/2019 ∙ by Yuwen Li, et al. ∙ 0

• ### Generalized Prager-Synge Inequality and Equilibrated Error Estimators for Discontinuous Elements

The well-known Prager-Synge identity is valid in H^1(Ω) and serves as a ...
01/24/2020 ∙ by Cuiyu He, et al. ∙ 0

• ### Optimal maximum norm estimates for virtual element methods

The maximum norm error estimations for virtual element methods are studi...
05/25/2021 ∙ by Wen-Ming He, et al. ∙ 0

• ### Adaptive virtual elements based on hybridized, reliable, and efficient flux reconstructions

We present two a posteriori error estimators for the virtual element met...
07/08/2021 ∙ by F. Dassi, et al. ∙ 0

• ### Hanging nodes in the Context of Algebraic Stabilizations for Convection-Diffusion Equations

In this work, we propose an investigation of hanging nodes in the contex...
07/16/2020 ∙ by Abhinav Jha, et al. ∙ 0

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1. Introduction

In this paper, we consider the following 2D toy model on ,

 (1.1) {−∇⋅(β∇u)=f, in Ω,u=0, on ∂Ω.

To approximate (1.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 [25, 47, 2] with lots of mature software packages [3, 1] as it provides preferable performance in the aspects of the accuracy and robustness [48] over simplicial ones. To guide the AMR, one possible way is through the a posteriori error estimation [39] 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 [50, 46, 4, 38].

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.1) as example, given that the data , the flux is in , which has less continuity constraint than the aforementioned ones which are vertex-patch based with the recovered flux being -conforming. The -flux recovery shows more robustness than vertex-patch based ones (e.g., [15, 44, 34, 14]), some of them even have reliability constant being 1 [11, 45] and robust to the polynomial degrees [10, 30].

However, these -flux recovery techniques work mainly on conforming mesh. For nonconforming discretizations on nonmatching grids, some simple treatment of hanging nodes exists by recovering the flux on a conforming mother mesh (see e.g., [29]). To our best knowledge, there is no literature about the local -flux recovery on the multilevel irregular quadtree meshes. One major difficulty is that it is hard to construct a robust recovered flux locally to satisfy the -continuity constraint.

More recently, a new class of methods called the virtual element methods (VEM) were introduced in [6, 5, 12]

, 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, while the functions are represented by merely the degrees of freedom 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 -conforming finite dimensional approximation space on meshes with multilevels of irregularities. We shall follow this approach to perform flux recovery for conforming discretization of problem (1.1). Recently, arbitrary level of irregular quadtree meshes have been studied in [26, 42, 21]. Analyses of the residual-based error estimator on 1-irregular (balanced) quadtree mesh can be found in [19, 49, 33, 41]. In the virtual element context, vertex-patch recovery techniques are studied for linear elasticity in [24], and for diffusion problems in [32].

The major ingredient in our study is an -conforming virtual element space modified from the ones used in [12, 7] (Section 2.2). Afterwards, an

-conforming flux is recovered by a robust weighted averaging of the numerical flux, in which some special properties of the tensor-product type element 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 (Section 5).

## 2. Preliminaries

### 2.1. Discretization and notations

Given a shape-regular rectangular partition of , is assumed to be a piecewise, positive constant with respect to . The weak form of problem (1.1) is then discretized in a tensor-product finite element space as follows,

 (2.1) (β∇u\raisebox−0.5pt$T$,∇v\raisebox−0.5pt$T$)=(f,v\raisebox−0.5pt$T$),∀v\raisebox−0.5pt$T$∈Qk(T)∩H10(Ω),

in which the standard notation is opted. denotes the inner product on , and , with the subscript omitted when . The discretization space is

 Qk(T):={v∈H1(Ω):v|K∈Qk(K),∀K∈T}.

and on

 Qk(K):=Pk,k(K)={p(x)q(y),p∈Pk([a,b]),q∈Pk([c,d])}.

Henceforth, we shall simply denote when no ambiguity arises.

On , the sets of 4 vertices, as well as 4 edges of the same generation with , are denoted by and , respectively. The sets of nodes and edges in are denoted by and . A node is called a hanging node if it is on but is not counted as a vertex of , and we denote the set of hanging nodes as

 NH:={z∈N:∃K∈T,z∈∂K∖NK}

Otherwise the node is a regular node. If an edge contains at most hanging nodes, the partition , as well as the element these hanging nodes lie on, is called -irregular.

For each edge

, a unit normal vector

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 , then this element is denoted by , otherwise it is denoted by , i.e., is pointing from to . For any function or distribution well-defined on the two elements, the intersection of whose closures is an edge . We note that it is possible that or vice versa if and . Define on an edge , in which and are defined in the limiting sense for . If is a boundary edge, the function is extended by zero outside the domain to compute . Furthermore, the following notation denotes a weighted average of on edge ,

 {v}γe:=γv−+(1−γ)v+.

### 2.2. Virtual element spaces

In this subsection, the quadtree mesh of interest is embedded into a polygonal mesh . On any given quadrilateral element , for example we consider a , it has 4 degrees of freedom associated with 4 nodes . Its normal flux is well-defined on the 4 edges , 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 is linear function in -variable.

For the embedded element , which geometrically coincides with , 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 as . Using Figure 1 as example, it is possible to define a flux on with piecewise linear normal component on which now consists of three edges on .

Subsequently, shall be denoted by simply in the context of flux recovery, and the notion denotes an edge on the boundary of , which takes into account of the edges formed with one end point or both end points as the hanging nodes.

On , we consider the following Brezzi-Douglas-Marini-type virtual element modification inspired by the ones used in [12, 7]. The local space on a is defined as for

 (2.2) Vk(K):={ τ∈H(div;K)∩H(rot;K): ∇⋅τ∈Pk−1(K),∇×τ=0, τ⋅ne∈Pk(e),∀e⊂∂K}.

An -conforming global space for recovering the flux is then

 (2.3) Vk:={τ∈H(div):τ|K∈Vk(K), on K∈Tpoly}.

Next we turn to define the degrees of freedom (DoFs) of this space. An edge parametrized it by , where is the starting point of , and is the unit tangential vector of . The basis set for is chosen as the set of scaled monomials :

 (2.4) Pk(e):=span{1,s−mehe,(s−mehe)2,…,(s−mehe)k},

where representing the midpoint. Similar to the edge case, ’s basis set is chosen as follows:

 (2.5) Pk(K):=span{mα(x):=(x−xKhK)α,|α|≤k}

The degrees of freedom (DoFs) are then set as follows for a :

 (2.6) (e)k≥1 ∫e(τ⋅ne)mds,∀m∈Pk(e), on e⊂Epoly. (i)k≥2 ∫Kτ⋅∇mdx,∀m∈Pk−1(K)/R on K∈Tpoly.
###### Remark 2.1.

We note that in our construction, the degrees of freedom to determine the curl of a VEM function originally in [12] are replaced by a curl-free constraint thanks to the flexibility to virtual element. The unisolvency of the set of DoFs (2.6) including the curl-part can be found in [12], while for the modified space (2.2), a simplified argument is in the proof of Lemma 7.5.

## 3. Flux recovery

As the data , the true flux . Consequently, we shall seek a postprocessed flux in by specifying the DoFs in (2.6). Throughout this section, whenever considering an element , we treat it a polygon as .

### 3.1. Virtual element-based flux recovery

Consider the numerical flux on . We note that . The normal flux on each edge is in as and on vertical edges, and on horizontal edges. Therefore, the edge-based DoFs can be computed by simple averaging thanks to the matching polynomial degrees of the numerical flux to the functions in .

On each , define

 (3.1) {−β∇u\raisebox−0.5pt$T$}γee⋅ne:=(γe(−βK−∇u\raisebox−0.5pt$T$|K−)+(1−γe)(−βK+∇u\raisebox−0.5pt$T$|K+))⋅ne,

where

 (3.2) γe:=β1/2K+β1/2K++β1/2K−.

For the lowest order case , the normal component of the recovered flux is set as

 (3.3) σ\raisebox−0.5pt$T$⋅ne={−β∇u\raisebox−0.5pt$T$}γee⋅ne.

By (2.2), is a constant on , thus it can be determined by edge DoFs in (2.6):

 (3.4) |K|∇⋅σ\raisebox−0.5pt$T$=∫K∇⋅σ\raisebox−0.5pt$T$dx=∫∂Kσ\raisebox−0.5pt$T$⋅nds=∑e⊂∂K∫eσ\raisebox−0.5pt$T$⋅neds.

For , after the normal component (3.3) is set, furthermore on each , denote stands for the -projection to , and we let

 (3.5) ∇⋅σ\raisebox−0.5pt$T$=Πk−1f+cK.

Here is a constant added to ensure the compatibility of locally,

 (3.6) cK=1|K|(−∫KΠk−1fdx+∑e⊂∂K∫e{−β∇u\raisebox−0.5pt$T$}γee⋅neds),

Consequently, the set of DoFs can be set as:

 (3.7)

### 3.2. Locally projected flux

To the end of constructing a computable local error indicator, inspired by the novel idea of VEM framework, the recovered flux is projected to a space with a much simpler structure. A local oblique projection is defined as follows:

 (3.8) (Πτ,∇p)K=(τ,∇p)K,∀p∈Pk(K)/R.

Next we are gonna show that this projection operator can be straightforward computed for vector fields in .

#### 3.2.1. k=1

When , we can compute the right hand side of (3.8) as follows:

 (3.9) (τ,∇p)K=−(∇⋅τ,p)K+(τ⋅n,p)∂K.

By definition of the space (2.2) when , is a constant on and can be determined by edge DoFs in (2.6) similar to (3.4). Moreover, , thus the boundary term can be evaluated using DoFs in (2.6).

#### 3.2.2. k≥2

When , the right hand side of (3.8) can be evaluated following a similar procedure as (3.9), if we exploit the fact that , we have

 (3.10) (τ,∇p)K =−(∇⋅τ,Πk−1p)K+(τ⋅n,p)∂K =(τ,∇Πk−1p)K+(τ⋅n,p−Πk−1p)∂K,

which can be evaluated using both DoF sets and .

## 4. A posteriori error estimation

Given the recovered flux in Section 3, the recovery-based local error indicator and the element residual as follows:

 (4.1) ηf,K:=∥∥β−1/2(σ\raisebox−0.5pt$T$+β∇u\raisebox−0.5pt$T$)∥∥K, and ηr,K:=∥∥β−1/2(f−∇⋅σ\raisebox−0.5pt$T$)∥∥K,

then

 (4.2)

A computable is defined as:

 (4.3) ˆηf,K:=∥∥β−1/2KΠ(σ\raisebox−0.5pt$T$+βK∇u\raisebox−0.5pt$T$)∥∥K,

with the oblique projection defined in (3.8). The stabilization part is

 (4.4) ηs,K:=∣∣β−1/2K(I−Π)(σ\raisebox−0.5pt$T$+βK∇u\raisebox−0.5pt$T$)∣∣S,K.

Here is seminorm induced by the following stabilization

 (4.5) SK(v,w):=∑e⊂∂Khe(v⋅ne,w⋅ne)e+∑α∈Λ(v,∇mα)K(w,∇mα)K,

where is the index set for the monomial basis of with cardinality , i.e., the second term in (4.5) is dropped in the case.

The computable error estimator is then

 (4.6) ˆη2=⎧⎪⎨⎪⎩∑K∈T(ˆη2f,K+ˆη2s,K)when k=1,∑K∈T(ˆη2f,K+ˆη2s,K+η2r,K)when k≥2.

### 4.1. Efficiency

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

###### Theorem 4.1.

Let be the solution to problem (2.1), and be the error indicator in (4.6). On , can be locally bounded by the residual-based ones:

 (4.7) ˆη2f,K≲osc(f;K)2+η2R,K+η2J,K,

where

 osc(f;K) =β−1/2KhK∥∥f−Πk−1f∥∥K, ηR,K :=β−1/2KhK∥∥f+∇⋅(β∇u\raisebox−0.5pt$T$)∥∥K, and ηJ,K :=(∑e⊂∂KheβK+βKe∥∥[[β∇u\raisebox−0.5pt$T$⋅ne]]\raisebox−2.0pt$$∥∥2e)1/2. In the edge jump term, is the element on the opposite side of with respect to an edge . The constant depends on and the number of edges on . ###### Proof. Let on , then and we have  (4.8) ˆη2f,K =(Π(σ\raisebox−0.5ptT+βK∇u\raisebox−0.5ptT),∇p)K=(σ\raisebox−0.5ptT+βK∇u\raisebox−0.5ptT,∇p)K =−(∇⋅(σ\raisebox−0.5ptT+βK∇u\raisebox−0.5ptT),p)K+∑e⊂∂K∫e(σ\raisebox−0.5ptT+βK∇u\raisebox−0.5ptT)⋅nepds. By (3.3), without loss of generality we assume (the local orientation of agrees with the global one), and which is the element opposite to with respect to , and , we have on edge  (4.9) (σ\raisebox−0.5ptT+βK∇u\raisebox−0.5ptT)⋅ne =((1−γe)βK∇u\raisebox−0.5ptT|K−(1−γe)βKe∇u\raisebox−0.5ptT|Ke)⋅ne =β1/2Kβ1/2K+β1/2Ke[[β∇u\raisebox−0.5ptT⋅ne]]\raisebox−2.0pte. The boundary term in (4.8) can be then rewritten as  (4.10) ∫e(σ\raisebox−0.5ptT+βK∇u\raisebox−0.5ptT)⋅nepds = ≲ 1(βK+βKe)1/2h1/2e∥∥[[β∇u\raisebox−0.5ptT⋅ne]]\raisebox−2.0pt$$∥∥eβ1/2Kh−1/2e∥p∥e.

By a trace inequality on an edge of a polygon (Lemma 7.2), and the Poincaré inequality for , we have,

 h−1/2e∥p∥e≲h−1K∥p∥K+∥∇p∥K≲∥∇p∥K.

As a result,

 ∑e⊂∂K∫e(σ\raisebox−0.5pt$T$+βK∇u\raisebox−0.5pt$T$)⋅nepds≲ηJ,Kβ1/2K∥∇p∥e=ηJ,Kˆηf,K.

For the bulk term in (4.8), when , by (3.4), the representation in (4.10), and the Poincaré inequality for again with , we have

 −(∇⋅(σ\raisebox−0.5pt$T$+βK∇u\raisebox−0.5pt$T$),p)K≤∣∣∇⋅(σ\raisebox−0.5pt$T$+βK∇u\raisebox−0.5pt$T$)∣∣|K|1/2∥p∥K ≤ 1|K|1/2∣∣∣∫K∇⋅(σ\raisebox−0.5pt$T$+βK∇u\raisebox−0.5pt$T$)dx∣∣∣∥p∥K = 1|K|1/2∣∣ ∣∣∑e⊂∂K∫e(σ\raisebox−0.5pt$T$+βK∇u\raisebox−0.5pt$T$)⋅neds∣∣ ∣∣∥p∥K ≤ ⎛⎜⎝∑e⊂∂K1β1/2K+β1/2Ke∥∥[[β∇u\raisebox−0.5pt$T$⋅ne]]\raisebox−2.0pt$$∥∥eβ1/2Khe⎞⎟⎠∥∇p∥ ≲ ηJ,Kˆηf,K. When , by (3.5),  (4.11) −(∇⋅(σ\raisebox−0.5ptT+βK∇u\raisebox−0.5ptT),p)K=−(Πk−1f+cK+∇⋅(βK∇u\raisebox−0.5ptT),p)K ≤ (∥∥f−Πk−1f∥∥K+∥∥f+∇⋅(β∇u\raisebox−0.5ptT)∥∥K+|cK||K|1/2)∥p∥K. The first two terms can be handled by combining the weights and from . For , it can be estimated straightforwardly as follows  cK|K|1/2 =1|K|1/2(−∫K(Πk−1f−f)dx−∫K(f+∇⋅(β∇u\raisebox−0.5ptT))dx ≤∥∥f−Πk−1f∥∥K+∥∥f+∇⋅(β∇u\raisebox−0.5ptT)∥∥K +1|K|1/2∑e⊂∂K∫e(βK∇u\raisebox−0.5ptT−{β∇u\raisebox−0.5ptT}γee)⋅neds ≤∥∥f−Πk−1f∥∥K+∥∥f+∇⋅(β∇u\raisebox−0.5ptT)∥∥K +∑e⊂∂Kβ1/2Kβ1/2K+β1/2Ke∥∥[[β∇u\raisebox−0.5ptT⋅ne]]\raisebox−2.0pt$$∥∥e.

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

 −(∇⋅(σ\raisebox−0.5pt$T$+βK∇u\raisebox−0.5pt$T$),p)K≲(osc(f;K)+ηR,K+ηJ,K)β1/2K∥∇p∥

and the theorem follows. ∎

###### Theorem 4.2.

Under the same setting with Theorem 4.1, let as the estimator in (4.4), we have

 (4.12) ˆη2s,K≲osc(f;K)2+η2R,K+η2J,K,

The constant depends on and the number of edges on .

###### Proof.

This theorem follows directly the norm equivalence Lemma 7.5, thus

 ∣∣β−1/2K(I−Π)(σ\raisebox−0.5pt$T$+βK∇u\raisebox−0.5pt$T$)∣∣S,K≲∣∣β−1/2K(σ\raisebox−0.5pt$T$+βK∇u\raisebox−0.5pt$T$)∣∣S,K,

while evaluating the DoFs and using (3.3) and (3.7) reverts us back to the proof of Theorem 4.1. ∎

###### Theorem 4.3.

Under the same setting with Theorem 4.1, on any with defined as the collection of elements in which share at least 1 vertex with

 (4.13) ˆηK≲osc(f;K)+∥∥β1/2∇(u−u\raisebox−0.5pt$T$)∥∥ωK,

with a constant independent of , but dependent on and the maximum number of edges in .

###### 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 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 4.4 (T is l-irregular [19, 49, 33, 41, 39]).

Any given is always refined from a mesh with no hanging nodes by a quadsecting red-refinement. For any two neighboring elements in , the difference in their refinement levels is for a uniformly bounded constant , i.e., for any edge , it has at most hanging nodes.

By Assumption 4.4, we denote the father -irregular mesh of as . On , a subset of all nodes is denoted by , which includes the regular nodes on , as well as as the set of end points of edges with a hanging node as the midpoint. By [19, Theorem 2.1], there exists a set of bilinear nodal bases associated with , such that

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 4.5 (Quasi-monotonicity of β[40, 8, 15, 28]).

On , let be the bilinear nodal basis associated with , with . For every element , there exists a simply connected element path leading to , which is a Lipschitz domain containing the elements where the piecewise constant coefficient achieves the maximum (or minimum) on .

Denote

 (4.14) πzv=⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩∫ωz∩ωm(z)vϕz∫ωz∩ωm(z)ϕz if% z∈Ω,0 if z∈∂Ω.

We note that if is a constant on , . A quasi-interpolation can be defined as

 (4.15) Iv:=∑z∈N1(πzv)ϕz.
###### Lemma 4.6 (Estimates for πz and I).

Under Assumption 4.4 and 4.5, the following estimates hold for any

 (4.16)

and for

 (4.17) ∑K⊂ωzh−2z∥β1/2(v−πzv)ϕz∥2K≲∥∥β1/2∇v∥∥2ωz,

in which , and here denotes the union of elements in sharing at least a node (hanging or regular) with .

###### Proof.

The estimate for follows from [8, Lemma 2.8]. For , its error estimates and stability only rely on the partition of unity property of the nodal basis set (see e.g., [20, 43]), therefore the proof follows the same argument with the ones used on triangulations in [8, Lemma 2.8] or [15, Lemma 4.2]. ∎

Denotes the subset of nodes in (a) on the boundary as and (b) with high contrast coefficient on patch as . For the lowest order case, we need the following oscillation term for

 (4.18) osc(f;T)2:=∑z∈