In this paper, we consider the following 2D toy model 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  over simplicial ones. 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 [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., ). 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.
, 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 , and for diffusion problems in .
-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 (Section3). 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.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,
in which the standard notation is opted. denotes the inner product on , and , with the subscript omitted when . The discretization space is
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
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 vectoris 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 ,
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.
An -conforming global space for recovering the flux is then
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 :
where representing the midpoint. Similar to the edge case, ’s basis set is chosen as follows:
The degrees of freedom (DoFs) are then set as follows for a :
We note that in our construction, the degrees of freedom to determine the curl of a VEM function originally in  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 , 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
For the lowest order case , the normal component of the recovered flux is set as
For , after the normal component (3.3) is set, furthermore on each , denote stands for the -projection to , and we let
Here is a constant added to ensure the compatibility of locally,
Consequently, the set of DoFs can be set as:
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:
Next we are gonna show that this projection operator can be straightforward computed for vector fields in .
4. A posteriori error estimation
Given the recovered flux in Section 3, the recovery-based local error indicator and the element residual as follows:
A computable is defined as:
with the oblique projection defined in (3.8). The stabilization part is
Here is seminorm induced by the following stabilization
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
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.
Let on , then and we have
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
The boundary term in (4.8) can be then rewritten as
By a trace inequality on an edge of a polygon (Lemma 7.2), and the Poincaré inequality for , we have,
As a result,
When , by (3.5),
The first two terms can be handled by combining the weights and from . For , it can be estimated straightforwardly as follows
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
and the theorem follows. ∎
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
with a constant independent of , but dependent on and the maximum number of edges in .
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.
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.
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 .
We note that if is a constant on , . A quasi-interpolation can be defined as
Lemma 4.6 (Estimates for and ).
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