Gradient Recovery for the BEM‐based FEM and VEM

In this article, we propose new gradient recovery schemes for the BEM‐based Finite Element Method (BEM‐based FEM) and Virtual Element Method (VEM). Supporting general polytopal meshes, the BEM‐based FEM and VEM are highly flexible and efficient tools for the numerical solution of boundary value problems in two and three dimensions. We construct the recovered gradient from the gradient of the finite element approximation via local averaging. For the BEM‐based FEM, we show that, under certain requirements on the mesh, superconvergence of the recovered gradient is achieved, which means that it converges to the true gradient at a higher rate than the untreated gradient. Moreover, we propose a simple and very efficient a posteriori error estimator, which measures the difference between the unprocessed and recovered gradient as an error indicator. Since the BEM‐based FEM and VEM are specifically suited for adaptive refinement, the resulting adaptive algorithms perform very well in numerical examples.


Introduction
Gradient or stress recovery has a long tradition in the context of Finite Element Methods (FEM) [1][2][3][4]. In essence, gradient recovery is a post-processing applied to the gradient of the finite element solution and it has two major applications, namely the construction of superconvergent solutions [5] and a posterior error estimation [6]. In the past decades, it has been extensively studied and several different post-processing strategies have been proposed, for example the L 2 -recovery [7] or the popular patch recovery by Zienkiewicz and Zhu [8].
Yet, the application of gradient recovery is almost limited to standard FEM on triangular or quadrilateral discretisations. In [9], Guo, Xie and Zhao propose a Zienkiewicz-Zhu-type recovery scheme for Virtual Element Method (VEM), which is a new FEM-like method for the numerical solution of partial differential equations [10,11]. The VEM belongs to the family of Galerkin methods based on polytopal grids and, as such, features a great flexibility with handling complex geometries and allows for easy adaptive refinement and coarsening. Despite its newness, it has already been applied to a wide range of problems [12][13][14]. Besides the VEM, the hybridised discontinuous Galerkin method [15], the mimetic finite difference method [16] and the BEM-based Finite Element Method (BEM-based FEM) [17] are prominent examples for numerical methods on polytopal grids. The latter has been introduced in [18] and has been studied, in particular, for adaptive FEM strategies involving residual [19,20] and goal-oriented error estimators [21]. Other fields of application of the BEM-based FEM include, but are not limited to, convection dominated problems [22], anisotropic discretisations [23] and Nyström-based formulations [24].
In this work, we formulate gradient recovery schemes by averaging for the lowest order BEM-based FEM and VEM. In Section 3, we show that for the BEM-based FEM the centroids of regular k-gons are points of extraordinary accuracy, which we use to construct superconvergent recovered gradients in Section 4. Thereafter, we propose an a posteriori error estimator based on this recovery scheme in Section 5.

Preliminaries
Let Ω ⊂ R 2 be a bounded polygonal domain that admits a finite decomposition T h into open non-overlapping polygonal elements E with maximal diameter h > 0 such that Ω = E∈T h E. The boundary ∂E of each element E is assumed to be not self-intersecting. Moreover, we denote by x i , i = 1, . . . , N , the set of interior nodes of T h .
We denote by L 2 (Ω) the space of square-integrable functions and by H k (Ω) the Sobolev space of order k ∈ N with corresponding norms · L 2 (Ω) and · H k (Ω) . Furthermore, we define H k 0 (Ω) to be the closure of C ∞ 0 (Ω) in H k (Ω). For simplicity, we consider the Poisson equation with right-hand side f ∈ L 2 (Ω) and zero Dirichlet conditions, i.e., find In the following, we introduce the shape functions of lowest order used by both methods. Let P p (E) and P p (e) be the spaces of polynomials of degree p on the element E and edge e respectively. We define the local space of shape functions by Consequently, the global finite element space and degrees of freedoms are given by The two methods depart in the discretisation of (1). The BEM-based FEM uses the theory of boundary integral operators, whereas the VEM reduces the problem to the polynomial component We refer to [10,17] for details. Nonetheless, in both cases, we end up with a discrete formulation of the form: find u h ∈ V h such that By introducing the Lagrangian basis we can rewrite the discrete formulation (2) as a system of linear equations

Superconvergent Points
Originally, the term "gradient superconvergence" describes the phenomenon that for certain types of elements there exist points at which the gradient of the finite element solution converges to the true solution at a higher rate than that encountered globally. Such points of extraordinary accuracy are known as stress points and were first discovered by Barlow [1]. To the best of our knowledge, there are no results on superconvergent points for the FEM on polygonal elements. In the following, we extend the strategy of locating stress points by Strang and Fix [2] to the BEM-based FEM on regular polygonal elements. To this end, let Ω ⊂ R 2 be a regular convex or star polygon with k vertices. Note that the number of degrees of freedom equals k for convex and 2k for star k-gons. Given a quadratic harmonic polynomial m ∈ P 2 (Ω), we definem to be the BEM-interpolant in terms of the shape functions ϕ i . Then, we may characterise stress points as those points, wherem coincides with m up to a BEM error. We check if the centroid x c of Ω is superconvergent, i.e., if the interpolantm is almost exact at x c .

Superconvergent Gradient Recovery
The primary use of superconvergent points lies in superconvergent gradient recovery. To be more precise, we apply a postprocessing technique that incorporates superconvergent points, with the intention that the recovered gradient is superconvergent not only at certain points but throughout subdomains or even the whole domain.  In the following, we assume the mesh to be made of regular hexagons. We define the recovery operator G : Since this translates to local averaging, the recovery process is fairly simple and localised and practically does not increase numerical costs. Note that ∇u h is only given implicitly and therefore approximated by means of the BEM. Moreover, we observe superconvergence of the recovered gradient in numerical experiments. To demonstrate this, we solve the Laplace problem with Dirichlet conditions g(x, y) = exp(2π(x − 0.3)) cos(2π(y − 0.3)) on Ω = (−1, 1) 2 for a series of meshes made of regular hexagons and apply the averaging technique (3) on each level.
The error against the number of degrees of freedom for this example is depicted in Figure 1. We observe superconvergence, which is highlighted by an increase in the order of convergence, i.e., from O(N −1/2 ) for ∇u h to approximately O(N −3/4 ) for Gu h with N being the number of degrees freedom in the mesh.

Recovery-based error estimators
Another significant application of gradient recovery lies in a posteriori error estimation with the aim of adaptive mesh refinement techniques. In order to take general meshes into account, we modify our averaging scheme as follows Since u h is only given implicitly in the VEM, we apply the scheme to its projection Π ∇ 1 u h instead, compare [25].
www.gamm-proceedings.com Algorithm 1 We follow the basic concept of an adaptive FEM algorithm of the form 1. Solve: We solve on the current mesh level and compute ∇u h (BEM-based FEM) or ∇Π ∇ 1 u h (VEM). 2. Recovery: Then, we apply the post-processing by averaging (4) and obtain the recovered gradient Gu h (BEM-based FEM) or GΠ ∇ 1 u h (VEM). 3. Estimate: Subsequently, we compute the estimates Mark: Afterwards, we mark the elements for refinement. Here, we apply the Dörfler marking strategy [26].

5.
Refine: Finally, we refine the marked elements and start again. Here, we use the bisection algorithm introduced in [19].
We stop the algorithm if the maximum mesh level is reached or η is sufficiently small.
In the following, we test the performance of the recovery-based estimators for both methods. To this end, we consider the Laplace problem on the L-shaped domain Ω = (−1, 1) 2 \ [0, 1] 2 with Dirichlet conditions g(r, ϕ) = r 2/3 sin(2(ϕ − π/2)/3) in polar coordinates (r, ϕ). This example is a popular benchmark for adaptive algorithms, since the rate of convergence of the error in the energy norm is limited to O(N −1/3 ) for uniform refinement. As depicted in Figure 2, the recoverybased estimators recover the optimal convergence rate of O(N −1/2 ) and produce similar results compared to residual-based estimators [20,27]. In addition, we study the efficiency Ψ of the estimators, which is the quotient of estimated and true error. We see in Figure 3 that the recovery-based estimators operate nearly optimally and also outperforms the residual-based ones.
Overall, both estimators perform very well in this particular example. Due to the fact that the recovery-based scheme achieves slightly better results and estimates the error more accurately, we prefer it over the residual-based estimator.