# The eXtended Virtual Element Method for elliptic problems with weakly singular solutions

Jérôme Droniou*¹,², Gianmarco Manzini†³, and Liam Yemm‡²

¹IMAG, Univ. Montpellier, CNRS, Montpellier, France

²School of Mathematics, Monash University, Clayton, Australia

³T-5, Theoretical Division, Los Alamos National Laboratory, Los Alamos, NM, USA

June 17, 2024

## Abstract

This paper introduces a novel eXtended virtual element method, an extension of the conforming virtual element method. The X-VEM is formulated by incorporating appropriate enrichment functions in the local spaces. The method is designed to handle highly generic enrichment functions, including singularities arising from fractured domains. By achieving consistency on the enrichment space, the method is proven to achieve arbitrary approximation orders even in the presence of singular solutions. The paper includes a complete convergence analysis under general assumptions on mesh regularity, and numerical experiments validating the method's accuracy on various mesh families, demonstrating optimal convergence rates in the $L^2$- and $H^1$-norms on fractured or L-shaped domains.

**Key words:** X-VEM, polytopal method, corner singularities, fractured domains, error analysis, enriched method

## 1 Introduction

Introduced by [2, 30], the conforming virtual element method (VEM) is a modern numerical technique for the discretisation of partial differential equations. At its core, VEM transcends traditional finite element methods by allowing for a flexible and versatile discretisation strategy, wherein the computational domain is partitioned into polytopes (polygons in 2D, polyhedra in 3D) of arbitrary shapes. The VEM is one of several polytopal methods developed in the last couple of decades or so, which also include Discontinuous Galerkin methods [3, 14, 17], Hybrid Discontinuous Galerkin [15, 16], Hybrid High-Order (HHO) [18, 20], and Weak Galerkin [29], among others. We also note that some of these methods drew inspiration and/or turned out to be arbitrary-order methods of older low-order schemes, in particular from the Mimetic Finite

<footnote>*jerome.droniou@umontpellier.fr, jerome.droniou@monash.edu</footnote>

<footnote>†gmanzini@lanl.gov</footnote>

<footnote>‡liam.yemm@monash.edu</footnote>
Differences or Finite Volume families [21–23, 25, 31]. VEM’s distinctive feature lies in its use of a ‘virtual space’ for which basis functions are not known explicitly, avoiding the need for predefined shape functions.

Errors estimates of arbitrary-order methods, including VEM, are typically limited by the regularity of the exact solution. On non-smooth domains (such as regions with non-convex corners or those possessing cracks) it is expected that the exact solution to elliptic problems will contain weak singularities [24]. Singular behaviour can also stem from the problem data, for example, the boundary conditions in the Motz problem [28]. This lack of regularity is well documented in the finite element literature and is typically overcome through enriched approximations based on a partition of unity method [6, 26]. The extended finite element method [8, 27] is one such method, originally designed to handle discontinuities in crack growth models. In particular, by enriching the local spaces with basis functions that are discontinuous across the crack, the method allows for optimal approximation without the need for mesh refinement near the discontinuity. Finite element methods, whether standard or enriched, are restricted to specific meshes (mostly composed of triangles/tetrahedra or quadrangles/hexahedra), which limits their flexibility. For instance, preserving mesh regularity during refinement necessitates a large number of elements and layers, hindering local and adaptive refinements. Similarly, mesh agglomeration (e.g., for multigrid solvers) and mesh cutting (e.g., for domains with interfaces) result in admissible finite element meshes only in very specific, often trivial, cases. For these reasons, polytopal methods have received a lot of attention since the early 2000s.

To blend the benefits of polytopal methods (mesh flexibility) together with enriched finite elements (optimal convergence for solutions with singularities), a trend in the literature recently arose around designing and analysing *enriched polytopal methods*. The work of [5] designs an enriched non-conforming virtual element method (NC-VEM) for the Poisson problem targeting singularities arising from the domain geometry (such as corners or cracks). However, the analysis there requires to assume a discrete trace inequality, which depends on the enrichment function. Following this approach, the same authors have defined an enriched NC-VEM for a plane elasticity problem with corner singularities [4]. An extended HHO method has been designed for the Poisson problem [33] targeting general enrichments with (locally) square-integrable Laplacian and Neumann traces and avoids the use of discrete inequalities dependent on the enrichment function. We also make note of the recent work of [32] which designs an extended HHO scheme for the Stokes problem with generic enrichment functions, albeit with limited regularity disallowing singularities arising from non-convex corners. To date, only [9, 10] have attempted to design an enriched *conforming* virtual element method. However, both these methods only consider the lowest order VEM and lack a complete error analysis, mostly because of an approach based on the continuous enriched virtual space and the lack of inverse inequalities for enriched functions or degrees of freedom. In this paper, an eXtended virtual element method (X-VEM) is designed which is valid for highly generic enrichment functions, offers arbitrary approximation orders and is accompanied with a rigorous and complete analysis; to overcome the shortcomings of previous (partial) analyses mentioned above, we adopt a “fully discrete” approach closer to the approach used, e.g., for HHO schemes. Moreover, and contrary to the previously mentioned enriched NC-VEM, the method is not specific to harmonic singularities or those arising from cracks or corners in the domain, but is valid for a generic enrichment space satisfying Assumption 2.1 below. This assumption is still valid for singularities arising from fractured domains – something the extended HHO method fails to capture.

We perform the convergence analysis under suitable but still general assumptions on the mesh regularity, and we derive the estimates for the approximation error. We obtain error estimates that depend solely on the regular component of the continuous solution, both in a discrete $H^1$-like norm based on the degrees of
freedom and in the (computable) broken $H^1$-norm using the enriched elliptic reconstruction. To assess the
method's behaviour and confirm the theoretical expectations, we conduct a set of numerical experiments
on several mesh families, including meshes with convex and non-convex elements on the $L$-shaped domain
with a corner singularity. Our results assess the accuracy of the extended virtual element method and demon-
strate its optimal rates of convergence in the $L^2$- and $H^1$-norm. As is well-known for singularity-enriched
schemes, to avoid ill-conditioning of the system a local enrichment is preferable, through which only the
DOFs in a portion of the domain are enriched with the singular function; contrary to the case of HHO en-
riched methods [32, 33] our analysis does not cover local enrichment, which remains for a future work, but
we numerically demonstrate that it does work for the scheme developed here.

While an extension of the method designed in this paper to three-dimensions may be feasible, signif-
icant investigation into the design of virtual spaces on element boundaries and appropriate reconstruction
operators would be required.

The outline of the paper is as follows. In Section 2, we introduce the model problem and its extended
virtual element approximation (X-VEM). In Section 3, we present the convergence analysis of X-VEM. The
exposition is done in dimension 2 for simplicity, but can be generalised to dimension 3 in a straightforward
manner. In Section 4 we investigate the method’s performance on suitable numerical experiments. Finally,
in Section 5 we offer some concluding remarks.

**1.1 Notation and technicalities**

Throughout this paper, we adopt the notation of Sobolev spaces of [1]. Accordingly, we denote the space of
square-integrable functions defined on any open, bounded, connected domain $\mathcal{D} \subset \mathbb{R}^2$ with boundary $\partial \mathcal{D}$
by $L^2(\mathcal{D})$, and the Hilbert space of functions in $L^2(\mathcal{D})$ with all partial derivatives up to a positive integer $m$
also in $L^2(\mathcal{D})$ by $H^m(\mathcal{D})$. The norm of $L^2(\mathcal{D})$ is written $\|\cdot\|_{\mathcal{D}}$, while the norm and semi-norm in $H^m(\mathcal{D})$
are respectively denoted by $\|\cdot\|_{H^m(\mathcal{D})}$ and $|\cdot|_{H^m(\mathcal{D})}$; the latter is the sum of $L^2$-norms of derivatives of order
$m$.

Let $\Omega \subset \mathbb{R}^2$ be an open, bounded polygonal domain with boundary $\Gamma$. The virtual element method is
formulated on the mesh family $\{\Omega_h\}_h$, where each mesh $\Omega_h$ is a partition of the computational domain
$\Omega$ into non-overlapping polygonal elements $E$. A polygonal element $E$ is a compact subset of $\mathbb{R}^2$ with
boundary $\partial E$, area $|E|$ centre of mass $\mathbf{x}_E$, and diameter $h_E = \sup_{\mathbf{x},\mathbf{y} \in E} |\mathbf{x} - \mathbf{y}|$. The mesh elements $\Omega_h$
form a finite cover of $\Omega$ such that $\bar{\Omega} = \cup_{E \in \Omega_h} E$ and the mesh size labelling each mesh $\Omega_h$ is defined by
$h = \max_{E \in \Omega_h} h_E$. A mesh edge $e$ has centre $\mathbf{x}_e$ and length $h_e$; a mesh vertex $v$ has position vector $\mathbf{x}_v$. We
denote the set of mesh edges by $\mathcal{E}_h$ and the set of mesh vertices by $\mathcal{V}_h$. We denote the set of mesh edges by
$\mathcal{E}_h$ and the set of mesh vertices by $\mathcal{V}_h$.

The following regularity assumption is made on the mesh $\Omega_h$. We note that this assumption allows for
arbitrarily small (or numerous) edges.

**Assumption 1.1 (Regular mesh sequence)** There exists a constant $\rho > 0$ such that every $E \in \Omega_h$ is star-
shaped with respect to a ball of radius $\rho h_E$.

For any integer number $\ell \ge 0$, we let $\mathbb{P}_\ell(E)$ and $\mathbb{P}_\ell(e)$ denote the space of polynomials defined on the
element $E$ and the edge $e$, respectively; $\mathbb{P}_\ell(\Omega_h)$ denotes the space of piecewise polynomials of degree $\ell$ on
the mesh $\Omega_h$.
# 2 The eXtended virtual element method

## 2.1 Model problem

We consider the Poisson problem with homogeneous Dirichlet boundary conditions

$$-\Delta u = f \quad \text{in } \Omega, \qquad (1a)$$

$$u = 0 \quad \text{on } \Gamma, \qquad (1b)$$

for the scalar unknown $u$. The extension to non-homogeneous Dirichlet boundary conditions or to different types of boundary conditions such as Neumann or Robin conditions is straightforward.

Let $H_0^1(\Omega)$ be the subspace of the Sobolev space of functions $H^1(\Omega)$ with zero trace on $\Gamma$. The variational formulation of problem (1) reads as: Find $u \in H_0^1(\Omega)$ such that

$$a(u, v) = \int_{\Omega} f v \, dx \quad \forall v \in H_0^1(\Omega), \qquad (2)$$

where the bilinear form $a(\cdot, \cdot)$ is given by

$$a(u, v) = \int_{\Omega} \nabla u \cdot \nabla v \, dx. \qquad (3)$$

The continuity and coercivity of the bilinear form $a(\cdot, \cdot)$ imply the existence and uniqueness of the solution according to the Lax-Milgram lemma.

## 2.2 Enrichment space

We assume that the exact solution of the variational problem (2) is the sum of two terms

$$u = u^r + \psi,$$

where $u^r$ is sufficiently smooth and $\psi$ is the (weakly) singular part. We are interested in developing an extended virtual element method that incorporates $\psi$ in the design of the approximation space, in such a way that the approximation properties of the scheme only depend on the regular component $u^r$. The singular component $\psi$ is assumed to be an element of a finite dimensional enrichment space $\Psi(\Omega)$ satisfying the following regularity requirements.

**Assumption 2.1 (Regularity of the enrichment space)**

1. $\Psi(\Omega) \subset H^1(\Omega) \cap C^0(\bar{\Omega})$.

2. $\Delta\Psi(E) \subset L^2(E)$ for all $E \in \Omega_h$.

Here and in the following, if $U \subset \bar{\Omega}$ we set $\Psi(U) := \{w|_U : w \in \Psi(\Omega)\}$.
**Remark 2.2 (Conforming enrichment space and local enrichment)** The requirement that the enrichment space is globally conforming is due to the fact that the analysis is based on splitting the X-VEM interpolant into a singular component and a regular VEM interpolant (see Theorem 2.9 below). If the enrichment space is not contained in $H^1(\Omega)$, then by writing $u = u^r + \psi \in H^1(\Omega)$ it cannot hold that $u^r \in H^1(\Omega)$ and thus its regular VEM interpolant is not defined. This has the major drawback of not allowing for enrichment spaces defined piece-wise – a typical method of local enrichment, see e.g. [33] in the context of HHO methods. However, numerical tests suggest that locally enriching in this manner yields valid results (see Section 4).

An alternative, and perhaps more natural way to define the X-VEM interpolant, is as the element of the X-VEM space with the same degrees of freedom as the function being interpolated. This would no longer require splitting a function into its regular and singular components and interpolating them individually, thus allowing for the weaker assumption $\Psi(E) \subset H^1(E) \cap C^0(E)$ for all $E \in \Omega_h$. However, proving optimal approximation properties of an interpolant defined in this manner seems difficult; see Remark 2.10.

## 2.3 X-VEM spaces and elliptic projector

Consider an integer $k \ge 1$ and set $l = \max\{0, k - 2\}$. The following extended polynomial spaces are defined on each element $E \in \Omega_h$,

$$ \mathbb{P}_k^\Psi(E) := \mathbb{P}_k(E) + \Psi(E), \quad \mathbb{P}_l^\Delta(E) := \mathbb{P}_l(E) + \Delta\Psi(E), $$

and the following space is defined on each edge $e \in \mathcal{E}_h$,

$$ \mathbb{P}_k^\Psi(e) := \mathbb{P}_k(e) + \Psi(e). $$

The local X-VEM space is defined as

$$ V_{k,h}^\Psi(E) := \{ v_h \in H^1(E) : \Delta v_h \in \mathbb{P}_l^\Delta(E), \ v_{h|_{\partial E}} \in C^0(\partial E), \ v_{h|_e} \in \mathbb{P}_k^\Psi(e) \ \forall e \subset \partial E \}. $$

We note that $V_{k,h}^\Psi(E)$ contains both the regular local VEM space of [30] and the extended polynomial space $\mathbb{P}_k^\Psi(E)$. For each $e \in \mathcal{E}_h$, we denote by $\mathfrak{P}_e$ the $L^2$-orthogonal complement of $\mathbb{P}_k(e)$ in $\mathbb{P}_k^\Psi(e)$, so that $\mathbb{P}_k^\Psi(E) = \mathbb{P}_k(e) \oplus \mathfrak{P}_e$, and we set $\mathbb{P}_{k-2}^\Psi(e) = \mathbb{P}_{k-2}(e) \oplus \mathfrak{P}_e$. Heuristically, $\mathfrak{P}_e$ represents the “most singular” component of $\mathbb{P}_k^\Psi(e)$, obtained by only keeping the part of $\Psi(e)$ that is orthogonal to $\mathbb{P}_k(e)$ (that is, the part of $\Psi(e)$ which is the furthest away from polynomials of degree $\le k$).

The $L^2$-orthogonal projectors onto each of the spaces $\mathbb{P}_l^\Delta(E)$ and $\mathbb{P}_{k-2}^\Psi(E)$ are denoted by $\Pi_{l,E}^\Delta$ and $\Pi_{k-2,e}^\Psi$, respectively. We show in Lemmas 2.3 and 2.6 below that each $v_h \in V_{k,h}^\Psi(E)$ is uniquely characterised by the following degrees of freedom:

(D1) the values of $v_h$ at the vertices of $E$;

(D2) the $L^2$-orthogonal projection $\Pi_{k-2,e}^\Psi v_h$ of $v_h$ onto the space $\mathbb{P}_{k-2}^\Psi(e)$ for each edge $e \subset \partial E$;

(D3) the $L^2$-orthogonal projection $\Pi_{l,E}^\Delta v_h$ of $v_h$ onto the space $\mathbb{P}_l^\Delta(E)$.

**Lemma 2.3 (Unisolvence of boundary values)** The degrees of freedom (D1)–(D2) are unisolvent for the trace space of $V_{k,h}^\Psi(E)$ on $\partial E$.
*Proof.* We have to show that, for any $v_h \in V_{k,h}^{\Psi}(E)$, the degrees of freedom (D1)–(D2) entirely determine any $v_h|_{\partial E}$, and that for any choice of the values (D1)–(D2), we can find $w \in C^0(\partial E)$ with these degrees of freedom such that $w|_e \in \mathbb{P}_k^{\Psi}(e)$ for each $e \subset \partial E$.

Let $v_h \in V_{k,h}^{\Psi}(E)$. To show that its boundary value is entirely determined by (D1)–(D2), we only have to show that, for each $e \in \mathcal{E}_h$, $v_e := v_h|_e$ is uniquely determined by its values at the endpoints $\partial e$ of $e$ and by $\Pi_{k-2,e}^{\mathfrak{F}}v_e$. Since $v_e \in \mathbb{P}_k^{\Psi}(e) = \mathbb{P}_k(e) \oplus \mathfrak{P}_e$ we can write $v_e = q + z$ with $q \in \mathbb{P}_k(e)$ and $z \in \mathfrak{P}_e$. We first note that, by orthogonality of $\mathbb{P}_k(e)$ and $\mathfrak{P}_e$, $\Pi_{k-2,e}^{\mathfrak{F}}v_e = \Pi_{k-2}^{0,e}q + z$, where $\Pi_{k-2}^{0,e}$ is the $L^2$-orthogonal projection on $\mathbb{P}_{k-2}(e)$. Hence, letting $\Pi_e^{0,\mathfrak{F}}$ be the $L^2$-orthogonal projection on $\mathfrak{P}_e$, we have

$$ \Pi_e^{0,\mathfrak{F}}(\Pi_{k-2,e}^{\mathfrak{F}}v_e) = \Pi_e^{0,\mathfrak{F}}(\Pi_{k-2}^{0,e}q + z) = z, \quad (4) $$

where the conclusion follows from the orthogonality of $\mathbb{P}_{k-2}(e)$ and $\mathfrak{P}_e$.
We then recall that $v_e = q + z$ to write

$$ \Pi_{k-2}^{0,e}q = \Pi_{k-2}^{0,e}(v_e - z) = \Pi_{k-2}^{0,e}(\Pi_{k-2,e}^{\mathfrak{F}}v_e), \quad (5) $$

where the second equality is obtained using $\mathbb{P}_{k-2}(e) \subset \mathbb{P}_{k-2}^{\mathfrak{F}}(e)$ (which ensures that $\Pi_{k-2}^{0,e} = \Pi_{k-2}^{0,e} \circ \Pi_{k-2,e}^{\mathfrak{F}}$) and the orthogonality of $z$ and $\mathbb{P}_{k-2}(e)$. Moreover

$$ q(a) = v(a) - z(a) \quad \forall a \in \partial e. \quad (6) $$

The relation (4) shows that

$$ z \text{ is uniquely determined by } \Pi_{k-2,e}^{\mathfrak{F}}v_e. \quad (7) $$

Since $q \in \mathbb{P}_k(e)$, it is uniquely determined by $(q(a))_{a \in \partial e}$ and $\Pi_{k-2}^{0,e}q$, that is, according to (5)–(6), by $(v(a) - z(a))_{a \in \partial e}$ and $\Pi_{k-2,e}^{\mathfrak{F}}v_e$; by (7), we infer that

$$ q \text{ is uniquely determined by } (v(a))_{a \in \partial e} \text{ and } \Pi_{k-2,e}^{\mathfrak{F}}v_e. \quad (8) $$

Combining (7) and (8) shows that $v_e = q + z$ is uniquely determined by $(v(a))_{a \in \partial e}$ and $\Pi_{k-2,e}^{\mathfrak{F}}v_e$, which is what we needed to prove.

Conversely, take $(m_a)_{a \in \nu_h \subset \partial E}$ and $(r_e)_{e \in \mathcal{E}_h \cap \partial E}$ values for (D1)–(D2) (each $m_a$ belongs to $\mathbb{R}$, each $r_e$ belongs to $\mathbb{P}_{k-2}^{\mathfrak{F}}(e)$). Define $w : \partial E \to \mathbb{R}$ the following way; for each $e \subset \partial E$, set $z_e = \Pi_e^{0,\mathfrak{F}}r_e \in \mathfrak{P}_e$, define $q_e \in \mathbb{P}_{k-2}(e)$ such that $\Pi_{k-2}^{0,e}q_e = \Pi_{k-2}^{0,e}r_e$ and $q_e(a) = m_a - z_e(a)$ for each $a \in \partial e$, and set $w|_e = q_e + z_e$. Then it can easily be checked that $w \in C^0(\partial E)$ since its value at each vertex $a \in \partial E$ is $m_a$ (from either side of $a$), and the arguments developed above to establish (4)–(6) show that its degrees of freedom (D1)–(D2) match the chosen values. $\square$

**Remark 2.4 (Choice of boundary degrees of freedom)** A seemingly more natural choice for (D2) would be to consider the projection of $v_h$ on $\mathbb{P}_{k-2}(e) + \Psi(e)$ for each $e \subset \partial E$. However, such a choice does not seem to ensure the unisolvence stated in Lemma 2.3.
**Remark 2.5 (Computation of the boundary degrees of freedom)** A basis for the space $\mathfrak{P}_e$ can be generated by considering a canonical basis for $\mathbb{P}_k(e) + \Psi(e)$, orthonormalising, and retrieving the degrees of freedom corresponding to $\mathfrak{P}_e$. Consequently, a basis for the space $\mathbb{P}_{k-2}^\mathfrak{P}(e)$ follows. This procedure provides orthonormal bases of local spaces, which are known to positively impact the conditioning of the global algebraic system; see, e.g., the analysis in [7] for the HHO method.

**Lemma 2.6 (Unisolvence)** For all elements $E \in \Omega_h$, the values provided by the continuous linear functionals (D1), (D2), (D3) are unisolvent in the virtual element space $V_{k,h}^\Psi(E)$.

*Proof.* Let $V_{k,h}^\Psi(\partial E)$ be the trace space on $\partial E$ of $V_{k,h}^\Psi(E)$. The mapping

$$V_{k,h}^{\Psi}(\partial E) \times \mathbb{P}_{l}^{\Delta}(E) \mapsto V_{k,h}^{\Psi}(E) \\ (w, r) \rightarrow v_h \text{ such that } -\Delta v_h = r \text{ in } E \text{ and } v_h|_{\partial E} = w$$

is an isomorphism, by definition of $V_{k,h}^\Psi(E)$ and the well-posedness of the Poisson problem. Hence, the dimension of $V_{k,h}^\Psi(E)$ is equal to the dimension of $V_{k,h}^\Psi(\partial E) \times \mathbb{P}_l^\Delta(E)$ which, by Lemma 2.3, is identical to the number of degrees of freedom (D1)–(D3).

Therefore, it remains to prove that if all DOFs of $v_h \in V_{k,h}^\Psi(E)$ vanish, then $v_h = 0$. By Lemma 2.3, we already know that if (D1) and (D2) vanish, then $v_h = 0$ on $\partial E$. Therefore, an integration by parts reveals

$$\int_E \nabla v_h \cdot \nabla v_h \, d\boldsymbol{x} = - \int_E v_h \cdot \Delta v_h \, d\boldsymbol{x} = - \int_E \Pi_{l,E}^\Delta v_h \cdot \Delta v_h \, d\boldsymbol{x} = 0,$$

where the last equality follows from $\Pi_{l,E}^\Delta v_h = 0$ as a result of (D3) vanishing. Therefore, $v_h = \text{const}$, and this constant must be zero due to $v_h = 0$ on $\partial E$. $\square$

The extended elliptic projector $\Pi_{k,E}^{\nabla, \Psi}: V_{k,h}^{\Psi}(E) \to \mathbb{P}_k^{\Psi}(E)$ projects the extended VEM space $V_{k,h}^{\Psi}(E)$ onto the extended polynomial space $\mathbb{P}_k^{\Psi}(E)$. The projection $\Pi_{k,E}^{\nabla, \Psi} v_h$ is the solution to the variational problem:

$$\int_E \nabla \Pi_{k,E}^{\nabla, \Psi} v_h \cdot \nabla q \, d\boldsymbol{x} = \int_E \nabla v_h \cdot \nabla q \, d\boldsymbol{x} \quad q \in \mathbb{P}_k^{\Psi}(E), \qquad (9)$$

$$\int_E \Pi_{k,E}^{\nabla, \Psi} v_h \, d\boldsymbol{x} = \int_E v_h \, d\boldsymbol{x}. \qquad (10)$$

Clearly, $\Pi_{k,E}^{\nabla, \Psi} q = q$ for each $q \in \mathbb{P}_k^{\Psi}(E)$. For all mesh elements $E$ and all virtual element functions $v_h \in V_{k,h}^\Psi(E)$, the extended projection $\Pi_{k,E}^{\nabla, \Psi} v_h$ is computable from the DOFs as stated and proved in the following lemma.

**Lemma 2.7 (Computability of the elliptic projector)** Let $E$ be an element of $\Omega_h$ and $v_h \in V_{k,h}^\Psi(E)$. Then, the extended elliptic projection $\Pi_{k,E}^{\nabla, \Psi} v_h \in \mathbb{P}_k^{\Psi}(E)$ is computable using only the degrees of freedom (D1)–(D3) of $v_h$.
Proof. It holds, by an integration by parts, that for all $q \in \mathbb{P}_k^\Psi(E)$

$$ \int_E \nabla \Pi_{k,E}^{\nabla, \Psi} v_h \cdot \nabla q \, d\boldsymbol{x} = - \int_E v_h \Delta q \, d\boldsymbol{x} + \langle \boldsymbol{n}_E \cdot \nabla q, v_h \rangle_{-\frac{1}{2}, \frac{1}{2}, \partial E} \quad (11) $$

where $\langle \cdot, \cdot \rangle_{-\frac{1}{2}, \frac{1}{2}, \partial E}$ denotes the duality product between $H^{-1/2}(\partial E)$ and $H^{1/2}(\partial E)$. The moments of $v_h$ against $\Delta q \in \mathbb{P}_l^\Delta(E)$ are known from (D3) and $v_h$ is known entirely on the boundary $\partial E$ from (D1)-(D2) (see Lemma 2.3. Therefore, all the terms in the right-hand side of (11) are computable, and so is $\nabla \Pi_{k,E}^{\nabla, \Psi} v_h$. Since $\mathbb{P}_l^\Delta(E)$ contains the constant functions, the integral of $v_h$ over $E$ (and thus that of $\Pi_{k,E}^{\nabla, \Psi} v_h$) is computable from (D3). The entire elliptic projector is therefore computable from the DOFs. $\square$

The global X-VEM space is constructed by patching the local spaces:

$$ V_{k,h}^{\Psi} := \{ v_h \in H^1(\Omega) : v_{h|E} \in V_{k,h}^{\Psi}(E) \quad \forall E \in \Omega_h \}. $$

Each $v_h \in V_{k,h}^\Psi$ is uniquely defined by the following degrees of freedom:

(D1) the values of $v_h$ at each vertex of $\mathcal{V}_h$;

(D2) the $L^2$-orthogonal projection $\Pi_{k-2,e}^{\mathfrak{F}} v_h$ of $v_h$ onto the space $\mathbb{P}_{k-2}^{\mathfrak{F}}(e)$ for each edge $e \in \mathcal{E}_h$;

(D3) the $L^2$-orthogonal projection $\Pi_{l,E}^{\Delta} v_h$ of $v_h$ onto the space $\mathbb{P}_l^\Delta(E)$ for each element $E \in \Omega_h$.

## 2.4 Formulation of the scheme and main results

At this point, the construction of the X-VEM is straightforward and follows the usual procedural steps [30]. We define the discrete bilinear form $a_h : V_{k,h}^\Psi \times V_{k,h}^\Psi \to \mathbb{R}$ by assembling all elemental contributions

$$ a_h(u_h, v_h) := \sum_{E \in \Omega_h} \left[ a^E(\Pi_{k,E}^{\nabla, \Psi} u_h, \Pi_{k,E}^{\nabla, \Psi} v_h) + S^E(u_h, v_h) \right], \quad (12) $$

where

$$ a^E(u, v) = \int_E \nabla u \cdot \nabla v \, d\boldsymbol{x} $$

and the stabilisation term $S^E(\cdot, \cdot) : V_{k,h}^\Psi(E) \times V_{k,h}^\Psi(E) \to \mathbb{R}$ can be any computable (from the DOFs), symmetric, positive semi-definite bilinear form satisfying:

* Coercivity and boundedness on $\ker(\Pi_{k,E}^{\nabla, \Psi})$. For all $v_h \in V_{k,h}^\Psi(E) \cap \ker(\Pi_{k,E}^{\nabla, \Psi})$:

$$ h_E^{-2} \| \Pi_{l,E}^{\Delta} v_h \|_E^2 + h_E^{-1} \| v_h \|_{\partial E}^2 \lesssim S^E(v_h, v_h) \lesssim h_E^{-2} \| \Pi_{l,E}^{\Delta} v_h \|_E^2 + h_E^{-1} \| v_h \|_{\partial E}^2. \quad (13) $$

* Consistency on $\mathbb{P}_k^\Psi(E)$. For all $w_h \in \mathbb{P}_k^\Psi(E)$ and $v_h \in V_{k,h}^\Psi(E)$

$$ S^E(w_h, v_h) = 0. \quad (14) $$
Here and throughout the rest of this paper, the notation $a \lesssim b$ is used if there exists a constant $C > 0$ independent of the quantities $a, b$ and the mesh size $h$ such that $a \le Cb$. The hidden constant depends only on $\Omega, k$ and the mesh regularity $\rho$ (see Assumption 1.1 below).

Defining the local seminorm $\|\cdot\|_{1,E} : V_{k,h}^{\Psi}(E) \to \mathbb{R}$ via

$$ \|v_h\|_{1,E}^2 := \|\nabla \Pi_{k,E}^{\nabla, \Psi} v_h\|_E^2 + h_E^{-2} \|\Pi_{l,E}^{\Delta} (v_h - \Pi_{k,E}^{\nabla, \Psi} v_h)\|_E^2 + h_E^{-1} \|v_h - \Pi_{k,E}^{\nabla, \Psi} v_h\|_{\partial E}^2, \quad (15) $$

and the global seminorm $\|\cdot\|_{1,h} : V_{k,h}^{\Psi} \to \mathbb{R}$ via

$$ \|v_h\|_{1,h}^2 := \sum_{E \in \Omega_h} \|v_h\|_{1,E}^2, \quad (16) $$

we infer from the definition (12) of $a_h$ and (13) that

$$ \|v_h\|_{1,h}^2 \lesssim a_h(v_h, v_h) \lesssim \|v_h\|_{1,h}^2. \quad (17) $$

We prove in Section 3 that $\|\cdot\|_{1,h}$ (and, thus, $a_h(\cdot, \cdot)^{1/2}$) defines a norm on the homogeneous subspace

$$ V_{k,h,0}^{\Psi} := \{v_h \in V_{k,h}^{\Psi} : v_h|_{\partial \Omega} = 0\}. $$

**Remark 2.8 (Example of stabilisation and choice of local semi-norm)** An example of stabilisation bilinear form satisfying (13) and (14) is given by:

$$ S^E(u_h, v_h) = h_E^{-2} \int_E \left( \Pi_{l,E}^{\Delta} (u_h - \Pi_{k,E}^{\nabla, \Psi} u_h) \right) \left( \Pi_{l,E}^{\Delta} (v_h - \Pi_{k,E}^{\nabla, \Psi} v_h) \right) dx \\ + h_E^{-1} \int_{\partial E} \left( u_h - \Pi_{k,E}^{\nabla, \Psi} u_h \right) \left( v_h - \Pi_{k,E}^{\nabla, \Psi} v_h \right) ds. \quad (18) $$

Indeed, if $v_h \in \ker(\Pi_{k,E}^{\nabla, \Psi})$ then $S^E(v_h, v_h) = h_E^{-2} \|\Pi_{l,E}^{\Delta} v_h\|_E^2 + h_E^{-1} \|v_h\|_{\partial E}^2$ so (13) holds and, if $w_h \in \mathbb{P}_k^{\Psi}(E)$, then $\Pi_{k,E}^{\nabla, \Psi} w_h = w_h$ and $S^E(w_h, \cdot) = 0$, which proves (14).

As mentioned in the introduction, to circumvent the issues of discrete inverse/trace inequalities on non-polynomial spaces (involving the singular function), we adopt an analysis that is inspired by those developed for fully discrete methods, such as the Hybrid High-Order or Discrete De Rham methods. For such methods, two (equivalent) norms are classically introduced in the analysis: one based on all the degrees of freedom (internal and boundary), and one based on the available (higher-order) local polynomial reconstruction plus terms penalising the difference between this reconstruction and the boundary degrees of freedom; see, e.g., [20, Eqs. (2.7) and (2.15)] and [19, Sections 4.4 and 4.5]. Here, we bypass this need for two norms and only consider the latter one: (15) is a discrete $H^1$-norm built from the gradient of the local reconstruction $\Pi_{k,E}^{\nabla, \Psi} v_h$, plus penalisation terms between this reconstruction and the degrees of freedom $\Pi_{l,E}^{\Delta} v_h$ and $(v_h)|_{\partial E}$. The stabilisation (18) is then chosen to simply embed these penalisation terms, and as a consequence, for this particular choice, the inequalities in (17) are equalities.

To approximate the volumetric term on the right-hand side of (2), the forcing term $f$ is replaced with its orthogonal projection $f_h|_E = \Pi_{l,E}^{\Delta} f$ for all $E \in \Omega_h$, which shows that

$$ \int_{\Omega} f_h v_h \, dx := \sum_{E \in \Omega_h} \int_E (\Pi_{l,E}^{\Delta} f) v_h \, dx = \sum_{E \in \Omega_h} \int_E f (\Pi_{l,E}^{\Delta} v_h) \, dx \quad (19) $$
is computable from the DOFs of $v_h$ since $f$ is known. The X-VEM scheme then reads: Find $u_h \in V_{k,h,0}^{\Psi}$ such that

$$a_h(u_h, v_h) = \int_{\Omega} f_h v_h \quad \forall v_h \in V_{k,h,0}^{\Psi}. \qquad (20)$$

**Theorem 2.9 (Discrete Energy Error)** Let $u_h \in V_{k,h,0}^{\Psi}$ denote the solution to the X-VEM scheme (20) and $u = u^r + \psi \in H_0^1(\Omega)$, with $u^r \in H^{k+1}(\Omega_h)$ and $\psi \in \Psi(\Omega)$, the solution to the continuous problem (2). Under Assumption 2.1, the following energy error estimate holds:

$$\|u_h - \mathcal{I}_{k,h}u\|_{1,h} \lesssim h^k |u^r|_{H^{k+1}(\Omega_h)}, \qquad (21)$$

where $\mathcal{I}_{k,h}u = \hat{\mathcal{I}}_{k,h}u^r + \psi \in V_{k,h,0}^{\Psi}$ with $\hat{\mathcal{I}}_{k,h}: H^1(\Omega) \to V_k^h$ the standard VEM interpolant of $u^r$, see [2].

We note in passing that splitting the interpolant $\mathcal{I}_{k,h}$ into a singular component and an element of the regular VEM space is crucial to the analysis and relies on the fact that the space $V_{k,h}^{\Psi}(E)$ contains the regular VEM space.

**Remark 2.10 (Alternate interpolator and local enrichment)** A more natural way to define the interpolant on $V_{k,h}^{\Psi}$ is to take $\mathcal{I}_{k,h}v$ as the unique virtual function that has $((\Pi_{l,E}^{\Delta}v)_{E \in \Omega_h}, (\Pi_{k-2,e}^{\Delta}v)_{e \in \mathcal{E}_h}, (v(\nu))_{\nu \in \mathcal{V}_h})$ as degrees of freedom (D1)–(D3). With this definition, $\mathcal{I}_{k,h}$ leaves the enrichment space $\Psi$ invariant, and allows for local enrichment as briefly discussed in Remark 2.2.

For the upcoming analysis, we would need to establish optimal approximation properties of $\Pi_{k,E}^{\nabla,\Psi}\mathcal{I}_{k,h}$, which essentially requires to prove the boundedness of this mapping (defined on $H^2(E)$) in a scaled $H^2$-seminorm (see [20, Lemma 1.43 and Section 5.5.6.2]). Proving this boundedness is however not a trivial matter for the extended VEM, as discrete inverse and trace inequalities are not readily available in the non-polynomial spaces $\mathbb{P}_k^{\Psi}(E)$ and $\mathbb{P}_l^{\Delta}(E)$ (we note that it is already quite challenging for the regular VEM [12]). A more flexible approach to go in this direction would perhaps to use a fully discrete analysis (without direct usage of virtual functions, that are difficult to estimate), in the spirit of [33].

**Theorem 2.11 (H¹ Error)** Let $u_h \in V_{k,h,0}^{\Psi}$ and $u \in H_0^1(\Omega)$ be as in Theorem 2.9. Under Assumption 2.1, the following energy error estimate holds:

$$(\sum_{E \in \Omega_h} |\Pi_{k,E}^{\nabla,\Psi} u_h - u|_{H^1(E)}^2)^{\frac{1}{2}} \lesssim h^k |u^r|_{H^{k+1}(\Omega_h)}. \qquad (22)$$

# 3 Convergence analysis

We provide in this section a proof of Theorem 2.9. The proof hinges on the consistency of the scheme which we state in Theorem 3.5. We begin with the following lemma which guarantees the stability and well-posedness of the scheme, and is also crucial in validating Theorem 2.9. Given a mesh sequence satisfying Assumption 1.1, the following continuous trace inequality holds [13]:
**Lemma 3.1 (Continuous trace inequality)** For all $E \in \Omega_h$ and $v \in H^1(E)$, it holds,

$$h_E \|v\|_{\partial E}^2 \lesssim \|v\|_E^2 + h_E^2 \|\nabla v\|_E^2, \quad (23)$$

where the hidden constant depends only on the mesh regularity parameter $\rho$.

**Lemma 3.2** The mapping $\|\cdot\|_{1,h} : V_{k,h,0}^{\Psi} \to [0, \infty)$ defined by (16) is a norm.

*Proof.* As $\|\cdot\|_{1,h}$ is clearly a seminorm, we only have to prove that if $\|v_h\|_{1,h} = 0$ for some $v_h \in V_{k,h,0}^{\Psi}$, then $v_h = 0$. To this end, we note that $\|v_h\|_{1,h} = 0$ implies $\|\nabla \Pi_{k,E}^{\nabla, \Psi} v_h\|_E = 0$ for each $E \in \Omega_h$, and thus that $\Pi_{k,E}^{\nabla, \Psi} v_h$ is constant on each $E \in \Omega_h$. The condition $\|v_h - \Pi_{k,E}^{\nabla, \Psi} v_h\|_{\partial E} = 0$ shows that $v_h|_{\partial E} = (\Pi_{k,E}^{\nabla, \Psi} v_h)|_{\partial E}$ is also constant. Working from neighbour to neighbour and using the homogeneous condition $v_h|_{\partial \Omega} = 0$, we infer that those constants are all zero. Combining with the condition

$$\|\Pi_{l,E}^{\Delta}(v_h - \Pi_{k,E}^{\nabla, \Psi} v_h)\|_E = 0$$

we infer that for all $E \in \Omega_h$, $\Pi_{l,E}^{\Delta} v_h = 0$. Therefore all the DOFs of $v_h$ vanish, and thus $v_h = 0$. $\square$

**Lemma 3.3 (Consistency of $S^E$)** Let $S^E : V_{k,h}^{\Psi}(E) \times V_{k,h}^{\Psi}(E) \to \mathbb{R}$ be a stabilisation term satisfying equation (13). Then it holds for all $v = v^r + \psi \in H^1(E)$ with $v^r \in H^{k+1}(E)$ and $\psi \in \Psi(E)$ that

$$S^E(\mathcal{I}_{k,h}v, \mathcal{I}_{k,h}v) \lesssim [h_E^k |v^r|_{H^{k+1}(E)}]^2. \quad (24)$$

*Proof.* It follows from the definition of $\mathcal{I}_{k,h}$ and the $\mathbb{P}_k^{\Psi}(E)$-consistency (14) of $S^E$ together with the fact that $\psi + \Pi_{k,E}^{\nabla, \Psi} \hat{\mathcal{I}}_{k,h} v^r \in \mathbb{P}_k^{\Psi}(E)$ that

$$S^E(\mathcal{I}_{k,h}v, \mathcal{I}_{k,h}v) = S^E(\hat{\mathcal{I}}_{k,h}v^r + \psi, \hat{\mathcal{I}}_{k,h}v^r + \psi) = S^E(\hat{\mathcal{I}}_{k,h}v^r - \Pi_{k,E}^{\nabla, \Psi} \hat{\mathcal{I}}_{k,h}v^r, \hat{\mathcal{I}}_{k,h}v^r - \Pi_{k,E}^{\nabla, \Psi} \hat{\mathcal{I}}_{k,h}v^r).$$

Therefore, we infer from equation (13) that

$$S^E(\mathcal{I}_{k,h}v, \mathcal{I}_{k,h}v) \lesssim h_E^{-2} \|\Pi_{l,E}^{\Delta}(\hat{\mathcal{I}}_{k,h} - \Pi_{k,E}^{\nabla, \Psi} \hat{\mathcal{I}}_{k,h})v^r\|_E^2 + h_E^{-1} \|(\hat{\mathcal{I}}_{k,h} - \Pi_{k,E}^{\nabla, \Psi} \hat{\mathcal{I}}_{k,h})v^r\|_{\partial E}^2. \quad (25)$$

By the $L^2$-boundedness of $\Pi_{l,E}^{\Delta}$ and applying the continuous trace inequality (23) on the boundary term it follows that

$$h_E^{-2} \|\Pi_{l,E}^{\Delta}(\hat{\mathcal{I}}_{k,h} - \Pi_{k,E}^{\nabla, \Psi} \hat{\mathcal{I}}_{k,h})v^r\|_E^2 + h_E^{-1} \|(\hat{\mathcal{I}}_{k,h} - \Pi_{k,E}^{\nabla, \Psi} \hat{\mathcal{I}}_{k,h})v^r\|_{\partial E}^2 \\ \lesssim h_E^{-2} \|(\hat{\mathcal{I}}_{k,h} - \Pi_{k,E}^{\nabla, \Psi} \hat{\mathcal{I}}_{k,h})v^r\|_E^2 + \|\nabla(\hat{\mathcal{I}}_{k,h} - \Pi_{k,E}^{\nabla, \Psi} \hat{\mathcal{I}}_{k,h})v^r\|_E^2. \quad (26)$$

Substituting (26) into (25) and applying a Poincaré–Wirtinger inequality (due to the zero mean value of $(\hat{\mathcal{I}}_{k,h} - \Pi_{k,E}^{\nabla, \Psi} \hat{\mathcal{I}}_{k,h})v^r$ on $E$) yields

$$S^E(\mathcal{I}_{k,h}v, \mathcal{I}_{k,h}v) \lesssim \|\nabla(\hat{\mathcal{I}}_{k,h} - \Pi_{k,E}^{\nabla, \Psi} \hat{\mathcal{I}}_{k,h})v^r\|_E^2 \le \|\nabla(\hat{\mathcal{I}}_{k,h} - \Pi_{k,E}^{\nabla, \Psi} \hat{\mathcal{I}}_{k,h})v^r\|_E^2,$$
where the substitution of $\Pi_{k,E}^{\nabla, \Psi}$ with $\Pi_k^{\nabla, E}$ (the elliptic projector on $\mathbb{P}_k(E)$) is justified by the definition of $\Pi_{k,E}^{\nabla, \Psi}$, which ensures that $\nabla \Pi_{k,E}^{\nabla, \Psi} \hat{\mathcal{I}}_{k,h} v^r$ is the best $L^2$-approximation of $\nabla \hat{\mathcal{I}}_{k,h} v^r$ in $\nabla \mathbb{P}_k^\Psi(E) \supset \nabla \mathbb{P}_k(E)$, while $\nabla \Pi_k^{\nabla, E} \hat{\mathcal{I}}_{k,h} v^r \in \nabla \mathbb{P}_k(E)$. The proof of the consistency property (24) is complete by invoking the optimal approximation properties of $\Pi_k^{\nabla, E} \hat{\mathcal{I}}_{k,h}$ stated in [12, Lemma 2.23]. $\square$

**Lemma 3.4** It holds for all $v_h \in V_{k,h}^\Psi(E)$ that

$$ \|\nabla \Pi_{k,E}^{\nabla, \Psi} v_h\|_E + h_E^{-1} \|\Pi_{l,E}^{\Delta} v_h - \Pi_{k,E}^{\nabla, \Psi} v_h\|_E + h_E^{-\frac{1}{2}} \|v_h - \Pi_{k,E}^{\nabla, \Psi} v_h\|_{\partial E} + S^E(v_h, v_h)^{\frac{1}{2}} \lesssim \|v_h\|_{1,E}. \quad (27) $$

*Proof.* As each of $\|\nabla \Pi_{k,E}^{\nabla, \Psi} v_h\|_E$ and $h_E^{-\frac{1}{2}} \|v_h - \Pi_{k,E}^{\nabla, \Psi} v_h\|_{\partial E}$ appear in the definition of $\|\cdot\|_{1,E}$, their boundedness follows trivially. The boundedness of $S^E(v_h, v_h)^{\frac{1}{2}}$ is a direct result of (13) applied to $v_h - \Pi_{k,E}^{\nabla, \Psi} v_h$ and of the consistency (14) which gives $S^E(v_h - \Pi_{k,E}^{\nabla, \Psi} v_h, v_h - \Pi_{k,E}^{\nabla, \Psi} v_h) = S^E(v_h, v_h)$. It remains only to show that

$$ h_E^{-1} \|\Pi_{l,E}^{\Delta} v_h - \Pi_{k,E}^{\nabla, \Psi} v_h\|_E \lesssim \|v_h\|_{1,E}. $$

To see this, we add and subtract the term $\Pi_{l,E}^{\Delta} \Pi_{k,E}^{\nabla, \Psi} v_h$ and apply a triangle inequality to yield

$$ h_E^{-1} \|\Pi_{l,E}^{\Delta} v_h - \Pi_{k,E}^{\nabla, \Psi} v_h\|_E \le h_E^{-1} \|\Pi_{l,E}^{\Delta} v_h - \Pi_{l,E}^{\Delta} \Pi_{k,E}^{\nabla, \Psi} v_h\|_E + h_E^{-1} \|\Pi_{l,E}^{\Delta} \Pi_{k,E}^{\nabla, \Psi} v_h - \Pi_{k,E}^{\nabla, \Psi} v_h\|_E. $$

The bound

$$ h_E^{-1} \|\Pi_{l,E}^{\Delta} v_h - \Pi_{l,E}^{\Delta} \Pi_{k,E}^{\nabla, \Psi} v_h\|_E \le \|v_h\|_{1,E} $$

follows trivially from the definition of $\|\cdot\|_{1,E}$. As orthogonal projectors are the best approximations for their norm, and $\mathbb{P}_l(E) \subset \mathbb{P}_l^\Delta(E)$ it holds that

$$ h_E^{-1} \|\Pi_{l,E}^{\Delta} \Pi_{k,E}^{\nabla, \Psi} v_h - \Pi_{k,E}^{\nabla, \Psi} v_h\|_E \le h_E^{-1} \|\Pi_l^{0,E} \Pi_{k,E}^{\nabla, \Psi} v_h - \Pi_{k,E}^{\nabla, \Psi} v_h\|_E, $$

where $\Pi_l^{0,E}$ is the $L^2(E)$-orthogonal projector on $\mathbb{P}_l(E)$. Applying the approximation properties of $\Pi_l^{0,E}$ to yield

$$ h_E^{-1} \|\Pi_l^{0,E} \Pi_{k,E}^{\nabla, \Psi} v_h - \Pi_{k,E}^{\nabla, \Psi} v_h\|_E \lesssim \|\nabla \Pi_{k,E}^{\nabla, \Psi} v_h\|_E $$

completes the proof. $\square$

With $u$ the exact solution to (2), which we assume to satisfy the conditions in Theorem 2.9, the consistency error is given by the linear form $\mathcal{E}(u; \cdot) : V_{k,h,0}^\Psi \to \mathbb{R}$ defined for all $v_h \in V_{k,h,0}^\Psi$ as

$$ \mathcal{E}(u; v_h) := \int_{\Omega} f_h v_h \, d\mathbf{x} - a_h(\mathcal{I}_{k,h} u, v_h). \quad (28) $$

**Theorem 3.5 (Consistency error)** Recalling that the exact solution to (2) is written $u = u^r + \psi$, the consistency error satisfies the estimate

$$ |\mathcal{E}(u; v_h)| \lesssim h^k |u^r|_{H^{k+1}(\Omega_h)} \|v_h\|_{1,h}. \quad (29) $$
Proof. Applying the definition (12) of $a_h$ and the orthogonality properties of $\Pi_{k,E}^{\nabla, \Psi}$, it holds that

$$a_h(\mathcal{I}_{k,h}u, v_h) = \sum_{E \in \Omega_h} \left( \int_E \nabla \mathcal{I}_{k,h} u \cdot \nabla \Pi_{k,E}^{\nabla, \Psi} v_h \, d\mathbf{x} + S^E(\mathcal{I}_{k,h}u, v_h) \right). \quad (30)$$

Consider now, on each element $E \in \Omega_h$,

$$- \int_E \Delta u \Pi_{l,E}^{\Delta} v_h \, d\mathbf{x} = - \int_E \Delta u (\Pi_{l,E}^{\Delta} v_h - \Pi_{k,E}^{\nabla, \Psi} v_h) \, d\mathbf{x} + \int_E \nabla u \cdot \nabla \Pi_{k,E}^{\nabla, \Psi} v_h \, d\mathbf{x} \\ - \langle \nabla u \cdot \mathbf{n}, \Pi_{k,E}^{\nabla, \Psi} v_h \rangle_{-\frac{1}{2}, \frac{1}{2}, \partial E}, \quad (31)$$

which follows by adding and subtracting the term $\int_T \Delta u \Pi_{k,E}^{\nabla, \Psi} v_h \, d\mathbf{x}$ and integrating by parts. By the continuity of the virtual function $v_h$ and of the normal fluxes $\nabla u \cdot \mathbf{n}$ across the mesh edges and $v_h|_{\partial \Omega} = 0$, it holds that

$$0 = \langle \nabla u \cdot \mathbf{n}, v_h \rangle_{-\frac{1}{2}, \frac{1}{2}, \partial \Omega} = \sum_{E \in \Omega_h} \langle \nabla u \cdot \mathbf{n}_E, v_h \rangle_{-\frac{1}{2}, \frac{1}{2}, \partial E}. \quad (32)$$

Adding (32) to (31) and summing over all elements $E \in \Omega_h$, it follows that

$$- \sum_{E \in \Omega_h} \int_E \Delta u \Pi_{l,E}^{\Delta} v_h \, d\mathbf{x} = \sum_{E \in \Omega_h} \left[ - \int_E \Delta u (\Pi_{l,E}^{\Delta} v_h - \Pi_{k,E}^{\nabla, \Psi} v_h) \, d\mathbf{x} + \int_E \nabla u \cdot \nabla \Pi_{k,E}^{\nabla, \Psi} v_h \, d\mathbf{x} \right. \\ \left. \qquad + \langle \nabla u \cdot \mathbf{n}, v_h - \Pi_{k,E}^{\nabla, \Psi} v_h \rangle_{-\frac{1}{2}, \frac{1}{2}, \partial E} \right]. \quad (33)$$

Let $z \in \mathbb{P}_k^{\Psi}(E)$. Since $\Delta z \in \Delta \mathbb{P}_k^{\Psi}(E) \subset \mathbb{P}_l^{\Delta}(E)$ and $\Pi_{l,E}^{\Delta}$ is the $L^2$-projector on that space, we have

$$\int_E \Delta z (\Pi_{l,E}^{\Delta} v_h - \Pi_{k,E}^{\nabla, \Psi} v_h) \, d\mathbf{x} = \int_E \Delta z (v_h - \Pi_{k,E}^{\nabla, \Psi} v_h) \, d\mathbf{x}.$$

Hence, integrating by part and using the definition (9) of $\Pi_{k,E}^{\nabla, \Psi}$,

$$- \int_E \Delta z (\Pi_{l,E}^{\Delta} v_h - \Pi_{k,E}^{\nabla, \Psi} v_h) \, d\mathbf{x} + \langle \nabla z \cdot \mathbf{n}, v_h - \Pi_{k,E}^{\nabla, \Psi} v_h \rangle_{-\frac{1}{2}, \frac{1}{2}, \partial E} = \int_E \nabla z \cdot \nabla (v_h - \Pi_{k,E}^{\nabla, \Psi} v_h) \, d\mathbf{x} = 0.$$

Applying this to $z = \psi + \Pi_k^{\nabla, E} u^r \in \mathbb{P}_k^{\Psi}(E)$, we can cancel in the following expression the singular component $\psi$ of $u$ and introduce the term $\Pi_k^{\nabla, E} u^r$ to write

$$- \int_E \Delta u (\Pi_{l,E}^{\Delta} v_h - \Pi_{k,E}^{\nabla, \Psi} v_h) \, d\mathbf{x} + \langle \nabla u \cdot \mathbf{n}, v_h - \Pi_{k,E}^{\nabla, \Psi} v_h \rangle_{-\frac{1}{2}, \frac{1}{2}, \partial E} \\ = - \int_E \Delta (u^r - \Pi_k^{\nabla, E} u^r) (\Pi_{l,E}^{\Delta} v_h - \Pi_{k,E}^{\nabla, \Psi} v_h) \, d\mathbf{x} + \langle \nabla (u^r - \Pi_k^{\nabla, E} u^r) \cdot \mathbf{n}, v_h - \Pi_{k,E}^{\nabla, \Psi} v_h \rangle_{-\frac{1}{2}, \frac{1}{2}, \partial E}.$$
Substituting back into (33) yields

$$ - \sum_{E \in \Omega_h} \int_E \Delta u \Pi_{l,E}^{\Delta} v_h \, d\boldsymbol{x} = \sum_{E \in \Omega_h} \left[ \int_E \nabla u \cdot \nabla \Pi_{k,E}^{\nabla, \Psi} v_h \, d\boldsymbol{x} - \int_E \Delta(u^r - \Pi_k^{\nabla, E} u^r)(\Pi_{l,E}^{\Delta} v_h - \Pi_{k,E}^{\nabla, \Psi} v_h) \, d\boldsymbol{x} + \int_{\partial E} \nabla(u^r - \Pi_k^{\nabla, E} u^r) \cdot \boldsymbol{n} (v_h - \Pi_{k,E}^{\nabla, \Psi} v_h) \, ds \right], \quad (34) $$

where we have replaced the $\langle \cdot, \cdot \rangle_{-\frac{1}{2}, \frac{1}{2}, \partial E}$ duality product with the integral notation due to the regularity of $u^r$. Recalling the definition (28) of $\mathcal{E}(u; v_h)$ and invoking (19) together with $f = -\Delta u$, we infer

$$ \begin{aligned} \mathcal{E}(u; v_h) &= - \sum_{E \in \Omega_h} \int_E \Delta u \Pi_{l,E}^{\Delta} v_h \, d\boldsymbol{x} - a_h(\mathcal{I}_{k,h} u, v_h) \\ &= \sum_{E \in \Omega_h} \left[ \int_E \nabla(u - \mathcal{I}_{k,h} u) \cdot \nabla \Pi_{k,E}^{\nabla, \Psi} v_h \, d\boldsymbol{x} - S^E(\mathcal{I}_{k,h} u, v_h) \right. \\ & \qquad \left. - \int_E \Delta(u^r - \Pi_k^{\nabla, E} u^r)(\Pi_{l,E}^{\Delta} v_h - \Pi_{k,E}^{\nabla, \Psi} v_h) \, d\boldsymbol{x} + \int_{\partial E} \nabla(u^r - \Pi_k^{\nabla, E} u^r) \cdot \boldsymbol{n} (v_h - \Pi_{k,E}^{\nabla, \Psi} v_h) \, ds \right], \end{aligned} $$

where the second equality follows from (34) and (30). By applying Cauchy-Schwarz inequalities to each of the terms and (27) we obtain

$$ \mathcal{E}(u; v_h) \lesssim \sum_{E \in \Omega_h} \|v_h\|_{1,E} \left[ \|\nabla(u - \mathcal{I}_{k,h} u)\|_E + S^E(\mathcal{I}_{k,h} u, \mathcal{I}_{k,h} u)^{\frac{1}{2}} \right. \\ \left. \qquad + h_E \|\Delta(u^r - \Pi_k^{\nabla, E} u^r)\|_E + h_E^{\frac{1}{2}} \|\nabla(u^r - \Pi_k^{\nabla, E} u^r) \cdot \boldsymbol{n}\|_{\partial E} \right]. $$

It follows from the definition of $\mathcal{I}_{k,h}$ that $u - \mathcal{I}_{k,h} u = u^r - \hat{\mathcal{I}}_{k,h} u^r$ and thus applying the approximation properties of $\hat{\mathcal{I}}_{k,h}$ [12, Lemma 2.23] yields the bound

$$ \|\nabla(u - \mathcal{I}_{k,h} u)\|_E \lesssim h_E^k |u^r|_{H^{k+1}(E)}. $$

Combining with the consistency (24) of $S^E$ and the approximation properties of $\Pi_k^{\nabla, E}$, it is clear that

$$ \mathcal{E}(u; v_h) \lesssim \sum_{E \in \Omega_h} \|v_h\|_{1,E} h_E^k |u^r|_{H^{k+1}(E)}. $$

The result follows from a Cauchy-Schwarz inequality and the bound $h_E \le h$. $\square$

The proof of Theorem 2.9 now follows trivially.

*Proof of Theorem 2.9.* It follows from the coercivity property (17) and the Third Strang Lemma [20, Lemma A.7] that

$$ \|u_h - \mathcal{I}_{k,h} u\|_{1,h} \lesssim \sup_{\|v_h\|_{1,h}=1} |\mathcal{E}(u; v_h)|. $$
The proof then follows from the residual error (29). □

*Proof of Theorem 2.11.* It follows from introducing $\Pi_{k,E}^{\nabla,\Psi}u$ and $\Pi_{k,E}^{\nabla,\Psi}\mathcal{I}_{k,h}u$ on each element $E \in \Omega_h$ and applying a triangle inequality that

$$ \sum_{E \in \Omega_h} |\Pi_{k,E}^{\nabla,\Psi} u_h - u|_{H^1(E)}^2 \\ \lesssim \sum_{E \in \Omega_h} |\Pi_{k,E}^{\nabla,\Psi} (u_h - \mathcal{I}_{k,h}u)|_{H^1(E)}^2 + \sum_{E \in \Omega_h} |\Pi_{k,E}^{\nabla,\Psi} (\mathcal{I}_{k,h}u - u)|_{H^1(E)}^2 + \sum_{E \in \Omega_h} |\Pi_{k,E}^{\nabla,\Psi} u - u|_{H^1(E)}^2. \quad (35) $$

By definition of $\mathcal{I}_{k,h}$, we have $\mathcal{I}_{k,h}u - u = \hat{\mathcal{I}}_{k,h}u^r - u^r$ and thus a bound on the second term in the right-hand side of (35) follows from the $H^1(E)$-boundedness of $\Pi_{k,E}^{\nabla,\Psi}$ and the approximation properties [12, Lemma 2.23] of $\hat{\mathcal{I}}_{k,h}$. After noticing that $\Pi_{k,E}^{\nabla,\Psi}u - u = \Pi_{k,E}^{\nabla,\Psi}u^r - u^r$ since $\Pi_{k,E}^{\nabla,\Psi}\psi = \psi$, the last term is bounded due to the optimal approximation properties of $\Pi_{k,E}^{\nabla,\Psi}$ in the $H^1(E)$-seminorm. To bound the first term, note that

$$ \sum_{E \in \Omega_h} |\Pi_{k,E}^{\nabla,\Psi} (u_h - \mathcal{I}_{k,h}u)|_{H^1(E)}^2 \lesssim \|u_h - \mathcal{I}_{k,h}u\|_{1,h}^2 $$

and the result follows from Theorem 2.9. □

# 4 Numerical experiments

The XVEM method described in Section 2 is implemented using the open-source C++ library PolyMesh (https://github.com/liamyemm/PolyMesh), using linear algebra tools from the Eigen3 library (see http://eigen.tuxfamily.org). We focus here on testing the method on domains possessing fractures or re-entrant corners. For both these cases, the exact solution is not expected to be $H^2$, and so a classical virtual element method is expected to converge sub-optimally for all $k$.

Enriching elements and edges with singular functions far from the location of the singularity can cause severe ill-conditioning due to these functions being well-approximated by polynomials. This is typically mitigated through local enrichment [5, 33]. However, the analysis of the method described in this paper requires splitting the interpolant into a singular part, and an classical VEM interpolant on the regular part (see Remarks 2.2 and 2.10). This requires the enrichment function to be globally $H^1$-conforming. As such, a sufficiently smooth cut-off function would be required to facilitate local enrichment.

However, numerical tests suggest that removing the additional DOFs far from the singularity still leads to a valid scheme. In particular, we can enrich the local spaces on the elements and their edges only if the element intersects the disk of radius $\gamma > 0$ centred at the singularity. For all tests the local enrichment parameter is taken as $\gamma = 0.15$. While such an enrichment is not accounted for by the analysis in this paper, the tests show that locally enriching in this manner achieves optimal results.
## 4.1 Fractured domain

Consider the fractured domain $\Omega = (-1, 1)^2 \setminus \{(x, y) : y = 0, x > 0\}$. The solution to problem (1) posed in this domain is expected to contain a singularity at the fracture tip of the form [24]

$$ \psi(r, \theta) = r^{\frac{1}{2}} \sin(\frac{1}{2}\theta). $$

As such, we test with an exact solution

$$ u = \sin(\pi x) \sin(\pi y) + \psi $$

and enrich the local spaces with the function $\psi$.

A sequence of Cartesian meshes of the domain $\Omega$ is considered. The parameters of these meshes are displayed in Table 1 and two meshes are plotted in Figure 1 showing the fracture and the local enrichment scheme.

<table><caption>Table 1: Mesh data used for the fractured domain test</caption><thead><tr><th>Mesh #</th><th>h</th><th>Nb. Elements</th><th>Nb. Edges</th><th>Nb. Vertices</th></tr></thead><tbody><tr><td>1</td><td>0.7071</td><td>16</td><td>40</td><td>25</td></tr><tr><td>2</td><td>0.3536</td><td>64</td><td>144</td><td>81</td></tr><tr><td>3</td><td>0.1768</td><td>256</td><td>544</td><td>289</td></tr><tr><td>4</td><td>0.0884</td><td>1,024</td><td>2,112</td><td>1,089</td></tr><tr><td>5</td><td>0.0442</td><td>4,096</td><td>8,320</td><td>4,225</td></tr><tr><td>6</td><td>0.0221</td><td>16,384</td><td>33,024</td><td>16,641</td></tr></tbody></table>

The solution to the discrete problem (20) is denoted by $u_h$ and the solution to the continuous problem (1) is denoted by $u$. The accuracy of the scheme is determined by the following relative errors (respectively measuring the error in an $L^2$-like norm and an $H^1$-like seminorm):

$$ E_{0,h}^2 := \frac{\sum_{E \in \Omega_h} \|\Pi_{k,E}^{\nabla, \Psi}(u_h - u)\|_E^2}{\sum_{E \in \Omega_h} \|\Pi_{k,E}^{\nabla, \Psi} u\|_E^2} \quad \text{and} \quad E_{1,h}^2 := \frac{\sum_{E \in \Omega_h} |\Pi_{k,E}^{\nabla, \Psi}(u_h - u)|_{H^1(E)}^2}{\sum_{E \in \Omega_h} |\Pi_{k,E}^{\nabla, \Psi} u|_{H^1(E)}^2}. $$

By Theorem 2.11, for the X-VEM the error $E_{1,h}$ is expected to decay as $O(h^k)$ despite the presence of the singularity in $u$; we have not studied the convergence in $L^2$-norm but, classically, we would expect $E_{0,h}$ to decay as $O(h^{k+1})$. On the contrary, for the non-enriched VEM these errors will be strongly limited by the lack of regularity of $u$.

As enriching the virtual element method increases the number of degrees of freedom, it is most appropriate to plot the error against the degrees of freedom. After static condensation is performed, and the boundary degrees of freedom are fixed, the only remaining degrees of freedom are those corresponding to the internal edges and vertices. As the number of DOFs grows linearly with the number of elements, which is of order $\sim h^{-2}$ for our meshes, the $O(h^k)$ error estimate predicted for the error in $(H^1$-like) energy norm by (21) translates to $O((\#DOFs)^{-k/2})$ in our graphs.
<figure><img src="image_1.png" /><figcaption>Figure 1: Plots of meshes 2 and 3 used in the fractured domain test. The fracture is coloured blue and enriched elements are coloured red.</figcaption></figure>

In Figure 2, the error of the scheme is plotted against these remaining DOFs for a non-enriched, globally enriched and locally enriched scheme. The ill-conditioning of the globally enriched method is apparent due to the scheme failing on Meshes 5 & 6 when $k = 2, 3$. Moreover, for a given mesh and polynomial degree $k$, the globally enriched scheme typically has a worse $H^1$ error than the locally enriched scheme, suggesting that ill-conditioning is plaguing the scheme. When $k = 1, 2$, the non-enriched scheme converges optimally in $L^2$ error, however the convergence rate is sub-optimal in $H^1$ error. For $k = 2$ the non-enriched scheme converges sub-optimally in both $L^2$ and $H^1$ error and is significantly outperformed by both enriched schemes. The locally enriched scheme converges optimally in all cases and is clearly the best performing scheme.

**Remark 4.1 (Rate of convergence in $L^2$ norm)** For $k = 1, 3$ we note an optimal convergence rate in $L^2$-norm, with an error decaying as $O((\|\mathrm{DOFs}\|)^{-(k+1)/2})$. The rate of convergence in $L^2$-norm for $k = 2$ however appears to be sub-optimal, and not better than the rate for the $H^1$-norm. Given our choice of (depleted) element unknowns in the VEM space, the VEM scheme for $k = 2$ can be compared with the (depleted) HHO scheme for $(k, \ell) = (1, 0)$ [20, Chapter 5.1], for which the sub-optimal convergence in $L^2$-norm is a well documented phenomenon (see [20, Remark 5.17 and Section 5.1.8]).

**Remark 4.2 (Managing ill-conditioning)** As discussed above, local enrichment leads to algebraic systems (both in computations of local operators, and for the final fully coupled system) that are much better conditioned than those obtained via global enrichment.

Another avenue to mitigate the ill-conditioning of the system obtained via global enrichment would be through the use of a pre-conditioner. Our current implementation solves the global linear system through a direct application of the BiCGStab solver from the Eigen3 library; it is expected that classical pre-conditioners, such as the BDDC method (see [11] for the application of this method to standard VEM), could lead to better numerical outcomes for the globally enriched scheme.
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /><figcaption>(a) $E_{0,h}$ vs DOFs, $k=1$</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>(b) $E_{1,h}$ vs DOFs, $k=1$</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>(c) $E_{0,h}$ vs DOFs, $k=2$</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>(d) $E_{1,h}$ vs DOFs, $k=2$</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(e) $E_{0,h}$ vs DOFs, $k=3$</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(f) $E_{1,h}$ vs DOFs, $k=3$</figcaption></figure>

<caption>Figure 2: Tests on fractured domain</caption>
## 4.2 L-shaped domain, hexagonal meshes

Consider the fractured domain $\Omega = (-1, 1)^2 \setminus [0, 1)^2$. The solution in this domain has a singularity located at the re-entrant corner of the form

$$ \psi(r, \theta) = r^{\frac{2}{3}} \sin(\frac{2}{3}(\theta - \frac{\pi}{2})). $$

Analogous to the fractured domain test, we consider an exact solution of the form

$$ u = \sin(\pi x) \sin(\pi y) + \psi $$

and enrich the local spaces with the function $\psi$.

A sequence of hexagonal meshes of the domain $\Omega$ is considered. The parameters of these meshes are displayed in Table 2 and two meshes are plotted in Figure 3 showing the enrichment scheme considered for the locally enriched tests.

<table><caption>Table 2: Mesh data used for the L-shaped domain test</caption><thead><tr><th>Mesh #</th><th>h</th><th>Nb. Elements</th><th>Nb. Edges</th><th>Nb. Vertices</th></tr></thead><tbody><tr><td>1</td><td>0.7454</td><td>21</td><td>76</td><td>56</td></tr><tr><td>2</td><td>0.4177</td><td>65</td><td>224</td><td>160</td></tr><tr><td>3</td><td>0.2992</td><td>133</td><td>444</td><td>312</td></tr><tr><td>4</td><td>0.1665</td><td>481</td><td>1,536</td><td>1,056</td></tr><tr><td>5</td><td>0.1019</td><td>1,281</td><td>4,000</td><td>2,720</td></tr><tr><td>6</td><td>0.0514</td><td>4,961</td><td>15,200</td><td>10,240</td></tr><tr><td>7</td><td>0.0257</td><td>19,521</td><td>59,200</td><td>39,680</td></tr></tbody></table>

Figure 4 illustrates the error of the method plotted against the remaining internal DOFs after static condensation for a non-enriched, globally enriched, and locally enriched scheme. As with the fractured domain, the globally enriched approach exhibits evident ill-conditioning, as observed by its failure on Mesh 7 when $k = 1$, Mesh 6 & 7 when $k = 2$ and on Mesh 4, 5, 6 & 7 when $k = 3$. Indeed, for $k = 1, 2$, the error is greater for the globally enriched scheme than for either the non-enriched or locally enriched schemes. As with the fractured domain test, the locally enriched scheme is the best performing and converges optimally in all cases.

## 4.3 L-shaped domain, Cartesian meshes

In this section, we consider the test case of [5, Section 5.1.1]. The domain is still an L-shaped one, albeit with the bottom right corner removed: $\Omega = (-1, 1)^2 \setminus [0, 1) \times (-1, 0]$. The exact solution and singular function are

$$ u = \sin(\pi x) \sin(\pi y) + \psi \quad \text{and} \quad \psi(r, \theta) = r^{\frac{2}{3}} \sin(\frac{2}{3}\theta). $$

We run the simulations using uniform Cartesian meshes (see their characteristics in Table 3). The convergence graphs for the globally, locally and non-enriched schemes are presented in Figure 5; note that there
<figure><img src="image_9.png" /><figcaption>Figure 3: Plots of meshes 2 and 4 used in the L-shaped domain test of Section 4.2. The enriched elements are coloured red.</figcaption></figure>

is no data for the non-enriched scheme with $k = 1$ on the coarsest mesh since this mesh only has boundary vertices and, as a consequence, there are zero DOFs for that scheme on this mesh.

As in the previous tests, the locally enriched version outperforms the other two, with the non-enriched version displaying a sub-optimal convergence rate while the globally enriched version shows optimal convergence but larger magnitudes of error. The initial error growths in Figure 5 are probably due to the use in the denominator of $E_{1,h}$ of the mesh-dependent quantity $(\sum_E |\Pi_{k,E}^\nabla, \Psi u|_{H^1}^2)^{1/2}$, which fluctuates a lot on the first coarse meshes before stabilising to a number close to $|u|_{H^1(\Omega)}$ on finer meshes. The ill-conditioning issues of X-VEM mentioned above are visible in Figure 5 in the sense that the linear solver fails to converge on the finest mesh for the globally enriched scheme with degrees $k = 2, 3$; here too, we notice that local enrichment leads to a better conditioned scheme.

These results are comparable to the ones obtained with the enriched NCVEM in [5, Fig. 11]. A direct quantitative comparison is however not feasible, because the errors in VEM and NCVEM are measured in different norms (through the use of scheme-specific reconstruction operators). The ill-conditioning of the globally enriched NCVEM is visible in [5, Fig. 11] as a blow-up of the error. A mitigation was proposed there using new degrees of freedom corresponding to orthonormal polynomial basis functions on the edges. We also use such edge basis functions, but there is a fundamental difference between conforming VEM and NCVEM (and HHO), which could explain the worse ill-conditioning encountered in the former method: in NCVEM/HHO, the edge degrees of freedom directly give the boundary contributions in the local equations solved to compute the elliptic projectors/potential reconstructions (see [5, Eqs. (68) and following] or [33, Eq. (2.8)]), so the orthonormalisation of the bases corresponding to these degrees of freedom has a direct beneficial aspect in the computation of this boundary contribution; on the contrary, in VEM, the boundary contributions to the elliptic projectors involves the whole boundary value of the virtual function (see (11)), which do not directly correspond to edge degrees of freedom but has to be reconstructed from vertex and edge degrees of freedom; this reconstruction probably reduces the positive impact of having chosen
<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /><figcaption>(a) $E_{0,h}$ vs DOFs, $k=1$</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(b) $E_{1,h}$ vs DOFs, $k=1$</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(c) $E_{0,h}$ vs DOFs, $k=2$</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(d) $E_{1,h}$ vs DOFs, $k=2$</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(e) $E_{0,h}$ vs DOFs, $k=3$</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(f) $E_{1,h}$ vs DOFs, $k=3$</figcaption></figure>

<caption>Figure 4: Tests on L-shaped domain, hexagonal meshes</caption>

orthonormal bases for polynomial edge spaces.
<table><caption>Table 3: Uniform Cartesian meshes for the L-shaped domain</caption><thead><tr><th><i>h</i></th><th>Nb. Elements</th><th>Nb. Edges</th><th>Nb. Vertices</th></tr></thead><tbody><tr><td>1.4142</td><td>3</td><td>10</td><td>8</td></tr><tr><td>0.7071</td><td>12</td><td>32</td><td>21</td></tr><tr><td>0.3536</td><td>48</td><td>112</td><td>65</td></tr><tr><td>0.1768</td><td>192</td><td>416</td><td>225</td></tr><tr><td>0.0884</td><td>768</td><td>1,600</td><td>833</td></tr></tbody></table>

<figure><img src="image_17.png" /><figcaption>Figure 5: Tests on L-shaped domain, Cartesian meshes</figcaption></figure>

# 5 Conclusions

In this paper, a novel extension of the virtual element method has been introduced, which achieves consistency on highly generic enrichment spaces. This is applicable to problems posed in complex geometries, where singularities are known to exist near corners and fractures. However, the method also has applicability
to any scenario where some component of the solution can be asymptotically obtained. This could include, for example, singular behaviour in the source term, or a highly oscillatory component of the solution. While the method was formulated for the Poisson problem in two dimensions, it has natural extensions to higher dimensions and more general linear elliptic problems.

We performed a complete analysis of the proposed method, proving its capacity to achieve optimal convergence rate in discrete and continuous energy norm. Numerical tests support the theoretical findings in the presence of corner and fracture singularities in both $L^2$- and $H^1$-norm. However, ill-conditioning issues plague the globally enriched scheme – something also observed by [5, 33]. The analysis in this paper does not account for the method of local enrichment considered in Section 4, however, numerical tests show optimal convergence. Providing a robust analysis for the X-VEM with discontinuous enrichment spaces is a potential avenue of future research.

## Acknowledgement

The work of JD was partially supported by the Australian Government through the Australian Research Council's Discovery Projects funding scheme (project number DP210103092) and by the European Union (ERC Synergy, NEMESIS, project number 101115663); views and opinions expressed are however those of the authors only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency. Neither the European Union nor the granting authority can be held responsible for them.

The work of GM was partially supported by the Laboratory Directed Research and Development (LDRD) program of Los Alamos National Laboratory under project number 20220129ER. Los Alamos National Laboratory is operated by Triad National Security, LLC, for the National Nuclear Security Administration of U.S. Department of Energy (Contract No. 89233218CNA000001). GM is a member of the Gruppo Nazionale Calcolo Scientifico-Istituto Nazionale di Alta Matematica (GNCS-INdAM).

The work of LY was supported by a Monash University Postgraduate Publications Award.

## References

[1] R. A. Adams and J. J. F. Fournier. Sobolev spaces. 2nd ed. Pure and Applied Mathematics. Academic Press, 2003.

[2] B. Ahmad, A. Alsaedi, F. Brezzi, L. D. Marini, and A. Russo. "Equivalent projectors for virtual element methods". In: *Computers & Mathematics with Applications* 66 (2013), pp. 376–391.

[3] P. F. Antonietti, F. Brezzi, and L. D. Marini. "Bubble stabilization of discontinuous Galerkin methods". In: *Comput. Methods Appl. Mech. Engrg.* 198.21-26 (2009), pp. 1651–1659. DOI: 10.1016/j.cma.2008.12.033.

[4] E. Artioli and L. Mascotto. "Enriched virtual elements for plane elasticity with corner singularities". In: *Comput. Mech.* (2023). DOI: 10.1007/s00466-023-02418-4.

[5] E. Artioli and L. Mascotto. "Enrichment of the nonconforming virtual element method with singular functions". In: *Computer Methods in Applied Mechanics and Engineering* 385 (2021), p. 114024.
[6] I. Babuška and J. M. Melenk. "The partition of unity method". In: *International journal for numerical methods in engineering* 40.4 (1997), pp. 727–758.

[7] S. Badia, J. Droniou, and L. Yemm. "Conditioning of a hybrid high-order scheme on meshes with small faces". In: *J. Sci. Comput.* 92.2 (2022), Paper No. 71, 23. ISSN: 0885-7474. DOI: 10.1007/s10915-022-01913-9. URL: https://arxiv.org/abs/2109.09983.

[8] T. Belytschko and T. Black. "Elastic crack growth in finite elements with minimal remeshing". In: *International journal for numerical methods in engineering* 45.5 (1999), pp. 601–620.

[9] E. Benvenuti, A. Chiozzi, G. Manzini, and N. Sukumar. "Extended virtual element method for the Laplace problem with singularities and discontinuities". In: *Computer Methods in Applied Mechanics and Engineering* 356 (2019), pp. 571–597.

[10] E. Benvenuti, A. Chiozzi, G. Manzini, and N. Sukumar. "Extended virtual element method for two-dimensional linear elastic fracture". In: *Computer Methods in Applied Mechanics and Engineering* 390 (2022), p. 114352.

[11] S. Bertoluzza, M. Pennacchio, and D. Prada. "BDDC and FETI-DP for the virtual element method". In: *Calcolo* 54 (2017), pp. 1565–1593.

[12] S. C. Brenner, Q. Guan, and L.-Y. Sung. "Some estimates for virtual element methods". In: *Computational Methods in Applied Mathematics* 17.4 (2017), pp. 553–574.

[13] S. C. Brenner and L.-Y. Sung. "Virtual element methods on meshes with small edges or faces". In: *Mathematical Models and Methods in Applied Sciences* 28.07 (2018), pp. 1291–1336.

[14] A. Cangiani, Z. Dong, E. H. Georgoulis, and P. Houston. "hp-version discontinuous Galerkin methods for advection-diffusion-reaction problems on polytopic meshes". In: *ESAIM Math. Model. Numer. Anal.* 50.3 (2016), pp. 699–725. DOI: 10.1051/m2an/2015059.

[15] B. Cockburn, D. A. Di Pietro, and A. Ern. "Bridging the Hybrid High-Order and Hybridizable Discontinuous Galerkin methods". In: *ESAIM: Math. Model. Numer. Anal.* 50.3 (2016), pp. 635–650. DOI: 10.1051/m2an/2015051.

[16] B. Cockburn, J. Gopalakrishnan, and R. Lazarov. "Unified hybridization of discontinuous Galerkin, mixed, and continuous Galerkin methods for second order elliptic problems". In: *SIAM J. Numer. Anal.* 47.2 (2009), pp. 1319–1365. DOI: 10.1137/070706616.

[17] D. A. Di Pietro and A. Ern. *Mathematical aspects of discontinuous Galerkin methods*. Vol. 69. Mathématiques & Applications (Berlin) [Mathematics & Applications]. Springer, Heidelberg, 2012, pp. xviii+384. ISBN: 978-3-642-22979-4. DOI: 10.1007/978-3-642-22980-0.

[18] D. A. Di Pietro, A. Ern, and S. Lemaire. "An arbitrary-order and compact-stencil discretization of diffusion on general meshes based on local reconstruction operators". In: *Comput. Meth. Appl. Math.* 14.4 (2014), pp. 461–472. DOI: 10.1515/cmam-2014-0018.

[19] D. A. Di Pietro and J. Droniou. "An arbitrary-order discrete de Rham complex on polyhedral meshes: Exactness, Poincaré inequalities, and consistency". In: *Found. Comput. Math.* 23 (2023), pp. 85–164. DOI: 10.1007/s10208-021-09542-8.
[20] D. A. Di Pietro and J. Droniou. *The Hybrid High-Order Method for Polytopal Meshes: Design, Analysis, and Applications*. Vol. 19. Modeling, Simulation and Applications. Springer International Publishing, 2020, xxxi + 525p. ISBN: 978-3-030-37202-6 (Hardcover), 978-3-030-37203-3 (eBook). DOI: 10.1007/978-3-030-37203-3. URL: https://hal.archives-ouvertes.fr/hal-02151813.

[21] J. Droniou. "Finite volume schemes for diffusion equations: introduction to and review of modern methods". In: *Math. Models Methods Appl. Sci.* 24.8 (2014), pp. 1575–1619. DOI: 10.1142/S0218202514400041.

[22] J. Droniou, R. Eymard, T. Gallouët, and R. Herbin. "A unified approach to mimetic finite difference, hybrid finite volume and mixed finite volume methods". In: *Math. Models Methods Appl. Sci.* (M3AS) 20.2 (2010), pp. 1–31. DOI: 10.1142/S0218202510004222.

[23] R. Eymard, T. Gallouët, and R. Herbin. "Discretization of heterogeneous and anisotropic diffusion problems on general nonconforming meshes. SUSHI: a scheme using stabilization and hybrid interfaces". In: *IMA J. Numer. Anal.* 30.4 (2010), pp. 1009–1043. DOI: 10.1093/imanum/drn084.

[24] P. Grisvard. *Elliptic Problems in Nonsmooth Domains*. Boston, MA: Pitman Publishing, Inc, 1985.

[25] Y. Kuznetsov, K. Lipnikov, and M. Shashkov. "Mimetic finite difference method on polygonal meshes for diffusion-type problems". In: *Comput. Geosci.* 8 (2004), pp. 301–324.

[26] J. M. Melenk and I. Babuška. "The partition of unity finite element method: basic theory and applications". In: *Computer methods in applied mechanics and engineering* 139.1-4 (1996), pp. 289–314.

[27] N. Moës, J. Dolbow, and T. Belytschko. "A finite element method for crack growth without remeshing". In: *International journal for numerical methods in engineering* 46.1 (1999), pp. 131–150.

[28] M. Motz. "The treatment of singularities of partial differential equations by relaxation methods". In: *Quarterly of Applied Mathematics* 4 (1947), pp. 371–377.

[29] L. Mu, J. Wang, and X. Ye. "Weak Galerkin finite element methods on polytopal meshes". In: *Int. J. Numer. Anal. Model*. 12.1 (2015), pp. 31–53.

[30] L. Beirão da Veiga, F. Brezzi, A. Cangiani, G. Manzini, L. D. Marini, and A. Russo. "Basic principles of virtual element methods". In: *Mathematical Models & Methods in Applied Sciences* 23 (2013), pp. 119–214.

[31] L. Beirão da Veiga, K. Lipnikov, and G. Manzini. *The mimetic finite difference method for elliptic problems*. Vol. 11. MS&A. Modeling, Simulation and Applications. Springer, Cham, 2014, pp. xvi+392. ISBN: 978-3-319-02662-6; 978-3-319-02663-3. DOI: 10.1007/978-3-319-02663-3.

[32] L. Yemm. "An enriched hybrid high-order method for the Stokes problem with application to flow around submerged cylinders". In: (2024). URL: https://arxiv.org/abs/2311.13243.

[33] L. Yemm. "Design and analysis of the extended hybrid high-order method for the Poisson problem". In: *Advances in Computational Mathematics* 48.4 (2022), pp. 1–25. DOI: 10.1007/s10444-022-09958-y.