Recovered Finite Element Methods

TL;DR

This paper introduces recovered finite element methods (R-FEM), a new family of Galerkin methods that combine features of classical and discontinuous Galerkin approaches, offering stability, optimal error bounds, and adaptability.

Contribution

The paper presents a novel family of R-FEM that unifies and extends classical finite element and discontinuous Galerkin methods with proven stability and error estimates.

Findings

R-FEM can produce stable conforming approximations.

Optimal a priori error bounds are established.

Numerical experiments demonstrate good approximation properties.

Abstract

We introduce a family of Galerkin finite element methods which are constructed via recovery operators over element-wise discontinuous approximation spaces. This new family, termed collectively as recovered finite element methods (R-FEM) has a number of attractive features over both classical finite element and discontinuous Galerkin approaches, most important of which is its potential to produce stable conforming approximations in a variety of settings. Moreover, for special choices of recovery operators, R-FEM produces the same approximate solution as the classical conforming finite element method, while, trivially, one can recast (primal formulation) discontinuous Galerkin methods. A priori error bounds are shown for linear second order boundary value problems, verifying the optimality of the proposed method. Residual-type a posteriori bounds are also derived, highlighting the potential of R-FEM in the context of adaptive computations. Numerical experiments highlight the good approximation properties of the method in practice. A discussion on the potential use of R-FEM in various settings is also included.

Figures

Equations

$${\int }_{\Omega }\nabla \mathcal{E}(u_h)\cdot \nabla \mathcal{E}(v_h)\mathrm{d}x+s(u_h,v_h)={\int }_{\Omega }f\mathcal{E}(v_h)\mathrm{d}x,\text{ for all }{v}_h\in V_h,$$

$$-\nabla \cdot A\nabla u=f\text{in }\Omega ,$$

$$V_h^r:=\{v\in L^2(\Omega ):v{|}_T\in {\mathcal{R}}_r(T),T\in \mathcal{T}\},$$

$$[v]:=q^{+}\mbox{\boldmath$\rm{n}$}^{+}+q^{-}\mbox{\boldmath$\rm{n}$}^{-};$$

$$\mathcal{E}:V_h^r\to H_0^1(\Omega )\cap V_h^s,$$

$$\mathcal{E}:V_h^r\to V_{nc},$$

$$B(w_h,v_h):=a(\mathcal{E}(w_h),\mathcal{E}(v_h))+s_h(w_h,v_h)=\ell (\mathcal{E}(v_h)),\text{for all }{v}_h\in \tilde{V},$$

$$a(w,v):=\sum _{T\in \mathcal{T}}{\int }_{T}A\nabla w\cdot \nabla v\mathrm{d}x,\text{for all }w,v\in \prod _{T\in \mathcal{T}}{H}^1(T),$$

$$\ell (v):={\int }_{\Omega }fv\mathrm{d}x,\text{for all }v\in \prod _{T\in \mathcal{T}}{H}^1(T),$$

$$\mathcal{E}_{s}(v)(\nu):=\bigg{\{}\begin{array}[]{cc}\displaystyle|\omega_{\nu}|^{-1}\sum_{T\in\omega_{\nu}}v|_{T}({\nu}),&\nu\in\Omega;\\ 0,&\nu\in\partial\Omega,\end{array}$$

$$\sum_{T\in\mathcal{T}}|v-\mathcal{E}_{s}(v)|_{\alpha,T}^{2}\leq C_{KP,|\alpha|}\|\mbox{\boldmath$\rm{h}$}^{1/2-\alpha}[v]\|_{\Gamma}^{2},$$

$$s_h(w_h,v_h):={\int }_{\Gamma }\sigma [w_h]\cdot [v_h]\mathrm{d}s,$$

$$\Vert \sqrt{A}\nabla \mathcal{E}(w){\Vert }_{\Omega }^2+\Vert \sqrt{\sigma }[w]{\Vert }_{\Gamma }^2=B(w,w)\forall w\in V_h^r.$$

$$\widetilde{s_{h}}(w_{h},v_{h}):=\int_{\Omega}\widetilde{\sigma}\big{(}w_{h}-\mathcal{E}(w_{h})\big{)}\big{(}v_{h}-\mathcal{E}(v_{h})\big{)}\,\mathrm{d}x,$$

$$\|\sqrt{A}\nabla\mathcal{E}(w)\|_{\Omega}^{2}+\|\sqrt{\widetilde{\sigma}}\big{(}w-\mathcal{E}(w)\big{)}\|_{\Omega}^{2}=B(w,w)\quad\forall\ w\in V_{h}^{r}.$$

$$c_{0}\|\mbox{\boldmath$\rm{h}$}^{\alpha-1/2}[v]\|_{\Gamma}^{2}\leq s_{h}(v,v)\leq C_{0}\|\mbox{\boldmath$\rm{h}$}^{\alpha-1/2}[v]\|_{\Gamma}^{2},\quad\text{ for all }v\in V_{h}^{r},$$

$$s_{h}(w,v)\leq\widetilde{C}_{0}\big{(}s_{h}(w,w)\big{)}^{1/2}\big{(}s_{h}(v,v)\big{)}^{1/2},\quad\text{ for all }w,v\in V_{h}^{r},$$

$$a(u-{\mathcal{E}}_s(u_h),{\mathcal{E}}_s(v_h))=s_h(u_h,v_h),\text{for all }{v}_h\in V_h^r.$$

$$a(u-{\mathcal{E}}_s(u_h),v_h)=0,\text{for all }{v}_h\in V_h^r\cap H_0^1(\Omega ),$$

$$a(u_h^{FEM},v_h)=\ell (v_h),\text{for all }{v}_h\in V_h^r\cap H_0^1(\Omega ).$$

$$s_{h}(w_{h},v_{h})=\int_{\Gamma}\Big{(}\sigma[w_{h}]\cdot[v_{h}]-\{A\nabla w_{h}\}\cdot[v_{h}]-\theta\{A\nabla v_{h}\}\cdot[w_{h}]\Big{)}\mathrm{d}s,$$

$$c_{0}\|\mbox{\boldmath$\rm{h}$}^{\alpha-1/2}[v]\|_{\Gamma}^{2}-\frac{1}{2}\|\sqrt{A}\nabla\mathcal{E}(v)\|_{\Omega}^{2}\leq s_{h}(v,v)\leq C_{0}\|\mbox{\boldmath$\rm{h}$}^{\alpha-1/2}[v]\|_{\Gamma}^{2}+\frac{1}{2}\|\sqrt{A}\nabla\mathcal{E}(v)\|_{\Omega}^{2},\text{ for all }v\in V_{h}^{r},$$

$$a(u-\mathcal{E}(u_h),\mathcal{E}(v_h))=s_h(u_h,v_h),\text{for all }{v}_h\in V_h^r.$$

$$\Vert \sqrt{A}\nabla (u-\mathcal{E}(u_h)){\Vert }_{\Omega }^2+s_h(u_h,u_h)=$$

$$=$$

$$\le$$

$$\displaystyle+\widetilde{C}_{0}\big{(}s_{h}(u_{h},u_{h})\big{)}^{1/2}\big{(}s_{h}(v_{h},v_{h})\big{)}^{1/2},$$

$$\|\sqrt{A}\nabla(u-\mathcal{E}(u_{h}))\|_{\Omega}^{2}+s_{h}(u_{h},u_{h})\leq\inf_{v_{h}\in V_{h}^{r}}\Big{(}\|\sqrt{A}\nabla(u-\mathcal{E}(v_{h}))\|_{\Omega}^{2}+\widetilde{C}_{0}^{2}s_{h}(v_{h},v_{h})\Big{)}.$$

$$\displaystyle\quad\inf_{v_{h}\in V_{h}^{r}}s_{h}(v_{h},v_{h})\leq s_{h}(\Pi_{r}u,\Pi_{r}u)\leq C_{0}\|\mbox{\boldmath$\rm{h}$}^{\alpha-1/2}[\Pi_{r}u]\|_{\Gamma}^{2}=C_{0}\|\mbox{\boldmath$\rm{h}$}^{\alpha-1/2}[u-\Pi_{r}u]\|_{\Gamma}^{2}$$

$$\displaystyle\leq C\sum_{T\in\mathcal{T}}\|\mbox{\boldmath$\rm{h}$}^{\alpha-1}(u-\Pi_{r}u)\|_{T}^{2}+\|\mbox{\boldmath$\rm{h}$}^{\alpha}\nabla(u-\Pi_{r}u)\|_{T}^{2}\leq C\sum_{T\in\mathcal{T}}h_{T}^{2(\alpha+s)}|u|_{s+1,T}^{2},$$

$$\quad\|\sqrt{A}\nabla(u-\mathcal{E}(u_{h}))\|_{\Omega}^{2}+c_{0}\|\mbox{\boldmath$\rm{h}$}^{\alpha-1/2}[u_{h}]\|_{\Gamma}^{2}\leq\inf_{v_{h}\in V_{h}^{r}}\|\sqrt{A}\nabla(u-\mathcal{E}(v_{h}))\|_{\Omega}^{2}+C\sum_{T\in\mathcal{T}}h_{T}^{2(\alpha+s)}|u|_{s+1,T}^{2},$$

$$\sigma =c_{\sigma }\mathcal{A}{\mathbf{h}}^{-1},$$

Full text

1 Introduction

Galerkin procedures are extremely popular in numerical approximation of solutions to initial and/or boundary value problems for partial differential equations (PDEs). The most used families of Galerkin procedures are the finite element (FEM) and, more recently, discontinuous Galerkin (dG) families of methods. Roughly speaking, FEM are attractive for their simplicity and robustness, especially in structural mechanics and heat flow simulations, owing to their variational interpretation and origins; dG methods, on the other hand, are popular in fluid flow and fast convection/transport simulations, due to their superior numerical stability properties, stemming from the ability to incorporate general numerical flux functions seamlessly.

FEM typically incorporate (approximate) continuity/conformity of the state variable(s) and/or of some moments directly into the finite element space in order to imitate the respective properties of the underlying continuous problem. As a result, FEM’s approximation capabilities have to be assessed for each choice of finite element spaces. At the other end of the spectrum, dG methods typically employ element-wise discontinuous approximation spaces whose approximation properties are clear; the continuity/conformity of the state variable(s) and/or of some moments is enforced weakly via numerical flux functions incorporated in the variational formulation of the dG method. Consequently, the approximation capabilities of dG schemes is only dependent on the ability to show Céa-type quasi-optimality results, i.e., relies on the structure of the respective weak formulation of the dG method. This enables the construction of stable dG methods in a number of settings, e.g., hyperbolic and/or convection-dominated problems [29, 26, 23, 15, 22], and locking-free approximations for elasticity [21, 31] (see also [6, 25] for similar ideas in the context of classical non-conforming methods).

In an effort to combine the “generic” approximation-space capabilities of dG methods while retaining the attractive conformity or near-conformity of the discretizations (enjoyed normally in FEM), we introduce a family of Galerkin finite element-type methods which are constructed via (conforming or classical non-conforming) recovery operators applied to element-wise discontinuous approximation spaces. This family will be termed collectively as recovered finite element methods (R-FEM) and has a number of attractive features over both FEM and dG approaches.

More specifically, R-FEM combines completely discontinuous local finite element-type spaces, resulting, nonetheless, to conforming FEM-like approximations, or classical non-conforming ones, e.g., of Crouzeix-Raviart type. To fix ideas, let $\mathcal{E}:V_{h}\to\tilde{V}_{h}\cap H^{1}_{0}(\Omega)$ an operator mapping a discontinuous piecewise polynomial space $V_{h}$ over a triangulation onto a space of continuous piecewise polynomial space $\tilde{V}_{h}\cap H^{1}_{0}(\Omega)$ over the same or a finer triangulation; such recovery operators $\mathcal{E}$ can be constructed locally, e.g., by (weighted) averaging of the nodal degrees of freedom [24, 13]. We can now consider the method: find $u_{h}\in V_{h}$, such that

[TABLE]

to solve the Poisson problem with homogeneous essential boundary conditions, $f\in H^{-1}(\Omega)$, for suitable stabilization $s(\cdot,\cdot):W_{h}\times W_{h}\to\mathbb{R}$, where $W_{h}\subset V_{h}$ such that $V_{h}=W_{h}\oplus(V_{h}\cap H^{1}_{0}(\Omega))$. We observe that, despite using element-wise discontinuous polynomial test space $V_{h}$, the method also produces simultaneously a conforming approximation $\mathcal{E}(u_{h})$. We note that the above method yields, in general, different numerical solutions to those one would get by postprocessing standard dG approximations via the recovery operator $\mathcal{E}$. On the other hand, as we shall see below, we can also retrieve classical conforming FEM approximations from R-FEM, by making specific choices of recovery operators $\mathcal{E}$. Therefore, in a sense R-FEM is both a generalisation and a variant of known finite element methods.

The above basic example is intended to highlight a number of attractive features for R-FEM: conformity is not hard-wired in the approximation spaces and there is considerable flexibility in the particular choice of: (a) the recovery operator $\mathcal{E}$; (b) the finite element space $\mathcal{E}(V_{h})$; and, (c) the stabilisation $s$ used. An important property of R-FEM is the extreme flexibility in the choice and nature of degrees of freedom to be recovered locally (e.g., nodal values in one element, normal fluxes or higher moments in another, etc.) while avoiding the often cumbersome proof of respective unisolvency of the global linear system. Moreover, a crucial practical attribute of R-FEM is that it can be implemented in a rather straightforward manner into existing conforming or non-conforming finite element implementations, as we shall discuss below.

The use of element-wise discontinuous polynomial spaces in standard dG methods is often found to increase the number of numerical degrees of freedom on a given mesh, without achieving better approximation per degree of freedom. Generally speaking, this is the case for low order approximation spaces. Recent work has shown that this is not necessarily the case for higher order polynomial degrees and/or $hp$-version dG methods [11, 10, 9], once appropriate choices of local basis are determined, e.g., using local basis of total polynomial degree on box-type or, generally, polygonal/polyhedral elements. Since the use of total degree bases on such meshes necessitates relaxation of conformity requirements for the approximation spaces, dG methods have been used in this context as the underlying discretization. It is evident that R-FEM can be naturally extended to accept such reduced approximation spaces, once a suitable choice of recovery is constructed; we refer to the forthcoming work [17] for the construction of R-FEM for general polygonal/polyhedral element shapes. Also, recalling that R-FEM is based on discontinuous approximation spaces also, R-FEM is expected to be able to produce stable conforming approximations in the context of convection-dominated problems. A numerical investigation in this direction is presented in the numerical experiments below.

The idea of variational methods involving some form of recovery, e.g., as (enhanced) gradient approximation, has appeared in various forms in the literature, especially in the context of dG methods, see, e.g., [18, 14, 8, 16]. All these approaches, however, recover the gradient of the numerical solution, through the solution of additional local (or global) weak problems. At the other end of the spectrum, the recent framework of virtual element methods (VEMs) [5] applied to meshes with polygonal/polyhedral element shapes, use non-polynomial basis functions to ensure conformity, which are not computed fully in practice. Indeed, the local VEM spaces are constructed so that they contain a, typically discontinuous, polynomial subspace ensuring optimal approximation; the approximate solution is then computed over this subspace and over the mesh skeleton.

As a first work introducing the recovered finite element method, we confine ourselves to its definition for linear elliptic problems on simplicial and/or box-type triangulations (Section 2). After discussing its relation on known, popular methods (Section 3), we prove a priori bounds (Section 4) showing the optimality of the proposed method with respect to meshsize $h$, as well as residual-type a posteriori error bounds (Section 5). We also investigate carefully the lowest order case, whereby the underlying discontinuous space $V_{h}$ consists of element-wise constant functions recovered into conforming linear elements on a triangular mesh; the resulting method is shown to converge optimally in the energy norm. We investigate a number of implementation and conditioning issues for the proposed method, assessing its competitiveness against FEM and dG in terms of complexity (Section 6). Numerical experiments highlighting the good approximation properties of the method in practice and comparisons against FEM and dG are performed in Section 7. Finally, we conclude by discussing a number of currently developed and future extensions of the R-FEM framework to complex problems.

Throughout this work we denote the standard Lebesgue spaces by $L^{p}(\omega)$, $1\leq p\leq\infty$, $\omega\subset\mathbb{R}^{d}$, $d=2,3$, with corresponding norms $\|\cdot\|_{L^{p}(\omega)}$; the norm of $L^{2}(\omega)$ will be denoted by $\|{\cdot}\|_{\omega}$ for brevity. Let also $H^{s}(\omega)$, be the Hilbertian Sobolev space of index $s\in\mathbb{R}$ of real-valued functions defined on $\omega\subset\mathbb{R}^{d}$, constructed via standard interpolation and/or duality procedures, along with the corresponding norm and seminorm $\|\cdot\|_{s,\omega}$ and $|\cdot|_{s,\omega}$, respectively. We also denote by $H^{1}_{0}(\omega)$ the space of functions in $H^{1}(\omega)$ with vanishing trace on $\partial\omega$.

For $\Omega$ a bounded open polygonal domain in $\mathbb{R}^{d}$, $d=2,3$, with $\partial\Omega$ denoting its boundary, we consider the elliptic problem

[TABLE]

where $f\in H^{-1}(\Omega)$, for some uniformly positive definite diffusion tensor $A\in[L^{\infty}(\Omega)]^{d\times d}$. For simplicity of the presentation only, we impose homogeneous essential boundary conditions $u=0$ on $\partial\Omega$, although this is by no means an essential restriction on what follows. Moreover, we shall also briefly consider a linear convection-diffusion problem in Section 7 to highlight the versatility of the proposed framework.

2 Recovered finite element method

Let $\mathcal{T}$ be a regular subdivision of $\Omega$ into disjoint simplicial or box-type (quadrilateral/hexahedral) elements $T\in\mathcal{T}$. We assume that the subdivision $\mathcal{T}$ is shape-regular (see, e.g., p.124 in [12]), that $\bar{\Omega}=\cup_{T\in\mathcal{T}}\bar{T}$ and that the elemental faces are straight line (for $d=2$) or planar (for $d=3$) segments; these will be, henceforth, referred to as facets. By $\Gamma$ we shall denote the union of all ($d-1$)-dimensional facets associated with the subdivision $\mathcal{T}$ including the boundary. Further, we set $\Gamma_{\text{\rm int}}:=\Gamma\backslash\partial\Omega$.

For a nonnegative integer $r$, we denote the set of all polynomials of total degree at most $r$ by $\mathcal{P}_{r}(T)$, while the set of all tensor-product polynomials on $T$ of degree at most $r$ in each variable is denoted by $\mathcal{Q}_{r}(T)$. For $r\geq 1$, we consider the finite element space

[TABLE]

where $\mathcal{R}_{r}(T)\in\{\mathcal{P}_{r}(T),\mathcal{Q}_{r}(T)\}$. We stress that, in this context, we can consider local bases in $\mathcal{P}_{r}(T)$ also for box-type elements; we shall return to this point below, cf., [11, 10, 9].

Further, let $T^{+}$, $T^{-}$ be two (generic) elements sharing a facet $e:=\partial T^{+}\cap\partial T^{-}\subset\Gamma_{\text{\rm int}}$ with respective outward normal unit vectors \mbox{\boldmath\rm{n}}^{+} and \mbox{\boldmath\rm{n}}^{-} on $e$. For a function $v:\Omega\to\mathbb{R}$ that may be discontinuous across $\Gamma_{\text{\rm int}}$, we set $v^{+}:=v|_{e\subset\partial T^{+}}$, $v^{-}:=v|_{e\subset\partial T^{-}}$, and we define the jump by

[TABLE]

if $e\in\partial T\cap\partial\Omega$, we set [v]:=v^{+}\mbox{\boldmath\rm{n}}. Also, we define $h_{T}:=\operatorname{diam}(T)$ and we collect them into the element-wise constant function ${\bf h}:\Omega\to\mathbb{R}$, with ${\bf h}|_{T}=h_{T}$, $T\in\mathcal{T}$, ${\bf h}|_{e}=(h_{T}+h_{T^{\prime}})/2$ for $e\subset\Gamma_{\text{\rm int}}$ and ${\bf h}|_{e}=h_{T}$ for $e\subset\partial T\cap\partial\Omega$. We assume that the families of meshes considered in this work are locally quasi-uniform.

For the definition of the proposed method, we require recovery operators of the form

[TABLE]

for some non-negative integer $r$, mapping element-wise discontinuous functions into functions in the solution space for the boundary value problem, for some $s,r\in\mathbb{N}\cup\{0\}$, or respective non-conforming recovery

[TABLE]

for some classical non-conforming finite element space $V_{nc}$, e.g., Crouzeix-Raviart elements. Evidently, there is considerable flexibility in the choice of such recovery operators; indeed, various choices of $\mathcal{E}$ may give rise to a different methods.

For recovery operator $\mathcal{E}:V_{h}^{r}\to\widetilde{V}$, with $\widetilde{V}\in\{H^{1}_{0}(\Omega)\cap V_{h}^{s},V_{nc}\}$, we consider the recovered finite element method reading: find $u_{h}\in V_{h}^{r}$ such that

[TABLE]

where

[TABLE]

and

[TABLE]

with $s_{h}(\cdot,\cdot):V_{h}^{r}\times V_{h}^{r}\to\mathbb{R}$ a symmetric bilinear form, henceforth referred to as the stabilisation; possible choices for $s_{h}$, resulting, in general, to different methods, will be given below. As we shall see below, specific choices of $\mathcal{E}$ and of $s_{h}$ can lead to the same recovered solutions produced by both classical FEM and dG methods.

Recovery operators of the form (3) have appeared in various settings in the theory of finite element methods, e.g., [13, 30, 27, 24, 7, 19]. They are typically used to recover a conforming function from a non-conforming one under minimal regularity requirements.

Classical examples include the construction of the so-called Clément or Scott-Zhang operators [13, 30], based on local averages of functions. Note that the classical Clément construction yields a quasi-interpolant, i.e., it does not necessarily preserve continuous functions.

Another popular example is the nodal averaging operator for which the following stability result was proven by Karakashian and Pascal in [24].

Lemma 2.1

Let $\mathcal{T}$ a locally quasi-uniform mesh. The operator $\mathcal{E}_{s}:V_{h}^{s}\to V_{h}^{s}\cap H^{1}_{0}(\Omega)$, defined on the conforming Lagrange nodes $\nu\in\mathcal{N}$, $\mathcal{N}$ denoting the set of all Lagrange nodes of $V_{h}^{s}$, by:

[TABLE]

with $\omega_{\nu}:=\bigcup_{T\in\mathcal{T}:\nu\in\bar{T}}T,$ the set of elements sharing the node $\nu\in\mathcal{N}$ and $|\omega_{\nu}|$ their number. Then, the following bound holds

[TABLE]

with $|\alpha|=0,1$, $C_{|\alpha|}\equiv C_{|\alpha|}(r)>0$ a constant independent of $\rm{h}$, $v$ and $\mathcal{T}$, but depending on the shape-regularity of $\mathcal{T}$ and on the polynomial degree $s$.

Proof.

See Karakashian and Pascal [24]. ∎

Note that, since $V_{h}^{r}\subset V_{h}^{s}$ for $r\leq s$, the recovery operator $\mathcal{E}_{s}$ from Lemma 2.1 is also an operator from $V_{h}^{r}$ into $V_{h}^{s}\cap H^{1}_{0}(\Omega)$ for $0\leq r\leq s$. The bound (7) shows, in particular, that \|\mbox{\boldmath\rm{h}}^{-1/2}[v]\|_{\Gamma}^{2} is a norm on the orthogonal complement $W_{h}^{s}$ of $V_{h}^{s}\cap H^{1}_{0}(\Omega)$ in $V_{h}^{s}$ with respect to the $a(\cdot,\cdot)$ inner product. This motivates the following choice for the stabilisation bilinear form:

[TABLE]

for some non-negative function $\sigma:\Gamma\to\mathbb{R}$, to be defined below, henceforth referred to as discontinuity-penalization parameter. In particular, since R-FEM (5) is defined on the discontinuous space $V_{h}^{r}$, one has to take into account the part of $V_{h}^{r}$ that is not in $V_{h}^{s}\cap H^{1}_{0}(\Omega)$: Lemma 2.1 suggests that, upon defining $s_{h}$ is in (8), the well-posedness of the method should be expected. Indeed, setting $w_{h}=v_{h}=w$ in (5) we immediately arrive at the coercivity identity

[TABLE]

With the help, now, of Lemma 2.1, we can show that the left-hand side of (9) is a norm in $V_{h}^{r}$. To this end, upon checking that $\|\sqrt{A}\nabla\mathcal{E}(w)\|_{\Omega}^{2}+\|\sqrt{\sigma}[w]\|_{\Gamma}^{2}=0$ implies $w=0$ (the other two norm properties being obvious). Indeed, we then have $\nabla\mathcal{E}(w)=0$ in $\Omega$ and $[w]=0$ on $\Gamma$. Since $\mathcal{E}(w)\in H^{1}_{0}(\Omega)$, this, in turn, implies $\mathcal{E}(w)=0$. Also, from Lemma 2.1, for $\alpha=0$, we have $w-\mathcal{E}(w)=0$, ensuring, thus, $w=0$.

Remark 2.2

The above argument may also motivate an alternative choice for the stabilisation $s$. In particular, choosing $s_{h}\equiv\widetilde{s_{h}}:V_{h}^{r}\times V_{h}^{r}\to\mathbb{R}$ with

[TABLE]

for some non-negative function $\widetilde{\sigma}:\Gamma\to\mathbb{R}$, as we we, then, have the coercivity identity

[TABLE]

To retain a level of generality in the R-FEM framework presented below, we make the following assumptions on the bilinear form $s_{h}$, rather than prescribe it specifically.

Assumption 2.3

We assume that the stabilisation $s_{h}\equiv s_{h,\alpha}$ satisfies the equivalence

[TABLE]

for some constant $\alpha\in\mathbb{R}$, and for $c_{0},C_{0}>0$, depending only on the local elemental polynomial degree and on the mesh regularity, topology and geometry, but not on the local mesh-size. (This means that any dependence on the local mesh-size $\rm{h}$ will be implicit in the definition of $s_{h}$ itself.)

Note that both choices of stabilisation (8) and (10) satisfy Assumption 2.3 with $\alpha=0$, while (8) also satisfies $c_{0}=C_{0}=1$.

Assumption 2.4

We assume that the stabilisation $s_{h}$ satisfies the bound

[TABLE]

for $\widetilde{C}_{0}>0$, depending only on the local elemental polynomial degree and on the mesh regularity, topology and geometry, but not on the local mesh-size.

It is evident that both choices of stabilisation (8) and (10) satisfy Assumption 2.4 with $\widetilde{C}_{0}=1$.

Remark 2.5

In the proof of a posteriori error estimates below, we shall make use of local versions of Assumptions 2.3 and 2.4, applied to subsets $\omega\subset\Omega$ of the computational domain. These can be easily shown to be valid by using the bilinearity of $s_{h}$ and selecting $supp(v)\subset\omega$, noting that (12) implies that $s_{h}(v,v)$ will vanish away from $\omega$.

3 Connections to known methods

Interestingly, employing the nodal averaging recovery operator $\mathcal{E}_{s}$ in the R-FEM formulation (5) on finite element spaces containing conforming subspaces, results to $\mathcal{E}_{s}(u_{h})$ being the approximate solution computed via classical conforming finite element method! Indeed, the boundary value problem in weak form reads: find $u\in H^{1}_{0}(\Omega)$ such that $a(u,v)=\ell(v)$ for all $v\in H^{1}_{0}(\Omega)$. Setting $v=\mathcal{E}_{s}(v_{h})$ for a $v_{h}\in V_{h}^{r}$, we deduce $a(u,\mathcal{E}_{s}(v_{h}))=\ell(\mathcal{E}_{s}(v_{h}))$, which, upon subtraction of (5), leads to

[TABLE]

Now, if $V_{h}^{r}$ is such that it contains the respective conforming subspace of the same order, setting $v_{h}\in V_{h}^{r}\cap H^{1}_{0}(\Omega)$, we have from (18)

[TABLE]

noting, furthermore, that $v_{h}=\mathcal{E}_{s}(v_{h})$ for $v_{h}\in V_{h}^{r}\cap H^{1}_{0}(\Omega)$. Therefore, (15) implies that $\mathcal{E}_{s}(u_{h})=u_{h}^{FEM}$ where, $u_{h}^{FEM}\in V_{h}^{r}\cap H^{1}_{0}(\Omega)$ is such that

[TABLE]

Thus, R-FEM for special choices of recovery operators $\mathcal{E}$ and of finite dimensional spaces $V_{h}^{r}$ can retrieve the classical FEM solutions (albeit in a wasteful fashion). This is the content of the following result, whose proof is essentially included in the above discussion.

Lemma 3.1

Let $V_{h}^{r}$ over a triangulation is such that it contains the conforming subspace $V_{h,c}^{r}:=V_{h}^{r}\cap H^{1}_{0}(\Omega)$ of the same polynomial degree. Assume also that the recovery operator used in (5) is such that it preserves conforming functions, i.e., we have $\mathcal{E}(v_{h})=v_{h}$, for all $v_{h}\in V_{h,c}^{r}$ and assume also that the stabilisation $s$ is such that $s(w_{h},v_{h})=0$. Then, $\mathcal{E}(u_{h})=u_{h}^{FEM}$, when $u_{h}\in V_{h}^{r}$ is computed via (5) and $u_{h}^{FEM}\in V_{h,c}^{r}$ is computed via (16).

We shall be explicit in the use any of the assumptions of Lemma 3.1 in the a priori and a posteriori analysis below.

At the other end of the spectrum, we can trivially obtain interior penalty discontinuous Galerkin methods by setting $\mathcal{E}$ to be the identity operator, i.e., no recovery and letting

[TABLE]

with $\{\cdot\}$ denoting the average operator on $\Gamma$, defined face-wise as $\{q_{h}\}|_{e}:=\frac{1}{2}(q_{h}|_{e\subset\partial T}+q_{h}|_{e\subset\partial T^{\prime}})$ for $e=\partial T\cap\partial T^{\prime}$ and $T,T^{\prime}$ neighbouring elements in $\mathcal{T}$; $\theta\in[-1,1]$ typically, and $\sigma$ denotes the usual discontinuity-penalization parameter. Note that Assumptions 2.3 and 2.4 are not satisfied verbatim in this case. However, provided the discontinuity penalization parameter $\sigma$ is chosen large enough, see, e.g., [1, 2, 22], it can be shown that

[TABLE]

and, respectively, for Assumption 2.4, allowing for the error analysis below to still hold. We refrained from using this, more general, version of Assumptions 2.3 and 2.4 in this work for simplicity of the presentation, as the benefit of generality does not seem to provide any significant new insight at this point.

Crucially, however, R-FEM offers significant flexibility in the choice of both the finite element spaces, of the recovery operators and of the stabilisation terms, thereby allowing also for new numerical methods. Indeed, many classical recovery operators, e.g., of Clément type do not satisfy the condition $\mathcal{E}(v_{h})=v_{h}$, for all $v_{h}\in V_{h,c}^{r}$, thereby giving rise to R-FEM with various properties even for the model elliptic problem considered above. A pertinent example will be given below where we consider recoveries from element-wise constant functions onto conforming linear elements or, in general, recovery between different finite element spaces. Furthermore, as we shall see below, for more general PDE problems, such as convection-diffusion equations the flexibility offered by the R-FEM construction allows for the construction stable methods in the convection-dominating regime. In this case, the use of $\mathcal{E}_{s}$ can result in non-standard/novel numerical methods.

4 A priori error analysis

The boundary value problem in weak form reads: find $u\in H^{1}_{0}(\Omega)$ such that $a(u,v)=\ell(v)$ for all $v\in H^{1}_{0}(\Omega)$. Setting $v=\mathcal{E}(v_{h})$ for a $v_{h}\in V_{h}^{r}$, we deduce $a(u,\mathcal{E}(v_{h}))=\ell(\mathcal{E}(v_{h}))$, which, upon subtraction of (5), leads to

[TABLE]

Using (18) and Assumption 2.4, we have, respectively,

[TABLE]

for all $v_{h}\in V_{h}^{r}$, with the last step following from the Cauchy-Schwarz inequality and Assumption 2.4. This implies the quasi-optimality bound

[TABLE]

For $\Pi_{r}:L^{2}(\Omega)\to V_{h}^{r}$ denoting the orthogonal $L^{2}$-projection operator onto the finite dimensional space $V_{h}^{r}$, we have the estimate

[TABLE]

for $0\leq s\leq\min\{r,l\}$, when $u\in\prod_{T\in\mathcal{T}}H^{l+1}(T)\cap H^{1}_{0}(\Omega)$, for $l\geq 0$, using the standard trace estimate and the best approximation properties of $\Pi_{r}$.

Now, using Assumption 2.3 and (21) on (20) already proves the following result.

Theorem 4.1

Let $u\in\prod_{T\in\mathcal{T}}H^{l+1}(T)\cap H^{1}_{0}(\Omega)$, for $l\geq 0$ be the solution to (1) and $u_{h}\in V_{h}^{r}$ its R-FEM approximation with the stabilisation term satisfying Assumptions 2.3 and 2.4. Then, we have the a priori bound

[TABLE]

for all $0\leq s\leq\min\{r,l\}$, with $C$ a positive constant, independent of $u,$ $u_{h}$, $\rm{h}$ and of $\mathcal{E}$.

Assuming now that the recovery operator $\mathcal{E}$ in the definition of R-FEM allows for optimal approximation of $u$ by $\mathcal{E}(v_{h})$ for some $v_{h}\in V_{h}^{r}$, one can recover optimal a priori error bounds. Below we highlight some interesting cases of $\mathcal{E}$ and polynomial order $r$.

4.1 Case $r=s\geq 1$.

We begin by considering the case of recovery into a conforming finite element space of the same polynomial degree. As discussed in Lemma 3.1, for certain choices of $\mathcal{E}$ this case results to $\mathcal{E}(u_{h})$ being the conforming FEM solution. Nonetheless, as we shall see below, this is not necessarily the case for more general elliptic operators involving lower order terms, (see Section 7.4 below,) or, indeed when the recovery takes place on a different triangulation than the underlying mesh $V_{h}^{r}$ is defined on; see Section 8.1 below.

When $r=s$, the discontinuous space $V_{h}^{r}$ contains sufficient approximation power to ensure $r$-th order convergence in the energy norm.

Corollary 4.2

Let the discontinuity-penalization parameter $\sigma$ be given by

[TABLE]

for any $c_{\sigma}>0$ and let the recovery operator $\mathcal{E}$ be such that $V_{h}^{r}\cap H^{1}_{0}(\Omega)\subset\mathcal{E}(V_{h}^{r})$. Further, assuming that $u|_{T}\in H^{k+1}(T)$, $T\in\mathcal{T}$, $k\geq 1$, we have the a priori bound

[TABLE]

with $0\leq q\leq\min\{k,r\}$, for a $C$ positive constant, independent of the meshsize and of $u$.

∎

4.2 Case $0\leq r<s$.

One may wish to recover from an element-wise discontinuous polynomial space of order $r$, into an $H^{1}$-conforming polynomial space of strictly higher order $s$. In this setting, a recovered element method of optimal order at least with respect to $r$ is expected. The analysis presented in Section 4.1 above, with the choice (22) for the discontinuity-penalization parameter, immediately gives the error bound (23) with $0\leq q\leq\min\{k,r\}<s$. It is clear that the above bound is of lower order than the respective finite element method involving polynomials of order $s$. This is, perhaps, not surprising as the approximation power of $V_{h}^{r}$ allows, in general, for optimal rates with respect to $r$ only.

However, in practice there are some advantages in recovering from lower into higher polynomial order spaces. In particular, as we shall see below when using $r=s-1$, for $s=1,2,\dots$, with appropriately chosen power $\alpha$ in the stabilisation $s_{h}$, one can recover $s$-th order convergent method on simplicial meshes! To highlight this, somewhat surprising feature, we focus on the case $0=r<s=1$, i.e., recovering element-wise constants $u_{h}\in V_{h}^{0}$ into conforming element-wise linear elements $\mathcal{E}(u_{h})\in V_{h}^{1}\cap H^{1}_{0}(\Omega)$, where first order convergence in the energy norm and second order convergence in the $L^{2}$-norm is observed for a suitable choice of the discontinuity-penalization parameter on shape-regular meshes; see Section 7 for a numerical illustration.

To show this, setting

[TABLE]

for any $c_{\sigma}>0$, (and, consequently we shall have $\alpha=1$ in Assumption 2.3,) we deduce from Theorem 4.1 that

[TABLE]

If one can show, e.g., that $\mathcal{E}(V_{h}^{0})=V_{h}^{1}\cap H^{1}_{0}(\Omega)$, (at least for some appropriately designed meshes,) then the first term on the right-hand side of (25) becomes

[TABLE]

from standard Bramble-Hilbert type approximation results [12], i.e., optimal linear convergence is retrieved. We shall now complete the proof of (26) for the case of a two-dimensional simplicial mesh.

Lemma 4.3

For $d=2$, let $\mathcal{T}$ be a regular simplicial mesh. Then, we have $\mathcal{E}(V_{h}^{0})=V_{h}^{1}\cap H^{1}_{0}(\Omega)$.

Proof.

Recalling that $\mathcal{N}$ denotes the set of all Lagrange nodes $\nu\in\Omega$, we set $\mathcal{N}_{0}$ to be the set of respective Lagrange nodes situated on $\partial\Omega$. (Since we are concerned with linear elements on a simplicial mesh, $\mathcal{N}$ is equal to the set of internal nodes in the mesh.) To prove the result, it is sufficient to prove that each Lagrange nodal basis $\phi_{\nu}$ can be constructed as a recovery of an element-wise constant function $\psi_{\nu}\in V_{h}^{0}$, i.e., $\phi_{\nu}=\mathcal{E}(\psi_{\nu})$.

To construct $\psi_{\nu}\in V_{h}^{0}$, we split $\mathcal{N}$ into a union of disjoint subsets as follows. Let $\mathcal{N}_{1}\subset\mathcal{N}$ denote the set of nodes $\nu_{1}$ for which there exists a triangle $T\in\mathcal{T}$ having $\nu_{1}$ as a node and the two remaining nodes in $\mathcal{N}_{0}$. For $i=2,\dots,r$, for some $r$, let $\mathcal{N}_{i}\subset\mathcal{N}\backslash\cup_{j=1}^{i-1}\mathcal{N}_{j}$ denote the set of nodes $\nu_{i}$ for which there exists a triangle $T\in\mathcal{T}$ having $\nu_{i}$ as a node and the two remaining nodes in $\mathcal{N}_{i-1}$. The existence of an $r\in\mathbb{N}$ such that $\mathcal{N}=\cup_{j=1}^{r}\mathcal{N}_{j}$ follows from the assumption that $\mathcal{T}$ is a regular simplicial mesh.

Having constructed a pairwise disjoint subdivision $\{\mathcal{N}_{j}\}_{j=1}^{r}$ of $\mathcal{N}$, we can now give an algorithm for the construction of $\psi_{\nu}$. Let $r_{\nu}\in\mathbb{N}$ be such that $\nu\in\mathcal{N}_{r_{\nu}}$ and set $\psi_{\nu}=1=(-1)^{0}$ on the element $T_{r_{\nu}}\in\mathcal{T}$ with vertex $\nu$ and the two other vertices $\nu_{1},\nu_{2}$ in $\mathcal{N}_{r_{\nu}-1}$. For each of the $\nu_{i},$ there exists a triangle $T_{r_{\nu}-1}^{i}\in\mathcal{T}$, $i=1,2$, with one vertex being $\nu_{i}$ and the other two vertices being in $\mathcal{N}_{r_{\nu}-2}$. We now have 2 possibilities: 1) if $T_{r_{\nu}-1}^{1}$ and $T_{r_{\nu}-1}^{2}$ have a common vertex, say $\nu_{3}$, then there exists a triangle $T_{r_{\nu}-1}\in\mathcal{T}$ having vertices $\nu_{1},\nu_{2}$ and $\nu_{3}$ on which we set $\psi_{\nu}=-1$ (notice that $T_{r_{\nu}-1}$ is necessarily different to the $T_{r_{\nu}-1}^{i}$’s); if $T_{r_{\nu}-1}^{1}$ and $T_{r_{\nu}-1}^{2}$ have no common vertex, we set $\psi_{\nu}=-1=(-1)^{1}$ on both $T_{r_{\nu}-1}^{1},T_{r_{\nu}-1}^{2}$. We continue the above algorithm by setting $\psi_{\nu}=(-1)^{j}$ on $T_{r_{\nu}-j}^{i}$ for an index set $\mathcal{I}_{j}\ni i$, until we reach the boundary nodes in $\mathcal{N}_{0}$ for which the recovery imposes the homogeneous boundary conditions strongly. Noticing that on each node $\nu_{i}$ there exist exactly two elements with non-zero values one $+1$ and the other $-1$, we conclude that $\phi_{\nu}=c\mathcal{E}(\psi_{\nu})$ for some $c>0$, which already proves the result. ∎

Thus, we have proved the following a priori error bound.

Corollary 4.4

For $d=2$, let $\mathcal{T}$ be a regular and shape-regular simplicial mesh. Further, assuming that $u|_{T}\in H^{2}(T)$, $T\in\mathcal{T}$, and that $\alpha=1$ in Assumption 2.3, we have that (23) holds with $q=1$.

So, we conclude that it is possible to have an optimally convergent finite element method for second order elliptic problems based on element-wise constant finite element spaces. This remarkable order increasing property of the method when recovering from element-wise constants to conforming linear elements appears to hold for greater values of $r$. In particular, we have observed numerically that R-FEM with $r=s-1$, for $s=1,2,3$ converges with order $s$ in the $H^{1}$-norm and with order $s+1$ in the $L^{2}$-norm when the stabilisation coefficient $\sigma$ in (8) is chosen as

[TABLE]

for $c_{\sigma}>0$ constant, i.e., $\alpha=s$ in (7). We refer to the numerical experiments in Section 7.1 for a numerical illustration. The proof of this interesting property will be discussed elsewhere. We conjecture that this property is due to dimensional considerations since \dim V_{h}^{r}>\dim\big{(}V_{h}^{r+1}\cap H^{1}_{0}(\Omega)\big{)} for $r=1,2,3$. Nonetheless, it is not clear at this point why the scaling in (27) appears to be necessary for observing optimal convergence.

Interestingly, such order-increasing behaviour is not observed for recoveries of the type $r=s-2$, $s=2,3,\dots$. Nonetheless, as we shall see in the numerical experiments below, it can be still beneficial for the constant of the convergence rate to recover lower order discontinuous spaces into higher order conforming ones for ‘compatible’ choices of $\alpha$.

5 An a posteriori error bound

To highlight further the flexibility offered by the proposed R-FEM framework, we also prove a basic residual-type a posteriori error bound for R-FEM with conforming recoveries.

Theorem 5.1

Let $u$ be the solution of (1) and $u_{h},\mathcal{E}(u_{h})$ be the R-FEM solution defined through (5) with $\mathcal{E}:V_{h}^{r}\to V_{h}^{s}\cap H^{1}_{0}(\Omega)$, $r\in\mathbb{N}_{0}$ and $s\in\mathbb{N}$. Assume also that the recovery operator $\mathcal{E}$ satisfies (7) and that the stabilisation $s_{h}$ can be decomposed into local contributions $s_{h,T}$, so that

[TABLE]

(cf. Remark 2.5). Then, we have the bound

[TABLE]

where

[TABLE]

and $\eta_{A,T}:=\|(A-\Pi A)\nabla\mathcal{E}(u_{h})\|_{T}$, for $C>0$ an estimable-from-above constant depending on $C_{0},\widetilde{C}_{0}$ and on the shape-regularity of the mesh; $\Pi\equiv\Pi_{t}:L^{2}(\Omega)\to V_{h}^{t}$ denotes the orthogonal $L^{2}$-projection onto $V_{h}^{t}$ for some $t\in\mathbb{N}\cup\{0\}$; when $\Pi$ is applied to tensors it will be understood as acting component-wise.

Proof.

We have, respectively,

[TABLE]

from (18) for $\chi\in V_{h}^{r}$ to be chosen precisely below, setting $v:=u-\mathcal{E}(u_{h})-\mathcal{E}(\chi)$ for brevity. Integration by parts and re-ordering of the terms yields

[TABLE]

Application of Cauchy-Schwarz inequality results in

[TABLE]

Selecting $\chi=\Pi_{r}(u-\mathcal{E}(u_{h}))$, (with $\Pi_{r}:L^{2}(\Omega)\to V_{h}^{r}$ the orthogonal $L^{2}$-projection operator onto $V_{h}^{r}$,) we have

[TABLE]

Now, we have the bound

[TABLE]

from the best approximation properties of the $L^{2}$-projection and its stability in the local $H^{1}$-seminorm. Completely analogously, we can show also that

[TABLE]

and the bound

[TABLE]

Next, using a similar line of argument, we also have

[TABLE]

Moreover, Assumptions 2.4, 2.3 and Remark 2.5 imply

[TABLE]

and we also have

[TABLE]

Combining the above estimates and applying them to estimate the respective terms on the right-hand side of (30), we arrive at the a posteriori bound. ∎

Trivially, Assumption 2.3 implies also the bound

[TABLE]

Remark 5.2

The estimators $\eta_{A,T}$ can be thought as data oscillation [28]. Indeed, setting $v_{h}=u_{h}$ in (5), using the Poincaré-Friedrichs inequality $\|\mathcal{E}(u_{h})\|_{\Omega}\leq C\|\nabla\mathcal{E}(u_{h})\|_{\Omega}$ on the right-hand side of (5), along with the uniform ellipticity of $A$, we can arrive at the stability estimate

[TABLE]

which, in turn, implies

[TABLE]

Remark 5.3

We remark that the hypothesis that $\mathcal{E}$ should satisfy (7) of the previous theorem does not imply that the result is applicable only when the averaging operator $\mathcal{E}_{s}$ is used in the R-FEM formulation. It merely states that for any recovery for which (7) holds, the above a posteriori bound also holds.

Using standard tools, a lower bound for the error $\|\sqrt{A}\nabla(u-\mathcal{E}(u_{h}))\|_{\Omega}^{2}$ can also be given.

Theorem 5.4

With the hypotheses of Theorem 5.1, for every $T\in\mathcal{T}$, we have the bound

[TABLE]

where $\omega_{T}:=\cup_{T^{\prime}\in\mathcal{T}:\partial T^{\prime}\cap\partial T\neq\emptyset}T^{\prime}$.

6 Implementation and Complexity

We shall now describe the implementation of the R-FEM and some of the properties of the resultant algebraic system.

Let $N_{C}:=\text{dim}\!\left({V_{h}^{s}\cap H^{1}_{0}(\Omega)}\right)$ and $N_{D}:=\text{dim}\!\left({V_{h}^{r}}\right)$, the dimensions of the conforming and discontinuous spaces, respectively. The operator $\mathcal{E}:V_{h}^{r}\to V_{h}^{s}\cap H^{1}_{0}(\Omega)$ can be represented algebraically as a matrix $\boldsymbol{E}\in\mathbb{R}^{N_{C}\times N_{D}}$ such that for $\boldsymbol{u}\in\mathbb{R}^{N_{D}}$, $\boldsymbol{E}\boldsymbol{u}\in\mathbb{R}^{N_{C}}$. Denote also by $\boldsymbol{K}_{\rm FEM}\in\mathbb{R}^{N_{C}\times N_{C}}$ the stiffness matrix associated to the classical conforming finite element method (15). Then the algebraic form of (5) is given by seeking $\boldsymbol{u}\in\mathbb{R}^{N_{D}}$ such that

[TABLE]

for $b\in\mathbb{R}^{N_{D}}$ given by $b_{i}=\int_{\Omega}f\psi_{i}\mathrm{d}x$, $i=1,\dots,N_{D}$ and $\psi_{i}\in V_{h}^{r}$ being a basis of $V_{h}^{r}$ and $\boldsymbol{S}$ being the algebraic representation of the stabilisation $s_{h}$. For comparison’s sake, we also introduce the dG stiffness matrix $\boldsymbol{K}_{IP}\in\mathbb{R}^{N_{D}\times N_{D}}$ associated to the interior penalty dG method with solution $u_{h}^{IP}\in V_{h}^{r}$.

Theorem 6.1** **(Condition number estimate)

For $d=2$, let $\mathcal{T}$ be a locally quasiuniform triangulation, and assume that the basis functions of $V_{h}^{r}$ used in the computation have local support. Assume also that the recovery operator $\mathcal{E}$ satisfies (7) and the stabilisation is such that Assumptions 2.3 and 2.4 hold for $0\leq\alpha\leq 1$. Then, the $\ell_{2}$-condition number $\kappa(\mathfrak{A})$ of the R-FEM stiffness matrix $\mathfrak{A}$ satisfies

[TABLE]

with \underline{\mbox{\boldmath\rm{h}}}:=\min_{\Omega}\mbox{\boldmath\rm{h}}, for $C$ a positive constant independent of the dimension of $\mathfrak{A}$. Also, under the same assumptions, for $d=3$, we have $\kappa(\mathfrak{A})\leq CN_{D}.$

Proof.

Without loss of generality, we assume \mbox{\boldmath\rm{h}}\leq 1. It is sufficient to find $\lambda,\Lambda>0$ such that $\lambda\boldsymbol{v}^{\boldsymbol{\intercal}}\boldsymbol{v}\leq\boldsymbol{v}^{\boldsymbol{\intercal}}\mathfrak{A}\boldsymbol{v}\leq\Lambda\boldsymbol{v}^{\boldsymbol{\intercal}}\boldsymbol{v}$ for all $\boldsymbol{v}\in\mathbb{R}^{N_{D}}$; then $\kappa(\mathfrak{A})\leq\Lambda/\lambda$. From (7) and Assumption 2.3, we have, respectively,

[TABLE]

since the basis functions used have local support.

Now, working analogously as in the proof of Theorem 4.1 from [3], for $p\in(2,\infty]$, we have

[TABLE]

using standard inverse estimates, Hölder’s inequality, (7), Sobolev Imbedding Theorem, and setting \bar{\mbox{\boldmath\rm{h}}}:=\max_{\Omega}\mbox{\boldmath\rm{h}}, respectively. Selecting now p=\max\{2,|\log(\underline{\mbox{\boldmath\rm{h}}}^{2}N_{D})|\}, the result follows.

The result for $d=3$ is completely analogous, but it is omitted for brevity, as it essentially follows the steps from [3]. ∎

The above result shows that R-FEM enjoys similar conditioning to standard FEM and dG methods for $0\leq\alpha\leq 1$.

For the case $r=s-1$, for $s=1,2,3$ discussed above, an inspection of the proof of Theorem 6.1 shows that the conditioning of $\mathfrak{A}$ deteriorates when $\alpha>1$. This is an interesting challenge for R-FEM implementation, which will be investigated further elsewhere. Crucially, however, $\mathfrak{A}$ is the respective stiffness matrix of the linear system solving for $u_{h}$ and not of the recovered solution $\mathcal{E}(u_{h})$, for which the respective linear system is expected to be better conditioned in practice. This could be a starting point for the development of preconditioning strategies for R-FEM.

6.1 Matrix structure

Due to the compact support of the finite element basis functions, linear systems resulting from the finite element discretisation of PDEs are typically sparse. Sparsity of the stiffness matrix is traditionally an indicator of complexity of the solution of the respective linear system. To this end, we begin by studing the sparsity of the R-FEM stiffness matrix and compare with FEM and dG counterparts.

The recovery procedure can be implemented through local transition matrices. These are defined by locally rewriting the basis functions of the underlying nonconforming space in terms of the space used for reconstruction.

Example 6.2** **($\mathcal{E}:V_{h}^{0}\to V_{h}^{1}\cap H^{1}_{0}(\Omega)$)

Let $M_{0}=\dim{\mathbb{P}^{0}(K)}=1$ and $M_{1}=\dim{\mathbb{P}^{1}(K)}=3$. The local transistion matrix is a mapping which takes the form $\boldsymbol{T}_{K}\in\mathbb{R}^{M_{0}\times M_{1}}$ where

[TABLE]

The global recovery matrix $\boldsymbol{E}$ can then be constructed by appropriately weighting and summing over all elements.

Example 6.3** **($\mathcal{E}:V_{h}^{1}\to V_{h}^{2}\cap H^{1}_{0}(\Omega)$)

For $M_{0}=\dim{\mathbb{P}^{1}(K)}=3$ and $M_{1}=\dim{\mathbb{P}^{2}(K)}=6$. The local transition matrix $\boldsymbol{T}_{K}\in\mathbb{R}^{M_{0}\times M_{1}}$ is given by

[TABLE]

When considering the resulting full linear systems, we observe that the R-FEM yield system matrices that have the same dimension but are less sparse than the respective dG counterparts on the same space $\tilde{V}$. Indeed, the R-FEM stiffness matrix appears to have the sparsity pattern of a (conforming or non-conforming) finite element method applied to $V_{h}^{s}$ but the dimension of $V_{h}^{r}$, which may in general be larger than the one of $\tilde{V}$ on standard simplicial meshes. Moreover, the higher the regularity of the recovered space $\tilde{V}$, the larger the bandwidth of the resultant matrix. We refer to Figure 3 for a visualisation of the sparsity structure of the recovery matrix $\boldsymbol{E}$ and to Figure 4 for a comparison of the sparsity patterns for $\mathfrak{A}$, $K_{\rm FEM}$ and $K_{\rm IP}$ over the same mesh.

7 Numerical Benchmarking

We now illustrate the performance of the R-FEM through a series of numerical experiments. The implementation of R-FEM used is available at [20].

7.1 Test 1 – Asymptotic behaviour approximating a smooth solution

As a first test, we consider $A=I$ and the domain $\Omega=[0,1]^{2}$. We fix $f$ such that the exact solution is given by

[TABLE]

and approximate $\Omega$ through a uniformly generated, criss-cross triangular type mesh to test the asymptotic behaviour of the numerical approximation. The results are summarised in Figure 5 (a) – (c), and confirm the theoretical findings in Section 4. More specifically, upon selecting $\sigma$ as in (27), we witness $s=r+1$ order of convergence in the energy norm error $u-\mathcal{E}(u_{h})$, for $r=0,1,2$ and respective $s+1=r+2$ order of convergence in the $L^{2}$-norm.

As a further test, we consider the case $r=1$, $s=3$ and $\alpha=3$, i.e., when recovering discontinuous linear elements into conforming cubic elements. In this case, only first order convergence rate in the $H^{1}$-norm and second order convergence in the $L^{2}$-norm are observed, i.e., there is no improvement in the convergence rate. However, this recovery can be still beneficial for the constant of the convergence rate: in Figure 5(d) we observe a reduction in the absolute energy norm error against numerical degrees of freedom along with a comparison to the (classical) conforming linear FEM.

Next, we compare the lowest order R-FEM method with $\mathcal{E}:V_{h}^{0}\to V_{h}^{1}\cap H^{1}_{0}(\Omega)$ and (8) with (24) against the dG interior penalty method for the same problem and perform a series of asymptotic benchmarks varying the penalty parameter $c_{\sigma}$. Recall that for R-FEM $c_{\sigma}>0$ can be chosen arbitrarily, as opposed to interior penalty dG methods where $\sigma$ should be given by (22) with $c_{\sigma}$ sufficiently large so as to ensure coercivity of the bilinear form. The results are given in Figure 6. We fix the dG penalty parameter to be \sigma_{\rm IP}=10\mbox{\boldmath\rm{h}}^{-1} and vary the R-FEM penalty parameter, $c_{\sigma}$. Note that R-FEM behaves comparably to the interior penalty method and there is a region of values of $c_{\sigma}$ where $\|u-\mathcal{E}(u_{h})\|_{L^{2}(\Omega)}<\|u-u_{h}^{IP}\|_{L^{2}(\Omega)}$. For small values of $c_{\sigma}$ we see $\|\nabla(u-\mathcal{E}(u_{h}))\|_{\Omega}\approx\|u-u_{h}^{IP}\|_{\rm dG}$ with $\|\cdot\|_{\rm dG}$ denoting energy-like dG-norm defined as

[TABLE]

7.2 Test 2 – Asymptotic behaviour approximating a singular solution

We choose $f$ so that

[TABLE]

where $r^{2}=x^{2}+y^{2}$ and $\theta=\arctan(y/x)$ are the associated polar coordinates, over the domain $\Omega=(-1,1)^{2}\backslash[0,1)\times(-1,0]$. We triangulate the domain with a criss-cross triangular mesh. Note the solution is not $H^{2}(\Omega)$, as the radial derivative is singular at the origin. We test the convergence of the R-FEM approximation under a uniform $h$-refinement. We also take the opportunity to test the a posteriori estimate given in Theorem 5.1. Results are summarised in Figure 7, where convergence plots for $r=0,s=1$ are presented.

7.3 Test 3 – Mesh adaptivity

We continue with testing the behaviour of the estimator within a standard adaptive algorithm of the form $\text{SOLVE}\to\text{ESTIMATE}\to\text{MARK}\to\text{REFINE}.$ In this context, over a coarse mesh we solve the R-FEM problem, estimate the error using the a posteriori error bound from Theorem 5.1, mark a subset of elements where the a posteriori indicator is large, and refine those elements. The algorithm is then iterated until a prescribed tolerance is achieved. For the marking step we make use of the maximum strategy, that is we mark a set of elements $\mathcal{M}\subset\mathcal{T}_{h}$ such that

[TABLE]

where $\theta\in(0,1]$ is the ratio of refined elements. In our experiments we choose $\theta=0.25$. The numerical results are summarised in Figure 8, where convergence plots for the lowest order ($r=0,s=1$) adaptive R-FEM approximation to the solution of the above problem are shown, indicating that the adaptive approximation is optimally convergent and its effectivity index, i.e., the ratio between the estimator and the exact error, is small and remains bounded.

7.4 Test 4 – R-FEM for convection-diffusion problems

To showcase the versatility of the R-FEM framework, we consider a convection-diffusion problem solved by a ‘hybrid’ R-FEM discretisation for the second order term and an upwinded dG method for the first order term. We envisage such a combination of R-FEM with stabilised methods for stiff problems as a typical area of applications for R-FEM, aiming to combined the advantages of conforming (or classical non-conforming) discretisations with considerable freedom in the treatment of lower order terms. To this end, we consider the convection-diffusion problem

[TABLE]

where $\mathbb{R}\ni\epsilon\ll 1$ represents viscosity and $\boldsymbol{w}\in C(\bar{\Omega})^{d}$ and $c\in L^{2}(\Omega)$ are some prescribed convection vector field and reaction coefficients, respectively. We supplement (37) with homogeneous Dirichlet (no-slip) boundary conditions on $\partial\Omega$.

The flexibility in the design of R-FEM described above, which employs element-wise discontinuous discrete spaces, allows naturally the treatment of the diffusion term in a conforming manner and, simultaneously, the convection term in an upwinded nonconforming fashion. This allows for the construction of stable methods in the convection-dominated regime. To this end, consider the R-FEM discretisation: find $u_{h}\in V^{0}_{h}$ (with $\mathcal{E}(u_{h})\in V^{1}_{h}\cap H^{1}_{0}(\Omega)$), such that:

[TABLE]

where \lfloor v\rfloor(x):=\lim_{\delta\to 0^{+}}\big{(}v(x+\delta\boldsymbol{w}(x))-v(x-\delta\boldsymbol{w}(x))\big{)}, for (almost every) $x\in\partial T$, is the upwind jump and

[TABLE]

denotes the inflow boundary of an element $T\in\mathcal{T}$; when $x\in\partial\Omega$, we set the respective downwind value to zero in the definition of $\lfloor\cdot\rfloor$.

We shall investigate the convergence and stability properties of the proposed method (38) for $\Omega=(0,1)^{2}$ by setting

[TABLE]

and vary the diffusion parameter $\epsilon$ throughout the tests. Using the convection field $(a)$, we test the method’s convergence properties for $\sigma$ given as in (24) and $\mathcal{A}=\epsilon$; the results are summarised in Figure 9.

To test the stability of the R-FEM approach, we consider (37) with the diagonal convection coefficient $(b)$ for small values of $\epsilon=10^{-2},10^{-3}$ on a quasi-uniform mesh with $h\approx 5\times 10^{-3}$. The exact solution to this problem admits a boundary layer of scale $\epsilon$ on the top and right boundaries. Therefore, the mesh is sufficiently fine to resolve the layer for $\epsilon=10^{-2}$ but not for $\epsilon=10^{-3}$. We compare the performance of R-FEM against a standard stable method, namely the (upwinded) interior penalty dG method which reads: find $u_{h}\in V^{1}_{h}$ such that

[TABLE]

with $\sigma$ given by (22). Some representative solutions are given in Figure 10, using the recovered finite element method (38) (left plots) and the dG method (41) (right plots), respectively. We observe that R-FEM produces stable solution even in the non-resolving regime $\epsilon=10^{-4}$ including smooth profiles on the boundary layer. On the other hand, as expected, the dG solution is stable away from the layer and oscillates in the vicinity of the boundary layer. We can see, therefore, that R-FEM has the potential in delivering stable, conforming discretisations for stiff PDE problems.

8 Outlook and further applications

In an effort to combine the advantages of discontinuous finite element spaces, yet produce conforming or classical non-conforming discretisations, we have introduced the framework of recovered finite element methods. We have demonstrated the convergence properties of R-FEM on standard simplicial and, in some instances, also on box-type meshes, we have proven an posteriori error bound and we have assessed the conditioning of the resulting stiffness matrices. As a first work on describing the R-FEM framework, we have refrained from describing advanced potential areas of applications in detail, in an effort to focus on the key ideas. Nonetheless, R-FEM is envisaged to be successfully applicable to a number of directions. Here, we shall briefly comment on some of these.

8.1 R-FEM for polytopic meshes

An important attribute of the definition of R-FEM (5) is that the spaces $V_{h}^{r}$ and $V_{h}^{s}$ can differ. This idea has been used above in the context of recovering into different local polynomial degrees, yet on the same mesh, aiming to improve on the convergence rate and/or the error per degree of freedom compared to respective conforming FEMs. However, it is by all means possible to consider spaces $V_{h}^{r}$ and $V_{\tilde{h}}^{s}$ subordinate to different meshes. One natural application of this idea is the construction of R-FEM on general polygonal/polyhedral (termed collectively henceforth as polytopic) meshes, based on a recovery operator $\mathcal{E}:V_{h}^{r}\to V_{\tilde{h}}^{s}\cap H^{1}_{0}(\Omega)$, whereby $V_{h}^{r}$ is a discontinuous space subordinate to a polytopic mesh and $V_{\tilde{h}}^{s}\cap H^{1}_{0}(\Omega)$ is a conforming space subordinate to a (perhaps simplicial) finer mesh. This way once can construct conforming approximations of elliptic problems on polytopic meshes where one has access to the whole of the approximate solution; this is in contrast to the recent virtual element framework [5] whereby only certain functionals of the approximate solution are available to the user. This is a significant topic in its own right and will be discussed in detail in [17].

8.2 R-FEM for high order PDEs

The availability of classical $H^{2}$-conforming (e.g., Argyris or Hsieh-Clough-Tocher elements [12])) has allowed for the construction of respective recovery operators $\mathcal{E}(w)\in H^{2}(\Omega)$ for $w\in V_{h}^{r}$ or $w\in V_{h}^{r}\cap H^{1}(\Omega)$ [7, 19]. Therefore, it is possible to construct R-FEM methods for biharmonic/plate problems. As this is also a significant topic in its own right, it will be discussed elsewhere.

8.3 Mixed R-FEM for problems with constraints

Upon constructing suitable recovery operators, it is interesting to revisit the classical problem of numerical approximation of PDE problems with constraints, e.g., (nearly) incompressible elasticity, Darcy or viscous incompressible flows. Indeed, the flexibility in the underlying degrees of freedom offered by R-FEM in conjunction with suitable recovery operators, e.g., recovering onto classical inf-sup stable mixed finite element pairs is a potentially substantial direction of further research. For some first results in this context, we refer to [4].

Acknowledgements

We wish to express our sincere gratitude to Gabriel R. Barrenechea (University of Strathclyde), Andrea Cangiani and Zhaonan (Peter) Dong (University of Leicester) for their insightful comments on an earlier version of this work, which lead to a clearer understanding of the method. EHG acknowledges funding by The Leverhulme Trust.