# A new paradigm for enriching virtual element methods

E. Artioli, L. Mascotto Dipartimento di Ingegneria Civile e Ingegneria
Informatica, Università di Roma Tor Vergata, 00133 Rome, Italy
(artioli@ing.uniroma2.it)Fakultät für Mathematik, Universität Wien, 1090
Vienna, Austria<EMAIL_ADDRESS>
###### Abstract

We construct a virtual element method (VEM) based on approximation spaces that
are enriched with special singular functions. This enriched VEM is tailored
for the approximation of solutions to elliptic problems, which have
singularities due to the geometry of the domain. Differently from the
traditional extended Galerkin method approach, based on the enrichment of
local spaces with singular functions, no partition of unity is employed.
Rather, the design of the method hinges upon the special structure of the
virtual element spaces. We discuss the full theoretical analysis of the method
and support it with several numerical experiments. We also present an
orthonormalization procedure drastically trimming the ill-conditioning of the
final system.

AMS subject classification: 65N12, 65N15, 65N30

Keywords: virtual element methods, extended Galerkin methods, enrichment with
singular functions, optimal convergence, polygonal meshes

##  1 Introduction

The virtual element method (VEM) is a recent generalization of the finite
element method (FEM) to very general polygonal/polyhedral meshes; see [12].
Differently from several other polytopal methods, virtual element spaces are
designed to mimic properties of the solution to the problem under
consideration. This is very much in the spirit of Trefftz methods and renders
the VEM extremely similar to the boundary element method-based FEM [42]. Even
on standard triangular and tetrahedral meshes, new elements can be
constructed.

The design of special virtual element spaces mimicking the continuous problem
has been exploited in various occasions. Amongst them, we mention the
approximation of solutions to the Stokes equation [17] with divergence free
spaces; polyharmonic problems [25, 1] with polyharmonic virtual element
spaces; problems with zero right-hand side tackled with the Trefftz VEM [34];
elasticity problems [3, 28] with symmetric stresses inserted in the virtual
element spaces; see also [11, 2, 43].

In this paper, we construct special virtual element spaces for the
approximation of solutions to elliptic problems, which have singularities due
to the geometry of the domain. Our approach falls within the broad family of
extended Galerkin methods based on the enrichment of local spaces with
singular functions, such as the extended finite element method (XFEM), see,
e.g., [37, 36], and the generalized finite element method (GFEM); see, e.g.,
[41]. These methods work as follows. Consider an elliptic problem on a
polygonal domain with smooth data. The solution to this problem has an a
priori known singular behaviour at the vertices of the domain; see, e.g., [32]
and the references therein. For this reason, a standard FEM converges to the
exact solution suboptimally. In order to cope with this suboptimality, in the
extended Galerkin methods, the approximation space is enriched with special
singular functions. A partition of unity is employed in order to patch the
local approximation spaces seamlessly. This enrichment permits to recover an
optimal convergence rate of the error of the method; see also [35].

In [18], the extended finite element setting of [37] is translated into the
virtual element one: local spaces consisting of polynomials plus singular
functions are patched with the aid of a virtual partition of unity, in the
spirit of [38].

Our approach is different and exploits the structure of virtual element
spaces. Instead of inserting the singular functions in the approximation
spaces explicitly and patching the local spaces with a partition of unity, we
proceed as follows. The singular functions, which are typically added to the
approximation spaces in the extended Galerkin methods, belong to the kernel of
the differential operator appearing in the problem under consideration. Such
singular functions can be inserted into local virtual element spaces, using
the fact that they are defined as solutions to local problems with data in the
finite dimensional spaces. By suitably tuning the boundary conditions in such
local spaces, we include the singular functions implicitly. Eventually, the
local spaces are patched in a nonconforming fashion.

A first advantage of our approach resides in the flexibility of using
polygonal meshes. On the other hand, the analysis and the implementation of
the method hinge upon a minor modification of what is done in the nonenriched
nonconforming VEM; see, e.g., [5]. Furthermore, the structure of nonconforming
spaces allows for the use of techniques suited to damp the ill-conditioning,
which typically arises in the extended Galerkin methods.

The method presented in this paper can be extended to more general problems,
such as linear and nonlinear elasticity problems. The extension to the three
dimension VEM is straightforward, thanks to the nonconforming structure of the
space; see [23].

The enriched virtual element method is based on two main ingredients: local
stabilizations and projections onto bulk and face enriched polynomial spaces.
We develop the analysis of the method for arbitrary polynomial order.

Importantly, we mention that there are two families of VEM: conforming and
nonconforming. For technical reasons, it appears that the latter family is
more suited to the enrichment we are going to present. This is a relevant
fact, for the nonconforming VEM is equivalent to Hybrid-High order methods and
Hybridizable Discontinuos Galerkin; see, e.g., [26, 30, 31]. Thus, the
proposed enrichment goes beyond the scope of the nonconforming VEM and could
be investigated in other settings as well.

##### Structure of the paper.

In Section 2, we present the model problem and recall regularity results for
elliptic partial differential equations on polygonal domain: we focus on the
case, where the singularities attain at the corner of the domain. We devote
Section 3 to the design of the enriched virtual element method. Its
theoretical analysis is the topic of Section 4. Here, we also discuss some
generalizations of the method. The theoretical results are validated by
several numerical experiments in Section 5, including an orthonormalization
procedure dramatically trimming the ill-conditioning of the final system. We
draw some conclusions in Section 6 and provide the implementation details in
Appendices A, B, and C.

##### Notation.

We employ a standard notation for Sobolev spaces. Given
$D\subset\mathbb{R}^{2}$ a domain and $s~{}\in~{}\mathbb{N}$, we denote the
standard Sobolev space of integer order $s$ over $D$ by $H^{s}(D)$. The case
$s=0$ is special: the Sobolev space $H^{0}(D)$ is the Lebesgue space
$L^{2}(D)$. We endow the Sobolev spaces with the standard inner products and
seminorms $(\cdot,\cdot)_{s,D}$ and $|\cdot|_{s,D}$. and denote the Sobolev
norm of order $s$ by

$\|\cdot\|_{s,D}^{2}:=\sum_{\ell=0}^{s}\|\cdot\|^{2}_{\ell,D}.$

For $s=1$, it is convenient to write

$a^{D}(\cdot,\cdot):=(\cdot,\cdot)_{1,D}.$

Fractional Sobolev spaces can be defined in several ways. We use the
definition of finiteness of the Aronszajin-Gagliardo-Slobodeckij norm; see,
e.g., [29] and the references therein. In particular, for any sufficiently
smooth $v$ on $\partial D$, set

$|v|_{\frac{1}{2},\partial D}^{2}:=\int_{\partial D}\int_{\partial
D}\frac{|v(\xi)-v(\eta)|^{2}}{|\xi-\eta|^{2}}d\xi
d\eta,\quad\quad\|v\|^{2}_{\frac{1}{2},\partial D}:=\|v\|^{2}_{0,\partial
D}+|v|^{2}_{\frac{1}{2},\partial D}.$

We define the fractional Sobolev space $H^{\frac{1}{2}}(\partial D)$ as

$H^{\frac{1}{2}}(\partial D):=\left\\{v\in L^{2}(\partial D)\text{ such that
}\|v\|_{\frac{1}{2},\partial D}\text{ is finite}\right\\}.$

Negative Sobolev spaces are defined via duality. In particular, we introduce
$H^{-\frac{1}{2}}(\partial D)$ as the dual space of $H^{\frac{1}{2}}(\partial
D)$. This space and its norm read

$H^{-\frac{1}{2}}(\partial D):=(H^{\frac{1}{2}}(\partial
D))^{*},\quad\quad\|v\|_{-\frac{1}{2},\partial
D}:=\sup_{\|w\|_{\frac{1}{2},\partial D}=1}(v,w)_{0,\partial D}.$

##  2 Model problem and regularity of the solution

Let $\Omega\subset\mathbb{R}^{2}$ be a polygonal domain with boundary
$\Gamma=\overline{\Gamma_{N}}\cup\overline{\Gamma_{D}}$, where $\Gamma_{D}$ is
a closed set in the topology of $\Gamma$ and
$\Gamma_{D}\cap\Gamma_{N}=\emptyset$. Denote the outward normal unit vector of
$\Gamma$ by $\mathbf{n}_{\Omega}$. Let $f$ be an analytic source term on
$\Omega$, and $g_{D}$ and $g_{N}$ be piecewise smooth functions on
$\Gamma_{D}$ and $\Gamma_{N}$. We allow for slit domains, see Figure 1 (left),
and domains with internal cuts (or cracks), see Figure 1 (right).

$\mathbf{A}_{1}$$\mathbf{A}_{2}$$\mathbf{A}_{3}$$\mathbf{A}_{4}$$\mathbf{A}_{5}=$$\mathbf{A}_{6}$$\,\,\mathbf{A}_{7}$

$\mathbf{A}_{1}$$\mathbf{A}_{2}$$\mathbf{A}_{3}$$\mathbf{A}_{4}$$\mathbf{A}_{5}$$\mathbf{A}_{6}$

Figure 1: _Left panel_ : a slit square domain. _Right panel_ : a domain with
an internal crack. We highlight the vertices and tips of the domain $\Omega$
in bold letters. If a cut starts from the boundary of $\Omega$, then two
vertices share the same coordinates.

Consider the following 2D Poisson problem on $\Omega$:

$\begin{cases}\text{find }u\text{ such that}\\\ -\Delta u=f&\text{in
}\Omega\\\ \mathbf{n}_{\Omega}\cdot\nabla u=g_{N}&\text{on }\Gamma_{N}\\\
u=g_{D}&\text{on }\Gamma_{D}.\\\ \end{cases}$ (1)

Define

$\begin{split}&V_{g_{D}}:=H^{1}_{g_{D}}(\Omega):=\\{v\in H^{1}(\Omega)\mid
v=g_{D}\text{ on }\Gamma_{D}\\},\\\ &V_{0}:=H^{1}_{0}(\Omega):=\\{v\in
H^{1}(\Omega)\mid v=0\text{ on }\Gamma_{D}\\},\\\ &a(u,v)=\int_{\Omega}\nabla
u\cdot\nabla v\quad\forall u,\,v\in H^{1}(\Omega).\end{split}$ (2)

In weak formulation, problem (1) reads

$\begin{cases}\text{find }u\in V_{g_{D}}\text{ such that}\\\
a(u,v)=(f,v)_{0,\Omega}+(g_{N},v)_{0,\Gamma_{N}}\quad\forall v\in V_{0}.\\\
\end{cases}$ (3)

Even if the right-hand side and the boundary conditions $g_{D}$ and $g_{N}$
are (piecewise) analytic, solution $u$ to problem (3) is not analytic over
$\overline{\Omega}$ in general. More precisely, $u$ is the combination of an
analytic function and a series of singular terms associated with the corners
of the domain and the tips of the cracks; see, e.g., [40, 8, 32, 27] and the
references therein.

We recall such expansion. Let $N_{\Omega}$ be the number of vertices and tips
of the cracks of $\Omega$. Denote the set of such vertices and tips by
$\\{\mathbf{A}_{i}\\}_{i=1}^{N_{\Omega}}$ and the associated angles by
$\\{\omega_{i}\\}_{i=1}^{N_{\Omega}}$. When no confusion occurs, we call
vertex both a vertex and a tip. If a crack has one of the two tips on the
boundary of the domain $\Omega$, then two vertices share the same coordinates;
see, e.g., the vertices $\mathbf{A}_{5}$ and $\mathbf{A}_{7}$ in Figure 1
(left). The boundary conditions are imposed on the two lips of the cut
separately.

We say that vertex $\mathbf{A}_{i}$ is a D (N) vertex if
$\mathbf{A}_{i}\in\Gamma_{D}$ ($\Gamma_{N}$). Otherwise, we say that
$\mathbf{A}_{i}$ is a D-N vertex. Introduce the singular exponents

$\alpha_{i,j}=\begin{cases}j\frac{\pi}{\omega_{i}}&\text{if
}\mathbf{A}_{i}\text{ is either D or N}\\\
{\left({j-\frac{1}{2}}\right)}\frac{\pi}{\omega_{i}}&\text{if
}\mathbf{A}_{i}\text{ is D-N}\\\
\end{cases}\quad\quad\forall\,i=1,\dots,N_{\Omega},\quad j\in\mathbb{N}.$ (4)

To each vertex $\mathbf{A}_{i}$, $i=1,\dots,N_{\Omega}$, we associate the two
(oriented counterclockwise) adjacent edges $\Gamma_{i(1)}$ and
$\Gamma_{i(2)}$, and the local set of polar coordinates

$\mathbf{A}_{i}\quad\xrightarrow{\hskip
28.45274pt}\quad{\left({r_{i},\,\theta_{i}}\right)}.$ (5)

Next, we introduce the so-called singular functions. For all
$i=1,\dots,N_{\Omega}$ and $j\in\mathbb{N}$, if the singular exponent
$\alpha_{i,j}$ in (4) does not belong to $\mathbb{N}$, then we set

$S_{i,j}(r_{i},\theta_{i})=\begin{cases}r_{i}^{\alpha_{i,j}}\sin{\left({\alpha_{i,j}\theta_{i}}\right)}&\text{if
}\Gamma_{i(2)}\subset\Gamma_{D}\\\
r_{i}^{\alpha_{i,j}}\cos{\left({\alpha_{i,j}\theta_{i}}\right)}&\text{otherwise}.\\\
\end{cases}$ (6)

Instead, if $\alpha_{i,j}\in\mathbb{N}$, $i=1,\dots,N_{\Omega}$,
$j\in\mathbb{N}$, then we set

$S_{i,j}(r_{i},\theta_{i})=\begin{cases}r_{i}^{\alpha_{i,j}}{\left({\log(r_{i})\sin{\left({\alpha_{i,j}\theta_{i}}\right)}+\theta_{i}\cos(\alpha_{i,j}\theta_{i})}\right)}&\text{if
}\Gamma_{i(2)}\subset\Gamma_{D}\\\
r_{i}^{\alpha_{i,j}}{\left({\log(r_{i})\cos{\left({\alpha_{i,j}\theta_{i}}\right)}+\theta_{i}\sin(\alpha_{i,j}\theta_{i})}\right)}&\text{otherwise}.\\\
\end{cases}$ (7)

If $\mathbf{A}_{i}$ is either a D or an N vertex, we can easily check that

$S_{i,j}\in
H^{1+j\frac{\pi}{\omega_{i}}-\varepsilon}(\Omega)\quad\quad\forall\varepsilon>0\quad\text{arbitrarily
small}.$ (8)

Moreover, we have

$\Delta S_{i,j}=0\quad\quad\text{pointwise in $\Omega$}.$

######  Theorem 2.1.

Given $s>0$, assume that $f\in H^{s-1}(\Omega)$ and $g_{D}$ and $g_{N}$ are
piecewise in $H^{s+\frac{1}{2}}(\Gamma_{D})$ and in
$H^{s-\frac{1}{2}}(\Gamma_{N})$. Then, the following decomposition of the
solution to problem (3) is valid:

$u=u_{0}+\sum_{i=1}^{N_{\Omega}}\sum_{\alpha_{i,j}<s}c_{i,j}S_{i,j}(r_{i},\theta_{i}),$
(9)

where $u_{0}\in H^{1+s}(\Omega)$ and $c_{i,j}\in\mathbb{R}$;

###### Proof.

see, e.g., [6, 7] and the references therein. ∎

Theorem 2.1 states that the solution to problem (3) is not analytic in
general, but rather has a known singular behaviour at the vertices of the
domain and at the tips of the crack. Such singular behaviour depends on the
magnitude of the angles associated with the vertices of the domain, regardless
of the smoothness of the data.

For ease of presentation, in the remainder of the paper, we assume that $u$,
the solution to problem (3), is such that the series of singular functions in
(9) reduces to a single term associated with a single vertex $\mathbf{A}$. In
other words, we assume that $u$ decomposes into

$u=u_{0}+\mathscr{S}_{\mathbf{A}}.$ (10)

In (10), $u_{0}$ denotes an analytic function over $\overline{\Omega}$,
whereas $\mathscr{S}_{\mathbf{A}}$ is a singular function of the form either
(6) or (7) with singularity centred at the corner/tip $\mathbf{A}$. From (6)
and (7), we have that $\mathscr{S}_{\mathbf{A}}(\lambda
r,\theta)=\lambda^{\alpha}\mathscr{S}_{\mathbf{A}}(r,\theta)$, where
$\alpha>0$ depends on $\mathscr{S}_{\mathbf{A}}$ and, consequently, on the
geometry of $\Omega$. Additionally, we assume that $\Gamma_{N}=\emptyset$ and
$g_{D}=0$ as well; see Remark 3 for further comments on more general cases.

Thus, in weak formulation, the problem we aim to solve reads

$\begin{cases}\text{find }u\in V:=H^{1}_{0}(\Omega)\text{ such that}\\\
a(u,v)=(f,v)_{0,\Omega}\quad\forall v\in H^{1}_{0}(\Omega).\\\ \end{cases}$
(11)

The analysis of this paper can be generalized: in Section 4.8, we discuss how
to cope with nonhomogeneous boundary conditions; multiple singularities; 3D
problems; general elliptic operators.

We exhibit a couple of examples falling in the setting of assumption (10).

######  Example 2.1.

Let $\Omega$ be the L-shaped domain, see Figure 2 (left),

$\Omega=(-1,1)^{2}\setminus\left([0,1)\times(-1,0]\right).$ (12)

According to (6)-(7), the expected strongest singularity is located at the re-
entrant corner $\mathbf{A}=(0,0)$. Assuming that only Dirichlet boundary
conditions are imposed, it is of our interest to consider the case when the
singular function is

$\mathscr{S}_{\mathbf{A}}(r,\theta)=r^{\frac{2}{3}}\sin\left(\frac{2}{3}\theta\right),$
(13)

where $(r,\theta)$ are the polar coordinates at $\mathbf{A}$.

######  Example 2.2.

Let $\Omega$ be the unit square with an internal crack, see Figure 2 (right),

$\Omega=(0,1)^{2}\setminus\\{(x,y)\in\mathbb{R}^{2}\mid
x\in[1/4,\,3/4],\,y=x\\}.$ (14)

According to (6)-(7), the expected strongest singularities are located at the
two tips of the internal crack. Denote one of the two tips by $\mathbf{A}$ and
its polar coordinates by $(r,\theta)$. It is of our interest to consider the
case when $u$ is singular at $\mathbf{A}$ only and the singular function is
given by

$\mathscr{S}_{\mathbf{A}}(r,\theta)=r^{\frac{1}{2}}\sin\left(\frac{1}{2}\theta\right).$

$\mathbf{A}$

$\mathbf{A}$

Figure 2: _Left panel_ : the L-shaped domain in (12). _Right panel_ : the unit
square with an internal crack (14). We denote the vertex/tip, where we assume
the singularity takes place, by $\mathbf{A}$.

##  3 The enriched virtual element method

We devote this section to the design of the novel enriched virtual element
method. To this aim, we introduce sequences of (regular) polygons in Section
3.1. We define the enriched virtual element spaces in Section 3.2, and the
discrete bilinear forms and right-hand side in Section 3.3. Eventually, we
exhibit the enriched VEM in Section 3.4.

###  3.1 Regular polygonal decomposition

Here, we introduce sequences of (regular) polygons and some geometric
assumptions.

Consider a sequence $\\{\mathscr{T}_{n}\\}_{n}$ of nonoverlapping polygons
partitioning $\Omega$. For all $n\in\mathbb{N}$ and $K\in\mathscr{T}_{n}$,
denote the size of $K$ by $h_{K}$ and the mesh size function of
$\mathscr{T}_{n}$, i.e., the maximum of such local diameter, by $h$. For each
$n\in\mathbb{N}$, denote the set of vertices and edges of $\mathscr{T}_{n}$ by
$\mathscr{V}_{n}$ and $\mathscr{E}_{n}$, and the set of internal and boundary
edges by $\mathscr{E}_{n}^{I}$ and $\mathscr{E}_{n}^{B}$. For all
$K\in\mathscr{T}_{n}$, let $\mathscr{E}_{n}^{K}$ be the set of edges of
element $K$ and $\mathbf{x}_{K}$ its barycenter. Given $e\in\mathscr{E}_{n}$
an edge, denote its length by $h_{e}$.

Henceforth, we demand that the following assumptions on $\mathscr{T}_{n}$, for
all $n\in\mathbb{N}$, are valid: there exists a positive constant
$\gamma\in(0,1)$, such that

  * (A0)

for all couples $K_{1}$ and $K_{2}\in\mathscr{T}_{n}$, $\gamma h_{K_{1}}\leq
h_{K_{2}}\leq\gamma^{-1}h_{K_{1}}$;

  * (A1)

for all $K\in\mathscr{T}_{n}$, $K$ is star-shaped with respect to a ball of
radius larger than or equal to $\gamma\,h_{K}$;

  * (A2)

for all $e\in\mathscr{E}^{K}$, the length $h_{e}$ of $e$ is larger than or
equal to $\gamma\,h_{K}$.

We employ assumptions (A0)-(A2) in the analysis of the method; see Section 4
below. Following, e.g., [16, 21, 24], they could be weakened. For the sake of
simplicity, we stick here to standard geometric assumptions.

We allow for elements with a crack having (at most) one endpoint on the
boundary and do not consider the case of completely cracked elements.

##### Notation on normal unit vectors.

Given $K\in\mathscr{T}_{n}$, we denote its outward normal unit vector by
$\mathbf{n}_{K}$. Besides, to each $e\in\mathscr{E}_{n}$, we associate the
normal unit vector $\mathbf{n}_{e}$ once and for all. In general,
$\mathbf{n}_{e}$ does not coincide with $\mathbf{n}_{K}{}_{|e}$ necessarily.
However, $\mathbf{n}_{e}\cdot\mathbf{n}_{K}{}_{|e}=\pm 1$. If
$e\in\mathscr{E}_{n}^{B}$, then we fix
$\mathbf{n}_{e}:=\mathbf{n}_{\Omega}{}_{|e}$.

###  3.2 Nonconforming enriched virtual element spaces

The core idea behind the design of the enriched VEM (EVEM) is that the
singular functions appearing in the expansions (9) and (10) belong to the
kernel of the Laplace operator.

Subdivide $\mathscr{T}_{n}$ into three layers. The first layer
$\mathscr{T}_{n}^{1}$ consists of the polygons abutting or containing the
singular vertex $\mathbf{A}$, as well as polygons that are sufficiently close
to $\mathbf{A}$: given $\widetilde{\gamma}>0$,

$\mathscr{T}_{n}^{1}:=\left\\{K\in\mathscr{T}_{n}\text{ such that
}\text{dist}(\mathbf{A},\mathbf{x}_{K})\leq\widetilde{\gamma}\operatorname{diam}(\Omega)\right\\}.$
(15)

Further comments on this definition will be provided in Remark 6 below.

The second layer $\mathscr{T}_{n}^{2}$ consists of the polygons sharing at
least one edge with elements in the first layer, i.e.,

$\mathscr{T}_{n}^{2}:=\left\\{K\in\mathscr{T}_{n}\setminus\mathscr{T}_{n}^{1}\text{
such that there exists~{}$\widetilde{K}\in\mathscr{T}_{n}^{1}$ with
}\text{card}(\mathscr{E}^{K}\cap\mathscr{E}^{\widetilde{K}})>0\right\\}.$

We set the third layer $\mathscr{T}_{n}^{3}$ as the remainder of the elements
in $\mathscr{T}_{n}$:

$\mathscr{T}_{n}^{3}:=\left\\{K\in\mathscr{T}_{n}\setminus(\mathscr{T}_{n}^{1}\cup\mathscr{T}_{n}^{2})\right\\}.$
(16)

We consider an analogous splitting for the set of edges $\mathscr{E}_{n}$. In
particular, we subdivide $\mathscr{E}_{n}$ into two layers of edges. The first
one, $\mathscr{E}_{n}^{1}$, is the set of all the edges belonging to the
boundary of elements in $\mathscr{T}_{n}^{1}$:

$\mathscr{E}_{n}^{1}=\left\\{e\in\mathscr{E}_{n}\text{ such that there
exists~{}$K\in\mathscr{T}_{n}^{1}$ with }e\subset\partial K\right\\},$

whereas second layer consists of the remainder of the edges, i.e.,

$\mathscr{E}_{n}^{2}=\left\\{e\in\mathscr{E}_{n}\setminus\mathscr{E}_{n}^{1}\right\\}.$

We exhibit a couple of graphical examples of such layers.

######  Example 3.1.

Let $\Omega$ be the L-shaped domain split into a Cartesian mesh of $48$
elements; see Figure 2 (left). The re-entrant corner is the singular vertex.
We show the distributions of element and edge layers in Figure 3 (left) and
(right). We pick $\widetilde{\gamma}=1/10$ in (15).

$\mathbf{A}$

$\mathbf{A}$

Figure 3: _Left panel_ : in blue, the elements in the first element layer
$\mathscr{T}_{n}^{1}$; in green, the elements in the second element layer
$\mathscr{T}_{n}^{2}$; in red, the elements in the third layer
$\mathscr{T}_{n}^{3}$. _Right panel_ : in red, the edges in the first edge
layer $\mathscr{E}_{n}^{1}$; in blue, the edges in the second edge layer
$\mathscr{E}_{n}^{2}$. The domain is the L-shaped domain defined in (12). We
assume that the solution to problem (11) is singular only at the re-entrant
corner $\mathbf{A}$. In red, we depict the circumference of radius $1/10$
centred at $\mathbf{A}$. The parameter $\widetilde{\gamma}$ in (15) is $1/10$.

######  Example 3.2.

Let $\Omega$ be the unit square domain with an internal crack split into a
Cartesian mesh of $9$ elements; see Figure 2 (right). The top-right tip is the
singular vertex. We show the distributions of element and edge layers, in
Figure 4 (left) and (right). We pick $\widetilde{\gamma}=1/10$ in (15).

$\mathbf{A}$

$\mathbf{A}$

Figure 4: _Left panel_ : in blue, the elements in the first element layer
$\mathscr{T}_{n}^{1}$; in green, the elements in the second element layer
$\mathscr{T}_{n}^{2}$; in red, the elements in the third layer
$\mathscr{T}_{n}^{3}$. _Right panel_ : in red, the edges in the first edge
layer $\mathscr{E}_{n}^{1}$; in blue, the edges in the second edge layer
$\mathscr{E}_{n}^{2}$. The domain is the unit square domain with an internal
crack defined in (14). We assume that the solution to problem (11) is singular
only at the right-upper tip of the crack $\mathbf{A}$. In red, we depict the
circumference of radius $1/10$ centred at $\mathbf{A}$. The parameter
$\widetilde{\gamma}$ in (15) is $1/10$.

The top-right corner contains a crack, which has a tip on the boundary.

Henceforth, we fix a $p\in\mathbb{N}$, which will denote the standard
polynomial order of accuracy of the method, and introduce auxiliary functions
and spaces. Given $K\in\mathscr{T}_{n}$, define the bulk-scaled enriching
function

$\mathscr{S}_{\mathbf{A}}^{K}(r,\theta):=\mathscr{S}_{\mathbf{A}}\left(\frac{r}{h_{K}},\theta\right).$
(17)

For example, if $\mathscr{S}_{\mathbf{A}}$ is given as in (13), then
$\mathscr{S}_{\mathbf{A}}^{K}(r,\theta)=\left(\frac{r}{h_{K}}\right)^{\frac{2}{3}}\sin\left(\frac{2}{3}\theta\right)$.

Given $e\in\mathscr{E}^{K}$, define the edge-scaled enriching function

$\mathscr{S}_{\mathbf{A}}^{e}(r,\theta):=\mathscr{S}_{\mathbf{A}}\left(\frac{r}{h_{e}},\theta\right).$
(18)

For example, if $\mathscr{S}_{\mathbf{A}}$ is given as in (13), then
$\mathscr{S}_{\mathbf{A}}^{e}(r,\theta)=\left(\frac{r}{h_{e}}\right)^{\frac{2}{3}}\sin\left(\frac{2}{3}\theta\right)$.

Observe that

$\mathscr{S}_{\mathbf{A}}^{K}(r,\theta)=\left(\frac{h_{e}}{h_{K}}\right)^{\alpha}\mathscr{S}_{\mathbf{A}}^{e}(r,\theta),$
(19)

where $\alpha>0$ depends on the definition of $\mathscr{S}_{\mathbf{A}}$. For
example, if $\mathscr{S}_{\mathbf{A}}$ is given as in (13), then $\alpha=2/3$.

We define the set of enriched polynomials over an element
$K\in\mathscr{T}_{n}$ as follows:

$\widetilde{\mathbb{P}}_{p}(K)=\begin{cases}\mathbb{P}_{p}(K)\oplus\mathscr{S}_{\mathbf{A}}^{K}&\text{if
}K\in\mathscr{T}_{n}^{1}\\\ \mathbb{P}_{p}(K)&\text{otherwise}.\\\
\end{cases}$

In other words, on the elements close to the singular vertex/tip $\mathbf{A}$,
we add the singular function $\mathscr{S}_{\mathbf{A}}$ to the nonenriched
polynomial space $\mathbb{P}_{p}(K)$.

Further, we define the set of enriched polynomials over edges:

$\widetilde{\mathbb{P}}_{p-1}(e)=\begin{cases}\mathbb{P}_{p-1}(e)\oplus\left(\mathbf{n}_{e}\cdot\nabla\mathscr{S}_{\mathbf{A}}^{e}{}_{|e}\right)&\text{if
}e\in\mathscr{E}_{n}^{1}\\\ \mathbb{P}_{p-1}(e)&\text{otherwise}.\\\
\end{cases}$ (20)

In other words, we consider nonenriched one dimensional polynomial spaces
$\mathbb{P}_{p-1}(e)$ on all edges except those belonging to the boundary of
the elements in the first layer $\mathscr{T}_{n}^{1}$. Here, we consider the
normal derivative of the scaled enriching function
$\mathscr{S}_{\mathbf{A}}^{e}$ as additional special function.

Next, we define the local enriched virtual element spaces: for all
$K\in\mathscr{T}_{n}$,

$V_{n}(K):={\left\\{{v_{n}\in H^{1}(K)\mid\Delta
v_{n}\in\mathbb{P}_{p-2}(K),\;\mathbf{n}_{e}\cdot\nabla
v_{n}{}_{|e}\in\widetilde{\mathbb{P}}_{p-1}(e)\quad\forall
e\in\mathscr{E}^{K}}\right\\}}.$ (21)

Functions in $V_{n}(K)$ are unknown in closed form both in the bulk and the
boundary of element $K$. This is the reason why the functions are referred to
as virtual. The space $V_{n}(K)$ contains the space of polynomials of degree
$p$. Furthermore, if $K\in\mathscr{T}_{n}^{1}$, then the singular function
$\mathscr{S}_{\mathbf{A}}$ belongs to $V_{n}(K)$ as well. This is the reason
why we regard space $V_{n}(K)$ as enriched.

We introduce additional notation. Let
${\left\\{{m_{\alpha}}\right\\}}_{|\boldsymbol{\alpha}|=0}^{\dim(\mathbb{P}_{p-2}(K))}$
be a basis of $\mathbb{P}_{p-2}(K)$. For instance, this basis consists of the
monomials introduced in [12] or some orthonormal basis as in [33]. We assume
that the elements $m_{\alpha}$ are centred in the barycenter of the element
and scaled according to the element diameter. It is known, that the
$m_{\alpha}$ basis can be $L^{2}$-orthonormalized for stability purposes; see,
e.g., [33]. For the sake of exposition, we stick here to the monomial basis.
Besides, let
${\left\\{{\widetilde{m}_{\alpha}^{e}}\right\\}}_{\alpha=0}^{p-1}(e)$ be a
basis of $\widetilde{\mathbb{P}}_{p-1}(e)$ defined in (20) for all edges
$e\in\mathscr{E}^{K}$. A possible choice of the basis is provided by the first
$p-1$ Legendre polynomials on the local system of coordinates over $e$, and
the global normal derivative of the singular function
$\mathscr{S}_{\mathbf{A}}^{e}$ over edge $e$. Such basis can be
orthonormalized, leading to a dramatic improvement of the performance of the
method; see Sections 5.3 and Appendix B below. Here, for the sake of
presentation, we stick to the former choice.

Consider the following set of linear functionals on $V_{n}(K)$: for all
$v_{n}\in V_{n}(K)$,

  * •

the internal moments:

$\frac{1}{|K|}\int_{K}v_{n}m_{\alpha}\quad\forall\alpha=1,\dots,\dim(\mathbb{P}_{p-2}(K));$
(22)

  * •

if $e\in\mathscr{E}_{n}^{1}$, the edge moments:

$\begin{cases}\frac{1}{h_{e}}\int_{e}v_{n}\widetilde{m}_{\alpha}^{e}\quad\forall\alpha=0,\dots,p-1,\,\forall
e\in\mathscr{E}^{K},\\\
\int_{e}v_{n}(\mathbf{n}_{e}\cdot\nabla\mathscr{S}_{\mathbf{A}}^{e}){}_{|e}\quad\alpha=p\,,\,\forall
e\in\mathscr{E}^{K}.\\\ \end{cases}$ (23)

If $e\in\mathscr{E}_{n}^{2}$, then the edge moments are the same, but there is
no special moment for $\alpha=p$.

######  Remark 1.

The enriched edge functionals in (23) are well posed for all possible singular
functions $\mathscr{S}_{\mathbf{A}}$. In order to see this, we first observe
that $(\mathbf{n}_{e}\cdot\nabla\mathscr{S}_{\mathbf{A}}^{e}){}_{|e}\in
L^{1}(e)$ for all $e\in\mathscr{E}_{n}$. In fact, in the worst case scenario,
i.e., the case of a slit domain with mixed Dirichlet and Neumann boundary
conditions, $\mathscr{S}_{\mathbf{A}}$ belongs to
$H^{\frac{5}{4}-\varepsilon}(\Omega)$, $\varepsilon>0$.

On the other hand, given $K\in\mathscr{T}_{n}$, any function $v_{n}$ in the
local virtual element space $V_{n}(K)$ solves a local elliptic problem. In
particular, $v_{n}\in H^{1+\varepsilon}(K)$, $\varepsilon>0$. Therefore, the
Sobolev embedding theorem in two dimensions yields
$v_{n}\in\mathscr{C}^{0}(\overline{K})$ and
$v_{n}{}_{|e}\in\mathscr{C}^{0}(e)$ for all $e\in\mathscr{E}^{K}$. The well
posedness of (23) follows.

This fact entails that all the forthcoming integrations by parts involving
functions $v_{n}\in V_{n}(K)$ are well defined; see, e.g., (24) and (44).

######  Lemma 3.1.

For all $K\in\mathscr{T}_{n}$, the set of linear functionals in (22)-(23) is a
unisolvent set of degrees of freedom for space $V_{n}(K)$.

###### Proof.

The dimension of $V_{n}(K)$ is equal to the number of linear functionals; see
[5]. Thence, it suffices to prove the unisolvence of such functionals. Observe
that

$|v_{n}|^{2}_{1,K}=-\int_{K}\underbrace{\Delta
v_{n}}_{\in\mathbb{P}_{p-2}(K)}\,v_{n}+\sum_{e\in\mathscr{E}^{K}}\int_{e}\underbrace{\mathbf{n}_{K}\cdot\nabla
v_{n}}_{\in\widetilde{\mathbb{P}}_{p-1}(e)}\,v_{n}.$ (24)

The first and second terms are zero, for the internal (22) and the edge
moments (23) are zero by assumption. Hence, $v_{n}$ is constant. This and the
fact that the average over $\partial K$ of $v_{n}$ is equal to zero entail the
assertion. ∎

Compared to the degrees of freedom (DOFs) in the nonenriched nonconforming VEM
[5], we consider the same internal DOFs (22). As for the edge DOFs (23), we
cope with additional moments related to the special functions on the edges of
the elements in the first edge layer $\mathscr{E}_{n}^{1}$. In Figures 5 and
6, we depict the DOFs for the nonenriched and enriched nonconforming virtual
element method, with order $p=1$ and $2$ on a pentagon.

$p=1$, nonenriched space

$p=1$, enriched space

Figure 5: Degrees of freedom on a pentagon for $p=1$. _Left panel_ :
nonenriched nonconforming VEM. _Right panel_ : enriched nonconforming VEM. The
blue circles represent polynomial moments on the edges. The green pentagons
represent the enriched edge moments.

$p=2$, nonenriched space

$p=2$, enriched space

Figure 6: Degrees of freedom on a pentagon for $p=2$. _Left panel_ :
nonenriched nonconforming VEM. _Right panel_ : enriched nonconforming VEM. The
blue circles represent the polynomial moments on the edges. The green
pentagons represent the enriched edge moments. The magenta triangle represents
the only polynomial bulk moment.

For future convenience, introduce the local canonical basis
$\\{\varphi_{i}\\}_{i=1}^{\dim(V_{n}(K))}$ defined as

$\text{dof}_{j}(\varphi_{i})=\delta_{i,j},$ (25)

where $\delta_{i,j}$ denotes the Kronecker delta.

Next, split the bilinear form $a(\cdot,\cdot)$ defined in (2) into local
contributions:

$a(u,v)=\sum_{K\in\mathscr{T}_{n}}a^{K}(u,v):=\sum_{K\in\mathscr{T}_{n}}\int_{K}\nabla
u_{|K}\cdot\nabla v_{|K}\quad\quad\forall u,\,v\in H^{1}(\Omega).$

The definition of the degrees of freedom allows for the computation the
enriched $H^{1}$-orthogonal projection
$\widetilde{\Pi}^{\nabla,K}_{p}:V_{n}(K)\rightarrow\widetilde{\mathbb{P}}_{p}(K)$:

$\begin{cases}a^{K}(v_{n}-\widetilde{\Pi}^{\nabla,K}_{p}v_{n},q_{p})=0\\\
\int_{\partial K}(v_{n}-\widetilde{\Pi}^{\nabla,K}_{p}v_{n})=0\\\
\end{cases}\quad\quad\forall v_{n}\in V_{n}(K),\quad\forall
q_{p}\in\widetilde{\mathbb{P}}_{p}(K).$ (26)

In fact, an integration by parts yields

$a^{K}(v_{n},q_{p})=(v_{n},\underbrace{-\Delta
q_{p}}_{\in\mathbb{P}_{p-2}(K)})_{0,K}+\sum_{e\in\mathscr{E}^{K}}(v_{n},\underbrace{\mathbf{n}_{K}\cdot\nabla
q_{p}}_{\in\widetilde{\mathbb{P}}_{p-1}(e)})_{0,e}.$ (27)

The first and second term on the right-hand side of (27) are computable from
(22) and (23). The second term can be approximated at any precision by a one
dimensional quadrature formula; see Remark 10 below for more comments on this
point. For all the elements $K\in\mathscr{T}_{n}^{1}$, the projection
$\widetilde{\Pi}^{\nabla,K}_{p}$ maps functions belonging to the virtual
element space into the space of bulk-enriched polynomial spaces
$\widetilde{\mathbb{P}}_{p-2}(K)$.

Further, for all $e\in\mathscr{E}_{n}$, we consider the possibly enriched
$L^{2}$ edge projector
$\widetilde{\Pi}^{0,e}_{p-1}:V_{n}(K){}_{|e}\rightarrow\widetilde{\mathbb{P}}_{p-1}(e)$,
defined as

$\int_{e}(v_{n}-\widetilde{\Pi}^{0,e}_{p-1}v_{n})\widetilde{m}_{\alpha}^{e}=0\quad\quad\forall
v_{n}\in
V_{n}(K),\quad\forall\widetilde{m}_{\alpha}^{e}\in\widetilde{\mathbb{P}}_{p-1}(e).$
(28)

The computability of such projector follows from the definition of the edge
degrees of freedom (23).

######  Remark 2.

Let $e\in\mathscr{E}_{n}$ and $K$ be such that $e\in\mathscr{E}^{K}$. The
projector $\widetilde{\Pi}^{0,e}_{p-1}$ in (28) can be computed if
$\mathscr{S}_{\mathbf{A}}\in H^{\frac{3}{2}+\varepsilon}(K)$ with
$\varepsilon>0$. In fact, we need
$\mathbf{n}_{K}\cdot\nabla\mathscr{S}_{\mathbf{A}}{}_{|e}\in L^{2}(e)$. This
appears as a partial limitation to design and applicability of the method.
However, we need to compute projector $\widetilde{\Pi}^{0,e}_{p-1}$ in the
following circumstances only: in the design of a theoretical stabilization for
the method; see Section 4.6 below; in the discretization of certain
nonhomogeneous Neumann boundary conditions; see Remark 3 and Appendix A.2
below; in the “orthonormalized” version of the method; see Appendix B below.

We also introduce the vector nonenriched polynomial projections
$\boldsymbol{\Pi}^{\nabla,K}_{p}:[H^{1}(K)]^{2}\rightarrow[\mathbb{P}_{p}(K)]^{2}$
and
$\boldsymbol{\Pi}^{0,e}_{p-1}:[L^{2}(e)]^{2}\rightarrow[\mathbb{P}_{p-1}(e)]^{2}$
similarly as in (26) and (28).

Eventually, we define the nonenriched $L^{2}$-bulk orthogonal projector
$\Pi^{0,K}_{p-2}:V_{n}(K)\rightarrow\mathbb{P}_{p-2}(K)$ as

$(v_{n}-\Pi^{0,K}_{p-2}v_{n},q_{p-2})_{0,K}=0\quad\quad\forall
q_{p-2}\in\mathbb{P}_{p-2}(K).$ (29)

The projector $\Pi^{0,K}_{p-2}$ is computable from the internal degrees of
freedom (22) and is used for the approximation of the Neumann boundary
conditions only; see Remark 3 below.

Next, we define the global nonconforming virtual element space $V_{n}$. Given
an internal edge $e\in\mathscr{E}_{n}^{I}$, denote its two adjacent elements
by $K^{+}$ and $K^{-}$. Instead, given a boundary edge
$e\in\mathscr{E}_{n}^{B}$, denote its adjacent element by $K$. Moreover,
denote the space of $L^{2}(\Omega)$ functions piecewise in $H^{1}$ over
$\mathscr{T}_{n}$ by $H^{1}(\mathscr{T}_{n})$, and define the broken Sobolev
norm

$|v|^{2}_{1,\mathscr{T}_{n}}:=\sum_{K\in\mathscr{T}_{n}}|v_{K}|^{2}_{1,K}\quad\quad\forall
v\in H^{1}(\mathscr{T}_{n}).$

Introduce the jump operator across an edge $e\in\mathscr{E}_{n}^{1}$: given
$v\in H^{1}(\mathscr{T}_{n})$, set

$\llbracket v\rrbracket_{e}=\llbracket
v\rrbracket:=\begin{cases}v_{|K^{+}}\mathbf{n}_{K^{+}}+v_{|K^{-}}\mathbf{n}_{K^{-}}&\text{if
}e\in\mathscr{E}_{n}^{I}\\\ v\mathbf{n}_{K}&\text{if
}e\in\mathscr{E}_{n}^{B}.\\\ \end{cases}$ (30)

Introduce the global nonconforming Sobolev space of order $p$, subordinated to
the mesh $\mathscr{T}_{n}$, including homogeneous boundary conditions in a
nonconforming sense:

$\begin{split}H^{1,nc}_{0}(\mathscr{T}_{n},p):=&\left\\{v\in
H^{1}(\mathscr{T}_{n})\middle|\int_{e}\llbracket
v\rrbracket\cdot\mathbf{n}_{e}\,\widetilde{m}_{\alpha}^{e}=0\quad\forall\widetilde{m}_{\alpha}^{e}\in\widetilde{\mathbb{P}}_{p-1}(e),\;\forall
e\in\mathscr{E}_{n}\right\\}.\\\ \end{split}$

We define the global test and trial nonconforming enriched virtual element
spaces as

$V_{n}:=\left\\{v_{n}\in H^{1,nc}_{0}(\mathscr{T}_{n},p)\mid v_{n}{}_{|K}\in
V_{n}(K)\;\,\forall K\in\mathscr{T}_{n}\right\\}.\\\ $ (31)

We construct space $V_{n}$ by a nonconforming coupling of the local edge
degrees of freedom (23). The global canonical basis is defined from its local
counterparts (25) accordingly.

###  3.3 The discrete bilinear form and right-hand side

Functions in the virtual element spaces are not available in closed form.
Thus, in order to design the numerical scheme for the approximation of
solutions to (11), we introduce a global bilinear form and right-hand side
that are computable in terms of the degrees of freedom. To this purpose, we
generalize the construction in [5] to the enriched setting.

##### The discrete bilinear form.

Using the orthogonality property of projector $\widetilde{\Pi}^{\nabla,K}_{p}$
in (26), we apply Pythagoras’ theorem in Hilbert spaces and get

$a^{K}(u,v)=a^{K}(\widetilde{\Pi}^{\nabla,K}_{p}u,\widetilde{\Pi}^{\nabla,K}_{p}v)+a^{K}((I-\widetilde{\Pi}^{\nabla,K}_{p})u,(I-\widetilde{\Pi}^{\nabla,K}_{p})v)\quad\forall
u,\,v\in H^{1}(K).$

The first term on the right-hand side is computable on $V_{n}(K)\times
V_{n}(K)$, see (27), whereas the second one is not. As standard in virtual
elements [12, 13], we introduce a symmetric bilinear form
$S^{K}:\ker(\widetilde{\Pi}^{\nabla,K}_{p})\times\ker(\widetilde{\Pi}^{\nabla,K}_{p})\rightarrow\mathbb{R}$
satisfying

$c_{*}(K)|v_{n}|^{2}_{1,K}\leq S^{K}(v_{n},v_{n})\leq
c^{*}(K)|v_{n}|^{2}_{1,K}\quad\quad\forall
v_{n}\in\ker(\widetilde{\Pi}^{\nabla,K}_{p}),\quad\forall
K\in\mathscr{T}_{n},$ (32)

where $c_{*}(K)$ and $c^{*}(K)$ are two positive constants, possibly depending
on the polynomial degree $p$, the geometric properties of $K$, and the
singular function $\mathscr{S}_{\mathbf{A}}$.

Having this at hand, define the local discrete bilinear forms

$a_{n}^{K}(u_{n},v_{n})=a^{K}(\widetilde{\Pi}^{\nabla,K}_{p}u_{n},\widetilde{\Pi}^{\nabla,K}_{p}v_{n})+S^{K}((I-\widetilde{\Pi}^{\nabla,K}_{p})u_{n},(I-\widetilde{\Pi}^{\nabla,K}_{p})v_{n})\quad\quad\forall
u_{n},\,v_{n}\in V_{n}(K),$

and the global discrete bilinear form

$a_{n}(u_{n},v_{n})=\sum_{K\in\mathscr{T}_{n}}a_{n}^{K}(u_{n}{}_{|K},v_{n}{}_{|K})\quad\quad\forall
u_{n},v_{n}\in V_{n}.$

We postpone explicit choices and further considerations about the
stabilization forms $S^{K}$ to Section 4.6 below.

For all $K\in\mathscr{T}_{n}$, the local discrete bilinear forms $a_{n}^{K}$
are coercive and continuous with respect to the $H^{1}$ seminorm, with
coercivity and continuity constants

$\alpha_{*}(K)=\min(1,c_{*}(K)),\quad\quad\alpha^{*}(K)=\max(1,c^{*}(K)).$

In other words, for all $K\in\mathscr{T}_{n}$, we have

$\alpha_{*}(K)|v_{n}|^{2}_{1,K}\leq
a_{n}^{K}(v_{n},v_{n})\leq\alpha^{*}(K)|v_{n}|^{2}_{1,K}\quad\quad\forall
v_{n}\in V_{n}(K).$ (33)

Moreover, the stabilization is symmetric and polynomially-enriched consistent:
for all $K\in\mathscr{T}_{n}$,

$a_{n}^{K}(q_{p},v_{n})=a^{K}(q_{p},v_{n})\quad\quad\forall
q_{p}\in\widetilde{\mathbb{P}}_{p}(K),\quad\forall v_{n}\in V_{n}.$ (34)

##### The discrete right-hand side.

Denote the number and the set of vertices of element $K$ by $N_{V}^{K}$ and
$\\{\nu_{i}\\}_{i=1}^{N_{V}^{K}}$, recall the definition of projector
$\Pi^{0,K}_{p-2}$ in (29), and define

$\langle f,v_{n}\rangle_{n}=\sum_{K\in\mathscr{T}_{n}}\langle
f_{|K},v_{n}{}_{|K}\rangle_{n,K}:=\begin{cases}\sum_{K\in\mathscr{T}_{n}}\int_{K}f_{|K}\Pi^{0,K}_{p-2}v_{n}{}_{|K}&\text{if
}p\geq 2\\\
\sum_{K\in\mathscr{T}_{n}}\frac{1}{N_{V}^{K}}\int_{K}f_{|K}(\sum_{i=1}^{N_{V}}v_{n}(\nu_{i}))&\text{if
}p=1.\\\ \end{cases}$

###  3.4 The method

The nonconforming enriched virtual element method for problem (11) reads

$\begin{cases}\text{find }u_{n}\in V_{n}\text{ such that}\\\
a_{n}(u_{n},v_{n})=\langle f,v_{n}\rangle_{n}\quad\quad\forall v_{n}\in
V_{n}.\end{cases}$ (35)

Method (35) is well posed thanks to the continuity and coercivity properties
detailed in Section 3.3. We devote Section 4 below to the analysis of (35),
whereas we provide the implementation details in Appendix A.

######  Remark 3.

Nonhomogeneous Neumann boundary conditions $g_{N}$ in (3) are approximated as
follows:

$\sum_{e\in\mathscr{E}_{n}^{B},\,e\subset\Gamma_{N}}\int_{e}g_{N}\widetilde{\Pi}^{0,e}_{p-1}v_{n}.$
(36)

As highlighted in Remark 2, for all $e\in\mathscr{E}_{n}$, projector
$\widetilde{\Pi}^{0,e}_{p-1}$ can be computed if $\mathscr{S}_{\mathbf{A}}\in
H^{\frac{3}{2}+\varepsilon}(K)$ with $\varepsilon>0$, where
$K\in\mathscr{T}_{n}$ is such that $e\in\mathscr{E}^{K}$. Further comments on
nonhomogenous Neumann boundary conditions and a generalization of (36) for
some special classes of $g_{N}$ are provided in Appendix A.2 below.

Nonhomogeneous Dirichlet boundary conditions are enforced through the degrees
of freedom on the Dirichlet edges in the trial and test spaces.

##### The role of the Heaviside function.

Our approach does not allow for the inclusion of the Heaviside in the test and
trial spaces as done in the XFEM, GFEM, and XVEM. In fact, the Laplacian of
the Heaviside function is not a function, rendering the construction of
$V_{n}(K)$ in (21) not feasible. However, the flexibility of polygonal meshes
renders the use of the Heaviside function practically useless. In fact,
elements cut into two by a crack can be re-meshed into distinct polygons,
whereas the setting in Section 3.2 allows for handling automatically elements
with internal cracks; see also [4].

##  4 Error analysis

In this section, we analyse the rate of convergence of method (35). We present
an abstract error result in Section 4.1. In Sections 4.2 and 4.3, we describe
the approximation properties of enriched polynomial and enriched virtual
element spaces. On the other hand, Sections 4.4 and 4.5 deal with the
approximation of the variational crimes perpetrated in the approximation of
the right-hand side and the nonconformity of the method. After introducing and
analysing possible stabilizations in Section 4.6, we collect all the above
estimates in Section 4.7. We discuss some extensions and generalization in
Section 4.8.

###  4.1 Abstract error analysis

Here, we present the abstract error analysis result for method (35). Recall
that jump operator $\llbracket\cdot\rrbracket_{e}$ across edge $e$ is defined
in (30), and define the bilinear form $\mathscr{N}_{n}:H^{1}(\Omega)\times
H^{1,nc}_{0}(\mathscr{T}_{n},p)\rightarrow\mathbb{R}$ as

$\mathscr{N}_{n}(u,v)=\sum_{e\in\mathscr{E}_{n}}\int_{e}\nabla
u\cdot\llbracket v\rrbracket_{e}.$ (37)

######  Theorem 4.1.

Let $u$ and $u_{n}$ be the solutions to (11) and (35). Recall that
$\alpha_{*}(K)$ and $\alpha^{*}(K)$ are the stability constants in (33), and
$\mathscr{N}_{n}$ is defined in (37). Then, the following a priori estimate is
valid: for all $u_{\pi}$ piecewise in $\widetilde{\mathbb{P}}_{p}(K)$ and for
all $u_{I}\in V_{n}$,

$\begin{split}|u-u_{n}|_{1,\mathscr{T}_{n}}\leq&\max_{K\in\mathscr{T}_{n}}\alpha_{*}^{-1}(K)\left\\{\sup_{v_{n}\in
V_{n}}\frac{\langle
f,v_{n}\rangle_{n}-(f,v_{n})_{0,\Omega}}{|v_{n}|_{1,\mathscr{T}_{n}}}+\sup_{v_{n}\in
V_{n}}\frac{\mathscr{N}_{n}(u,v_{n})}{|v_{n}|_{1,\mathscr{T}_{n}}}\right.\\\
&\quad\quad\left.+\left(1+\max_{K\in\mathscr{T}_{n}}\alpha^{*}(K)\right)\left(|u-u_{\pi}|_{1,\mathscr{T}_{n}}+|u-u_{I}|_{1,\mathscr{T}_{n}}\right)\right\\}.\\\
\end{split}$ (38)

###### Proof.

The proof follows along the same lines as that in [5, Theorem 4.1]. For the
sake of completeness, we carry out all the details.

For all $u_{I}\in V_{n}$, the triangle inequality entails, for
$\delta_{n}:=u_{n}-u_{I}$,

$|u-u_{n}|_{1,\mathscr{T}_{n}}\leq|u-u_{I}|_{1,\mathscr{T}_{n}}+|u_{I}-u_{n}|_{1,\mathscr{T}_{n}}=|u-u_{I}|_{1,\mathscr{T}_{n}}+|\delta_{n}|_{1,\mathscr{T}_{n}}.$
(39)

We estimate the second term on the right-hand side: for all $u_{\pi}$
piecewise in $\widetilde{\mathbb{P}}_{p}(K)$, apply (33) and (34) and get

$\begin{split}|\delta_{n}|_{1,\mathscr{T}_{n}}^{2}&=\sum_{K\in\mathscr{T}_{n}}|\delta_{n}|^{2}_{1,\mathscr{T}_{n}}\leq\sum_{E\in\mathscr{T}_{n}}\alpha_{*}^{-1}(K)a_{n}^{K}(\delta_{n},\delta_{n})\\\
&\leq\max_{K\in\mathscr{T}_{n}}\alpha_{*}^{-1}(K)\sum_{K\in\mathscr{T}_{n}}\left\\{a_{n}^{K}(u_{n},\delta_{n})-a_{n}^{K}(u_{I}-u_{\pi},\delta_{n})-a^{K}(u_{\pi}-u,\delta_{n})-a^{K}(u,\delta_{n})\right\\}.\\\
\end{split}$

Observe that

$\begin{split}\sum_{K\in\mathscr{T}_{n}}a^{K}(u,\delta_{n})&=\sum_{K\in\mathscr{T}_{n}}\left\\{\int_{K}-\Delta
u\,\delta_{n}+\int_{\partial K}\mathbf{n}_{K}\cdot\nabla
u\,\delta_{n}\right\\}\\\
&=(f,\delta_{n})_{0,\Omega}+\sum_{e\in\mathscr{E}_{n}}\int_{e}\nabla
u\cdot\llbracket\delta_{n}\rrbracket=:(f,\delta_{n})_{0,\Omega}+\mathscr{N}_{n}(u,\delta_{n}).\\\
\end{split}$

Deduce that

$\begin{split}|\delta_{n}|^{2}_{1,\mathscr{T}_{n}}&\leq\max_{K\in\mathscr{T}_{n}}\alpha_{*}^{-1}(K)\Big{\\{}\langle
f,\delta_{n}\rangle_{n}-(f,\delta_{n})_{0,\Omega}-\mathscr{N}_{n}(u,\delta_{n})\\\
&\quad\quad\quad\quad\quad\quad\quad\quad\left.+\max_{K\in\mathscr{T}_{n}}\alpha^{*}(K)|u_{I}-u_{\pi}|_{1,\mathscr{T}_{n}}|\delta_{n}|_{1,\mathscr{T}_{n}}+|u-u_{\pi}|_{1,\mathscr{T}_{n}}|\delta_{n}|_{1,\mathscr{T}_{n}}\right\\},\\\
\end{split}$

which implies

$\begin{split}|\delta_{n}|_{1,\mathscr{T}_{n}}&\leq\max_{K\in\mathscr{T}_{n}}\alpha_{*}^{-1}(K)\Big{\\{}\frac{\langle
f,\delta_{n}\rangle_{n}-(f,\delta_{n})_{0,\Omega}}{|\delta_{n}|_{1,\mathscr{T}_{n}}}-\frac{\mathscr{N}_{n}(u,\delta_{n})}{|\delta_{n}|_{1,\delta_{n}}}\\\
&\quad\quad\quad\quad\quad\quad\quad\quad\left.+(1+\max_{K\in\mathscr{T}_{n}}\alpha^{*}(K))|u-u_{\pi}|_{1,\mathscr{T}_{n}}+\max_{K\in\mathscr{T}_{n}}\alpha^{*}(K)|u-u_{I}|_{1,\mathscr{T}_{n}}\right\\}.\\\
\end{split}$

This, together with (39), yields the assertion. ∎

###  4.2 Best enriched polynomial approximation estimates

We show how to estimate from above the term $|u-u_{\pi}|_{1,\mathscr{T}_{n}}$
on the right-hand side of (38) for a specific choice of $u_{\pi}$ piecewise in
$\widetilde{\mathbb{P}}_{p}(K)$.

######  Lemma 4.2.

Let $u$ be the solution to problem (11), $u_{0}$ be as in (10), and
assumptions (A0)-(A2) be valid. Then, there exists $u_{\pi}$ piecewise in
$\widetilde{\mathbb{P}}_{p}(K)$ such that

$|u-u_{\pi}|_{1,\mathscr{T}_{n}}\leq
ch^{p}\left\\{\left(\sum_{K\in\mathscr{T}_{n}^{1}}|u_{0}|_{p+1,K}^{2}\right)^{\frac{1}{2}}+\left(\sum_{K\in\mathscr{T}_{n}^{2}\cup\mathscr{T}_{n}^{3}}|u|_{p+1,K}^{2}\right)^{\frac{1}{2}}\right\\},$

where $c$ is a positive constant depending on $p$ and $\gamma$, being $\gamma$
introduced in (A0)-(A2), but is independent of $h$ and $u$.

###### Proof.

For all $K\notin\mathscr{T}_{n}^{1}$, we have
$\widetilde{\mathbb{P}}_{p}(K)=\mathbb{P}_{p}(K)$. Therefore, we pick
$u_{\pi}$ as the best piecewise $H^{1}(K)$ polynomial approximant of $u$.
Deduce that

$|u-u_{\pi}|_{1,K}=\inf_{q_{p}\in\mathbb{P}_{p}(K)}|u-q_{p}|_{1,K}\leq
ch_{K}^{p}|u|_{p+1,K}.$ (40)

This is a consequence of the smoothness of $u$ on all
$K\notin\mathscr{T}_{n}^{1}$ and standard polynomial best approximation
estimates; see, e.g., [20]. The constant $c$ depends on the order of accuracy
of the method $p$ and on the shape of element $K$.

If $K\in\mathscr{T}_{n}^{1}$, then,
$\widetilde{\mathbb{P}}_{p}(K)\supsetneqq\mathbb{P}_{p}(K)$: the former space
is spanned by the latter _plus_ the singular function
$\mathscr{S}_{\mathbf{A}}^{K}$. Thus, a suitable choice of $u_{\pi}$ is given
by a combination of the singular function $\mathscr{S}_{\mathbf{A}}^{K}$ and
the best $H^{1}(K)$ polynomial approximant of $u_{0}$, being $u_{0}$
introduced in (10). This entails, for some constants
$\widetilde{c}\in\mathbb{R}$,

$|u-u_{\pi}|_{1,K}=\inf_{q_{p}\in\mathbb{P}_{p}(K)}|u_{0}-q_{p}+\widetilde{c}\mathscr{S}_{\mathbf{A}}^{K}-\widetilde{c}\mathscr{S}_{\mathbf{A}}^{K}|_{1,K}\leq
ch_{K}^{p}|u_{0}|_{p+1,K},$ (41)

where $c$ is a positive constant depending on $p$ and on $\gamma$, being
$\gamma$ introduced in (A0)-(A2), but is independent of $h$ and $u$.

Bound (41) is a consequence of the smoothness of $u_{0}$ on all
$K\in\mathscr{T}_{n}^{1}$ and standard polynomial best approximation
estimates; see, e.g., [20]. Collecting the local estimates (40) and (41) and
summing up over all the elements, we get the assertion. ∎

The name of the game in Lemma 4.2 is that the singular part of the solution is
approximated by the singular function in the virtual element spaces on the
elements close or containing the singular vertex.

###  4.3 Best interpolation estimates

Here, we show how to estimate from above the term
$|u-u_{I}|_{1,\mathscr{T}_{n}}$ on the right-hand side of (38) for a specific
choice of $u_{I}$ in $V_{n}$. In particular, we prove an upper bound on the
best interpolation error in nonconforming enriched virtual element spaces in
terms of a constant times an enriched polynomial best approximation term.

######  Lemma 4.3.

Let $u$ be any function in $H^{1}(\Omega)$. Then, for all $u_{\pi}$ piecewise
in $\widetilde{\mathbb{P}}_{p}(K)$, there exists $u_{I}\in V_{n}$ such that

$|u-u_{I}|_{1,\mathscr{T}_{n}}\leq|u-u_{\pi}|_{1,\mathscr{T}_{n}}.$

###### Proof.

The proof follows along the same lines as that of [34, Proposition 3.8]. For
the sake of completeness, we provide some details.

We define $u_{I}\in V_{n}$ by imposing the same degrees of freedom as $u$.
More precisely, set

$\begin{split}&\int_{K}(u-u_{I})q_{p-2}=0\quad\quad\forall
q_{p-2}\in\mathbb{P}_{p-2}(K),\quad\forall K\in\mathscr{T}_{n},\\\
&\int_{e}(u-u_{I})\widetilde{q}_{p-1}^{e}=0\quad\quad\forall\widetilde{q}_{p-1}^{e}\in\widetilde{\mathbb{P}}_{p-1}(e),\quad\forall
e\in\mathscr{E}_{n}.\\\ \end{split}$ (42)

Recall from the definition of the local virtual element spaces in (21) that,
for all $e\in\mathscr{E}_{n}$ and $K\in\mathscr{T}_{n}$,

$\mathbf{n}_{K}\cdot\nabla
u_{I}\in\widetilde{\mathbb{P}}_{p-1}(e),\quad\quad\Delta
u_{I}\in\mathbb{P}_{p-2}(K).$ (43)

We deduce

$\begin{split}|u-u_{I}|^{2}_{1,K}&=\int_{K}\nabla(u-u_{I})\cdot\nabla(u-u_{I})\\\
&=\int_{K}-\Delta(u-u_{I})(u-u_{I})+\int_{\partial
K}\mathbf{n}_{K}\cdot\nabla(u-u_{I})(u-u_{I})\\\
&\overset{\eqref{fixing:dofs},\,\eqref{nice:property}}{=}\int_{K}-\Delta(u-u_{\pi})(u-u_{I})+\int_{\partial
K}\mathbf{n}_{K}\cdot\nabla(u-u_{\pi})(u-u_{I})\\\
&=\int_{K}\nabla(u-u_{I})\cdot\nabla(u-u_{\pi})\leq|u-u_{\pi}|_{1,K}|u-u_{I}|_{1,K}.\\\
\end{split}$ (44)

The assertion follows dividing both sides by $|u-u_{I}|_{1,K}$ and summing
over all the elements. ∎

As a consequence, we have the following best interpolation result in
nonconforming enriched virtual element spaces.

######  Proposition 4.4.

Let $u$ be the solution to problem (11), $u_{0}$ be as in (10), and
assumptions (A0)-(A2) be valid. Then, there exists $u_{I}\in V_{n}$ such that

$|u-u_{I}|_{1,\mathscr{T}_{n}}\leq\,ch^{p}\left\\{\left(\sum_{K\in\mathscr{T}_{n}^{1}}|u_{0}|_{p+1,K}^{2}\right)^{\frac{1}{2}}+\left(\sum_{K\in\mathscr{T}_{n}^{2}\cup\mathscr{T}_{n}^{3}}|u|_{p+1,K}^{2}\right)^{\frac{1}{2}}\right\\},$

where $c$ is exactly the same constant appearing in the bound of Lemma 4.2.

###### Proof.

Combine Lemmas 4.2 and 4.3. ∎

###  4.4 Bound on the variational crime due to the right-hand side

Here, we show an upper bound on

$\sup_{v_{n}\in V_{n}}\frac{\langle
f,v_{n}\rangle-(f,v_{n})_{0,\Omega}}{|v_{n}|_{1,\mathscr{T}_{n}}},$

i.e., on the term representing the variational crime perpetrated in the
discretization of the right-hand side in (1).

######  Lemma 4.5.

Given $p\in\mathbb{N}$, let $f\in H^{p-1}(\Omega)$. Under assumptions
(A0)-(A2), the following bound is valid:

$\sup_{v_{n}\in V_{n}}\frac{\langle
f,v_{n}\rangle-(f,v_{n})_{0,\Omega}}{|v_{n}|_{1,\mathscr{T}_{n}}}\leq
ch^{p}\|f\|_{p-1,\Omega},$

where $c$ is a positive constant depending on $p$ and on $\gamma$, being
$\gamma$ introduced in (A0)-(A2).

###### Proof.

The proof is exactly the same as in the nonenriched VE conforming setting: no
special functions are used in the approximation of the right-hand side; see
[12, Section 4.7] for more details. ∎

###  4.5 Bound on the variational crime due to the nonconformity

Here, we prove an upper bound on the term

$\sup_{v_{n}\in
V_{n}}\frac{\mathscr{N}_{n}(u,v_{n})}{|v_{n}|_{1,\mathscr{T}_{n}}},$

i.e., the term representing the variational crime perpetrated when imposing
the nonconformity of trial and test spaces.

######  Lemma 4.6.

Let $u$ be the solution to problem (11), $u_{0}$ be as in (10),
$\mathscr{N}_{n}$ be defined in (37), and assumptions (A0)-(A2) be valid.
Then, we have

$\frac{\mathscr{N}_{n}(u,v_{n})}{|v_{n}|_{1,\mathscr{T}_{n}}}\leq
ch^{p}\left\\{\left(\sum_{K\in\mathscr{T}_{n}^{1}}|u_{0}|^{2}_{p,K}\right)^{\frac{1}{2}}+\left(\sum_{K\in\mathscr{T}_{n}^{2}\cup\mathscr{T}_{n}^{3}}|u|^{2}_{p,K}\right)^{\frac{1}{2}}\right\\},$

where $c$ is a positive constant depending on $p$ and on $\gamma$, being
$\gamma$ introduced in (A0)-(A2).

###### Proof.

We prove the bound edge by edge. Without loss of generality, we assume that
$e\in\mathscr{E}_{n}^{I}$, for the case $e\in\mathscr{E}_{n}^{B}$ can be
treated analogously, and $e\in\mathscr{E}_{n}^{1}$, for the case
$e\in\mathscr{E}_{n}^{2}$ follows as in [5, Lemma 4.1].

Let $K^{+}$ and $K^{-}$ be the two elements sharing edge $e$. We write

$\int_{e}\nabla u\cdot\llbracket
v_{n}\rrbracket_{e}=\int_{e}\mathbf{n}_{e}\cdot\nabla
u(v_{n}{}_{|K^{+}}-v_{n}{}_{|K^{-}}).$

Denote the $L^{2}(e)$ projector onto constant functions on $e$ by
$\Pi^{0,e}_{0}$. Assumption (10) (and notably the analiticity of $u_{0}$), the
definition of the nonconforming enriched virtual element space in (31), Remark
1, and the properties of orthogonal projectors entail

$\begin{split}\int_{e}\nabla u\cdot\llbracket v_{n}\rrbracket&=\int_{e}\nabla
u_{0}\cdot\llbracket
v_{n}\rrbracket+\int_{e}\nabla\mathscr{S}_{\mathbf{A}}\cdot\llbracket
v_{n}\rrbracket=\int_{e}\nabla u_{0}\cdot\llbracket v_{n}\rrbracket\\\
&=\int_{e}(\mathbf{n}_{e}\cdot\nabla
u_{0}-\mathbf{n}_{e}\cdot\boldsymbol{\Pi}^{0,e}_{p-1}\nabla
u_{0})(v_{n}{}_{|K^{+}}-v_{n}{}_{|K^{-}}-\Pi^{0,e}_{0}(v_{n}{}_{|K^{+}}-v_{n}{}_{|K^{-}}))\\\
&\leq\|\nabla u_{0}-\boldsymbol{\Pi}^{0,e}_{p-1}\nabla
u_{0}\|_{0,e}\|v_{n}{}_{|K^{+}}-v_{n}{}_{|K^{-}}-\Pi^{0,e}_{0}(v_{n}{}_{|K^{+}}-v_{n}{}_{|K^{-}})\|_{0,e}.\end{split}$
(45)

We estimate the two terms on the right-hand side of (45) separately. We begin
with the first one:

$\|\nabla u_{0}-\boldsymbol{\Pi}^{0,e}_{p-1}\nabla u_{0}\|_{0,e}\leq\|\nabla
u_{0}-\boldsymbol{\Pi}^{\nabla,K}_{p-1}\nabla u_{0}\|_{0,e}.$

Apply the trace inequality and the Poincaré-Wirtinger inequality [19, equation
(1.2)], which is valid due to the fact that each component of $\nabla
u_{0}-\boldsymbol{\Pi}^{\nabla,K}_{p-1}\nabla u_{0}$ has zero average on
$\partial K$ by definition of $\boldsymbol{\Pi}^{\nabla,K}_{p-1}$, in addition
to assumption (A2), and get

$\|\nabla u_{0}-\boldsymbol{\Pi}^{0,e}_{p-1}\nabla u_{0}\|_{0,e}\lesssim
h_{K}^{\frac{1}{2}}|\nabla u_{0}-\boldsymbol{\Pi}^{\nabla,K}_{p-1}\nabla
u_{0}|_{1,K}.$

Use standard polynomial approximation theory [20] and assumption (A1) to
arrive at

$\|\nabla u_{0}-\boldsymbol{\Pi}^{0,e}_{p-1}\nabla u_{0}\|_{0,e}\lesssim
h_{K}^{p-\frac{1}{2}}|u_{0}|_{p,K}$ (46)

Focus now on the second term on the right-hand side of (45). As proven in [5,
Lemma 4.1],

$\|v_{n}{}_{|K^{+}}-v_{n}{}_{|K^{-}}-\Pi^{0,e}_{0}(v_{n}{}_{|K^{+}}-v_{n}{}_{|K^{-}})\|_{0,e}\lesssim
h^{\frac{1}{2}}|v_{n}|_{1,K^{+}\cup K^{-}}.$ (47)

The assertion follows combining (45), (46), and (47), and summing over all the
edges. ∎

###  4.6 Stabilizations

Here, we exhibit explicit choices of stabilization $S^{K}(\cdot,\cdot)$
introduced in (32) and discuss their properties. More precisely, we exhibit a
theoretical stabilization, for which we are able to prove the bounds in (32)
explicitly, assuming that $\mathscr{S}_{\mathbf{A}}\in
H^{\frac{3}{2}+\varepsilon}(K)$ with $\varepsilon>0$. Moreover, we introduce a
practical stabilization, which we shall widely employ in the numerical
experiments in Section 5.

##### A theoretical stabilization.

For every $K\in\mathscr{T}_{n}$, define

$S^{K}_{T}(u_{n},v_{n})=h_{K}^{-2}(\Pi^{0,K}_{p-2}u_{n},\Pi^{0,K}_{p-2}v_{n})_{0,K}+h_{K}^{-1}\sum_{e\in\mathscr{E}^{K}}(\widetilde{\Pi}^{0,e}_{p-1}u_{n},\widetilde{\Pi}^{0,e}_{p-1}v_{n})_{0,e}\quad\forall
u_{n},\,v_{n}\in V_{n}(K).$ (48)

Recall that projector $\Pi^{0,K}_{p-2}$ is defined in (29), whereas projector
$\widetilde{\Pi}^{0,e}_{p-1}$ is defined in (28).

######  Proposition 4.7.

Let assumptions (A0)-(A2) be valid and $\mathscr{S}_{\mathbf{A}}\in
H^{\frac{3}{2}+\varepsilon}(K)$ with $\varepsilon>0$. Then, stabilization
$S^{K}_{T}(\cdot,\cdot)$ in $\eqref{theoretical:stabilization}$ satisfies
f(32).

###### Proof.

First, we show the lower bound in (32). For every
$v_{n}\in\ker(\widetilde{\Pi}^{\nabla,K}_{p})$, using the definition of the
local enriched spaces $V_{n}(K)$, we write

$|v_{n}|^{2}_{1,K}=\int_{K}\nabla v_{n}\cdot\nabla v_{n}=\int_{K}-\Delta
v_{n}\Pi^{0,K}_{p-2}v_{n}+\int_{\partial K}(\mathbf{n}_{K}\cdot\nabla
v_{n})\,\widetilde{\Pi}^{0,e}_{p-1}v_{n}.$

Use the Cauchy-Schwarz inequality to get

$|v_{n}|_{1,K}\leq\|\Delta
v_{n}\|_{0,K}\|\Pi^{0,K}_{p-2}v_{n}\|_{0,K}+\|\mathbf{n}_{K}\cdot\nabla
v_{n}\|_{0,\partial K}\|\widetilde{\Pi}^{0,e}_{p-1}v_{n}\|_{0,\partial K}.$
(49)

Using the hypothesis $\mathscr{S}_{\mathbf{A}}\in
H^{\frac{3}{2}+\varepsilon}(K)$ with $\varepsilon>0$, equivalence of norms in
finite dimensional spaces and a scaling argument give

$\|\mathbf{n}_{e}\cdot\nabla v_{n}\|_{0,\partial K}\lesssim
h_{K}^{-\frac{1}{2}}\|\mathbf{n}_{e}\cdot\nabla v_{n}\|_{-\frac{1}{2},\partial
K}.$ (50)

In order to prove (50), we have to be careful. In fact, space
$\widetilde{\mathbb{P}}_{p-1}(e)$ is not invariant by scaling. For each chosen
element in the mesh, this space scales to a different space. Then, we have to
deal with norm equivalences in an infinite number of finite-dimensional
spaces. However, equivalence of norms in finite dimensional spaces is still
valid and can be proven as follows.

The norms appearing in (50) can be seen as continuous functions of the
vertices of element $K$. Assume to fix the number of vertices of $K$ and
denote this number by $N_{V}^{K}$. Thanks to assumptions (A0)–(A1)–(A2), the
set $\mathscr{K}$ of admissible $N_{V}^{K}$ vertices is a bounded, closed set
in $\mathbb{R}^{N_{V}^{K}}$. Thus, the hidden constant is equal to the maximum
amongst all the constant for all possible configurations of the $N_{V}^{K}$
vertices in $\mathscr{K}$. Such maximum is finite, thanks to Weierstraß’
theorem.

Assumptions (A0)–(A1)–(A2) imply the existence of a uniform maximum number of
vertices. Thence, we take the maximum of the above constant for all possible
number of vertices. Importantly, this argument can be applied on all norms on
the boundary $\partial K$.

We are left with the proof of the scaling in (50). Again, this is not
immediate, since $\|\cdot\|_{-\frac{1}{2},\partial K}$ is the dual norm of
$\|\cdot\|_{\frac{1}{2},\partial K}$ containing an $L^{2}$-norm and an
$H^{\frac{1}{2}}$-seminorm components, which have different scalings. Denote
the piecewise quadratic bubble over $\partial K$ by $b_{\partial K}$ and
observe that, for any $\widetilde{q}_{p-1}$ piecewise in
$\widetilde{\mathbb{P}}_{p-1}(e)$,

$\|\widetilde{q}_{p-1}\|_{-\frac{1}{2},\partial K}=\sup_{\psi\in
H^{\frac{1}{2}}(\partial K)\setminus
0}\frac{(\widetilde{q}_{p-1},\psi)_{0,\partial
K}}{\|\psi\|_{\frac{1}{2},\partial
K}}\geq\frac{(\widetilde{q}_{p-1},\widetilde{q}_{p-1}b_{\partial
K})_{0,\partial K}}{\|\widetilde{q}_{p-1}b_{\partial
K}\|_{\frac{1}{2},\partial K}}=\frac{\|\widetilde{q}_{p-1}b_{\partial
K}^{\frac{1}{2}}\|_{0,\partial K}^{2}}{\|\widetilde{q}_{p-1}b_{\partial
K}\|_{\frac{1}{2},\partial K}}.$ (51)

Using equivalence of norms in finite dimensional spaces valid due to the above
arguments, for all edges $e\in\mathscr{E}^{K}$, we have

$\|\widetilde{q}_{p-1}b_{\partial
K}\|_{0,e}\leq\|\widetilde{q}_{p-1}b_{\partial
K}^{\frac{1}{2}}\|_{0,e},\quad\quad|\widetilde{q}_{p-1}b_{\partial
K}|_{1,e}\lesssim h_{e}^{-1}\|\widetilde{q}_{p-1}b_{\partial
K}^{\frac{1}{2}}\|_{0,e}\approx h_{K}^{-1}\|\widetilde{q}_{p-1}b_{\partial
K}^{\frac{1}{2}}\|_{0,e}.$

Summing over all edges and using interpolation theory yield

$\|\widetilde{q}_{p-1}b_{\partial K}\|_{\frac{1}{2},e}\lesssim
h_{K}^{-\frac{1}{2}}\|\widetilde{q}_{p-1}b_{\partial
K}^{\frac{1}{2}}\|_{0,e}.$

Inserting this in (51) gives

$\|\widetilde{q}_{p-1}\|_{-\frac{1}{2},\partial K}\gtrsim
h_{K}^{\frac{1}{2}}\|\widetilde{q}_{p-1}b_{\partial
K}^{\frac{1}{2}}\|_{0,\partial K}.$

Another equivalence equivalence of norms in finite dimensional spaces
eventually entails

$\|\widetilde{q}_{p-1}\|_{-\frac{1}{2},\partial K}\gtrsim
h_{K}^{\frac{1}{2}}\|\widetilde{q}_{p-1}\|_{0,\partial K},$

which is (50), picking $\widetilde{q}_{p-1}{}_{|e}=\mathbf{n}_{e}\cdot\nabla
v_{n}{}_{|e}$ on each edge $e\in\mathscr{E}^{K}$.

Recall also the following inverse inequality is valid; see [22, Lemma 10] and
[14, Theorem 2]:

$\|\Delta v_{n}\|_{0,K}\lesssim h_{K}^{-1}|v_{n}|_{1,K}.$ (52)

The hidden constant depends on the order of accuracy $p$ and on the shape of
element $K$. We point out that inverse estimate (52) was proven in the context
of conforming virtual elements. However, a careful inspection of [22, Lemma
10] and [14, Theorem 2] reveals that it is valid also in the nonconforming
setting with no modifications whatsoever in the proof.

The Neumann trace inequality is valid as well:

$\|\mathbf{n}_{K}\cdot\nabla v_{n}\|_{-\frac{1}{2},\partial
K}\lesssim|v_{n}|_{1,K}+h_{K}\|\Delta
v_{n}\|_{0,K}\overset{\eqref{inverse:VEM}}{\lesssim}|v_{n}|_{1,K}.$ (53)

The Neumann trace inequality is valid not only for polynomials or functions in
virtual element spaces, but for $H^{1}(K)$ functions with Laplacian in
$L^{2}(K)$; see, e.g., [40, Theorem A.33].

Collecting (49), (52), and (53) leads to

$|v_{n}|_{1,K}\lesssim
h_{K}^{-1}\|\Pi^{0,K}_{p-2}v_{n}\|_{0,K}+h_{K}^{-\frac{1}{2}}\|\widetilde{\Pi}^{0,e}_{p-1}v_{n}\|_{0,\partial
K},$

which is the lower bound.

Next, we show the upper bound in (32). We estimate from above the two terms on
the right-hand side of the following identity:

$S^{K}_{T}(v_{n},v_{n})=h_{K}^{-2}\|\Pi^{0,K}_{p-2}v_{n}\|^{2}_{0,K}+h_{K}^{-1}\|\widetilde{\Pi}^{0,e}_{p-1}v_{n}\|^{2}_{0,\partial
K}.$

As for the first term, we use the stability of orthogonal projections and the
Poincaré-Wirtinger inequality [19, equation (1.2)]:

$h_{K}^{-2}\|\Pi^{0,K}_{p-2}v_{n}\|^{2}_{0,K}\lesssim|v_{n}|^{2}_{1,K}.$

As for the second term, use the stability of orthogonal projections, the trace
inequality, the Poincaré-Wirtinger inequality [19, equation (1.2)] again, and
assumption (A2):

$h_{K}^{-1}\|\widetilde{\Pi}^{0,e}_{p-1}v_{n}\|_{0,\partial
K}^{2}\lesssim|v_{n}|^{2}_{1,K}.$

We can apply the Poincaré-Wirtinger inequality [19, equation (1.2)], because
$v_{n}$ belongs to $\ker(\widetilde{\Pi}^{\nabla,K}_{p})$, whence $v_{n}$ has
zero average on $\partial K$.

This concludes the proof. ∎

##### A practical stabilization.

Stabilization $S^{K}_{T}(\cdot,\cdot)$ introduced in (48) is computable in
terms of the degrees of freedom (22) and (23). Notwithstanding, it requires a
certain amount of work to implement.

Thence, we suggest to use the following practical stabilization, which is
defined on the local canonical basis (25) as follows: for all
$K\in\mathscr{T}_{n}$,

$S^{K}_{P}(\varphi_{i},\varphi_{j})=\max(1,a^{K}(\widetilde{\Pi}^{\nabla,K}_{p}\varphi_{i},\widetilde{\Pi}^{\nabla,K}_{p}\varphi_{j}))\quad\quad\forall
i,\,j=1,\dots,\dim(V_{n}(K)).$ (54)

Originally, such stabilization was introduced for the nonenriched 3D VEM in
[15], and its performance was analyzed in the 2D case in [33]. To the best of
our knowledge, such stabilization is amongst the most robust from the
numerical standpoint in the literature. Roughly speaking, this stabilization
keeps trace somehow of the true energy of the basis functions.

In Section 5 below, we perform the numerical experiments using stabilization
$S^{K}_{P}(\cdot,\cdot)$ in (54) in almost all experiments, for it is easier
to implement. We shall compare the performance of the method employing the two
stabilizations in Section 5.1.5 below.

######  Remark 4.

The practical stabilization in (54) is a weighted version of the original VEM
stabilization in [12]. It can be checked that an orthonormalization of the
nonenriched bulk polynomials in (22) and the enriched edge polynomials in (23)
leads to the equivalence of the two stabilizations.

###  4.7 Convergence of the method

In this section, we collect all the a priori bounds hitherto proven and show a
convergence result for the $h$-version of method (35).

######  Theorem 4.8.

Let $u$ and $u_{n}$ be the solutions to (11) and (35), $u_{0}$ be as in (10),
and assumptions (A0)-(A2) be valid. Then, we have the following a priori
$h$-convergence result:

$|u-u_{n}|_{1,\mathscr{T}_{n}}\leq
ch^{p}\left\\{\left(\sum_{K\in\mathscr{T}_{n}^{1}}\|u_{0}\|_{p+1,K}^{2}+\sum_{K\in\mathscr{T}_{n}^{2}\cup\mathscr{T}_{n}^{3}}\|u\|_{p+1,K}^{2}\right)^{\frac{1}{2}}+\|f\|_{p-1,\Omega}\right\\},$
(55)

where $c$ is a positive constant independent of $h$ and $u$, but which
possibly depends on the order of accuracy $p$, parameter $\gamma$ introduced
in (A0)-(A2), and the singular function $\mathscr{S}_{\mathbf{A}}^{K}$.

###### Proof.

Combine Theorem 4.1, Lemma 4.2, Proposition 4.4, Lemma 4.5, Lemma 4.6, and
Proposition 4.7. ∎

######  Remark 5.

In Theorem 4.8, we proved an upper bound on the error
$|u-u_{n}|_{1,\mathscr{T}_{n}}$, which is not computable explicitly. In
Section 5 below, we present a computable quantity based on enriched
projections, and prove that it scales as the exact error.

######  Remark 6.

Recall that elements $K$ in the first layer $\mathscr{T}_{n}^{1}$ are such
that

  * (i)

either $K$ abuts or contains $\mathbf{A}$;

  * (ii)

or $K$ is not too far apart from $\mathbf{A}$.

In particular,
$\text{dist}(\mathbf{A},\mathbf{x}_{K})\leq\widetilde{\gamma}\operatorname{diam}(\Omega)$,
where $\widetilde{\gamma}\in(0,1)$ is a given parameter. The reason why
$\mathscr{T}_{n}^{1}$ does not consist of elements abutting $\mathbf{A}$ only
is that a uniform $h$-refinement implies

$|\Omega|-\sum_{K\in\mathscr{T}_{n}^{2}\cup\mathscr{T}_{n}^{3}}|K|\rightarrow
0\quad\quad\text{as }h\rightarrow 0.$

This has to be avoided, for otherwise the second term on the right-hand side
of (55) would blow up.

###  4.8 Generalizations

Here, we discuss various generalizations of method (35).

##### Multiple singularities.

The first generalization is when expansion (10) is substituted by the general
case (9). We discuss two generalizations.

On the one hand, we can associate more singular functions at one singular
vertex/tip. The extension from the unique singular function is rather
straightforward: it suffices to add more singular functions to the local
spaces $V_{n}(K)$, for all $K\in\mathscr{T}_{n}$ close to $\mathbf{A}$.

On the other hand, we can consider more than one singular vertex/tip. In this
case, we define the first layer $\mathscr{T}_{n}^{1}$ as the set of all
elements close to one or more singular vertices-tips. Thus, the local spaces
are enriched with special functions having singularities at different points.

##### The 3D case.

In 3D, the singularities arising from the geometry of the domain are different
from those in the 2D case. In particular, vertex and edge singularities have
to be dealt with.

The advantage of employing nonconforming virtual element spaces over
conforming ones is that the definition of the local spaces is the same as in
the two dimensional case: fiven $\mathscr{T}_{n}$ a polyhedral decomposition
of the physical domain, for all $K\in\mathscr{T}_{n}$, we set

$V_{n}(K)=\left\\{v_{n}\in H^{1}(K)\mid\Delta
v_{n}\in\mathbb{P}_{p-2}(K),\;\mathbf{n}\cdot\nabla
v_{n}{}_{|F}\in\widetilde{\mathbb{P}}_{p-1}(F)\;\,\forall\,F\text{ face of
}K\right\\}.$

Here, $\widetilde{\mathbb{P}}_{p-1}(K)$ denotes the space of polynomials of
degree at most $p-1$ over a face $F$, possibly enriched of the normal
derivative of the singular functions that we add to the local spaces.

The construction of the method, as well as its analysis, follows exactly along
the same lines as those of the two dimensional case. Clearly, the big issue
here is to find explicit corner and edge singularities.

##### Other differential operators.

It is possible to generalize our setting to PDEs with more general elliptic
operators.

Let $\mathscr{L}$ be an elliptic differential operator of the second order,
$f$ a smooth datum, and $\gamma_{1}$ a trace operator. Consider the following
problem: find $u$ such that

$\begin{cases}\mathscr{L}u=f&\text{in }\Omega\\\ \gamma_{1}u=0&\text{on
}\partial\Omega.\\\ \end{cases}$ (56)

If $\mathscr{L}$ is an elliptic operator, then the solution to problem (56)
presents some singularities due to the geometry of the domain, and more
specifically at corners, edges, tips of cracks, etc. In case the singular
behaviour is known explicitly, the solution can be decomposed as in (9) into a
combination of a smooth and a singular part.

For simplicity, consider the case of a unique singular function $\mathscr{S}$.
Such $\mathscr{S}$ belongs to the kernel of operator $\mathscr{L}$, i.e.,
$\mathscr{L}(\mathscr{S})=0$. Thus, on every element $K\in\mathscr{T}_{n}$
sufficiently close to the singular corner, tip, edge etc., we define the local
space

$V_{n}(K):=\\{v_{n}\in\text{H}\mid\mathscr{L}v_{n}\in[\mathbb{P}_{p_{1}}(K)]^{\ell_{1}},\,\gamma_{2}(v_{n})_{|e}\in[\mathbb{P}_{p_{2}}(e)]^{\ell_{2}}\oplus\gamma_{2}(S)\,\forall
e\in\mathscr{E}^{K}\\}\quad\quad\text{$\ell_{1}$, $\ell_{2}\in\mathbb{N}$}.$

Here, H and $\gamma_{2}(\cdot)$ denote a suitable Sobolev space and a suitable
trace operator.

We postpone the generalizations presented in this section to future works. In
view of the ill-conditioning haunting extended methods, the above mentioned
generalizations might be not straightforward and require an orthonormalization
procedure as illustrated in Section 5.3 and Appendix B.

##  5 Numerical results

In this section, we present some numerical experiments validating the
theoretical results discussed in Section 4.

##### Errors.

In order to measure the convergence of the method, we cannot use the exact
relative error

$\frac{|u-u_{n}|_{1,\mathscr{T}_{n}}}{|u|_{1,\Omega}}.$ (57)

In fact, $u_{n}$ is not known in closed form, but only through its degrees of
freedom. Therefore, we measure the decay of the following _computable_
relative error: given $\widetilde{\Pi}^{\nabla,K}_{p}$ defined piecewise as
$\widetilde{\Pi}^{\nabla,K}_{p}$ over $\mathscr{T}_{n}$,

$\frac{|u-\widetilde{\Pi}^{\nabla,K}_{p}u_{n}|_{1,\mathscr{T}_{n}}}{|u|_{1,\Omega}}.$
(58)

Errors (57) and (58) have the same convergence rate in terms of the mesh size
$h$. In order to see this, on the one hand, we use the stability of the
$H^{1}$ projection to get

$|u-\widetilde{\Pi}^{\nabla}_{p}u_{n}|_{1,\mathscr{T}_{n}}\leq|u-\widetilde{\Pi}^{\nabla}_{p}u|_{1,\mathscr{T}_{n}}+|\widetilde{\Pi}^{\nabla}_{p}u-\widetilde{\Pi}^{\nabla}_{p}u_{n}|_{1,\mathscr{T}_{n}}\leq|u-\widetilde{\Pi}^{\nabla}_{p}u|_{1,\mathscr{T}_{n}}+|u-u_{n}|_{1,\mathscr{T}_{n}}.$

The convergence of the first term on the right-hand side is provided by Lemma
4.2.

On the other hand, we use Theorem 4.8 to get

$\begin{split}|u-u_{n}|_{1,\mathscr{T}_{n}}\lesssim&\max_{K\in\mathscr{T}_{n}}\alpha_{*}^{-1}(K)\left\\{\sup_{v_{n}\in
V_{n}}\frac{\langle
f,v_{n}\rangle_{n}-(f,v_{n})_{0,\Omega}}{|v_{n}|_{1,\mathscr{T}_{n}}}+\sup_{v_{n}\in
V_{n}}\frac{\mathscr{N}_{n}(u,v_{n})}{|v_{n}|_{1,\mathscr{T}_{n}}}\right.\\\
&\quad\quad\quad\quad\quad\quad\quad\left.+\left(1+\max_{K\in\mathscr{T}_{n}}\alpha^{*}(K)\right)|u-\widetilde{\Pi}^{\nabla}_{p}u_{n}|_{1,\mathscr{T}_{n}}\right\\}.\\\
\end{split}$

The convergence of the first two terms on the right-hand side is provided by
Lemmas 4.5 and 4.6.

##### Enriched polynomial basis functions.

We introduce the basis functions for the nonenriched and enriched bulk and
edge polynomial spaces. Consider the natural bijection between $\mathbb{N}$
and $\mathbb{N}_{0}\times\mathbb{N}_{0}$ given by

$1\leftrightarrow(0,0),\quad 2\leftrightarrow(1,0),\quad
3\leftrightarrow(0,1),\quad 4\leftrightarrow(2,0),\quad
5\leftrightarrow(1,1),\quad 6\leftrightarrow(0,2)\quad\dots$ (59)

Henceforth, given a positive integer scalar $\alpha$, we denote the
corresponding vector through the above bijection by
$\boldsymbol{\alpha}=(\alpha_{1},\alpha_{2})$.

Let $\mathbf{x}_{K}=(x_{K},y_{K})$ be the centroid of $K$, for all
$K\in\mathscr{T}_{n}$. For all $K\in\mathscr{T}_{n}$, as for space
$\mathbb{P}_{p}(K)$, we consider the basis elements defined as scaled and
centred monomials

$m_{\alpha}^{K}(\mathbf{x}):=\left(\frac{x-x_{K}}{h_{K}}\right)^{\alpha_{1}}\left(\frac{y-y_{K}}{h_{K}}\right)^{\alpha_{2}}\quad\quad\forall\alpha=1,\dots,\dim(\mathbb{P}_{p}(K)).$
(60)

As for the bulk-enriched polynomial space $\widetilde{\mathbb{P}}_{p}(K)$, we
consider the monomials in (60) _plus_ the function
$\mathscr{S}_{\mathbf{A}}^{K}$ defined in (17) as basis elements.

Denote the Legendre polynomial of degree $\alpha\in\mathbb{N}$ on $[-1,1]$ by
$\mathbb{L}_{\alpha}(x)$. Given $e\in\mathscr{E}_{n}$, denote the linear
transformation mapping the interval $[-1,1]$ to edge $e$ by $\Phi_{e}$. For
all $e\in\mathscr{E}_{n}$, as for space $\mathbb{P}_{p-1}(e)$, we consider the
basis elements defined as scaled and centred Legendre polynomials

$\mathbb{L}^{e}_{\alpha}(\mathbf{x}):=\mathbb{L}_{\alpha}(\Phi_{e}^{-1}(\mathbf{x}))\quad\quad\forall\alpha=1,\dots,\dim(\mathbb{P}_{p-1}(e)).$
(61)

As for the enriched polynomial space $\widetilde{\mathbb{P}}_{p-1}(e)$, we
consider the Legendre polynomials in (61) _plus_ the function
$\mathbf{n}_{e}\cdot\nabla\mathscr{S}_{\mathbf{A}}^{e}{}_{|e}$, where
$\mathscr{S}_{\mathbf{A}}^{e}$ is defined in (18), as basis elements.

For future convenience, the basis elements for the bulk-enriched polynomial
spaces are denoted by

$\widetilde{m}_{\alpha}^{K}\quad\quad\forall\alpha=1,\dots,\dim(\widetilde{\mathbb{P}}_{p}(K)),$

whereas for the edge-enriched polynomial spaces

$\widetilde{m}_{\alpha}^{e}\quad\quad\forall\alpha=1,\dots,\dim(\widetilde{\mathbb{P}}_{p-1}(e)).$

##### Stabilization.

We employ the stabilization defined in (54) throughout and check the
performance of the method using also the theoretical stabilization (48) in
Section 5.1.5 below.

######  Remark 7.

The nonconforming setting allows for an effective improvement of the condition
number of the final system. In order to get such improvement, we suggest to
orthonormalize the basis functions of the edge and bulk-enriched polynomial
basis functions. This is very much in the spirit of [33]; see Section 5.3 and
Appendix B.

Next, we present the test cases we shall analyse numerically. In all of them,
we consider solutions that are singular only at $(0,0)$.

##### Test case 1.

The first test case is defined on the L-shaped domain

$\Omega_{1}=(-1,1)^{2}\setminus[0,1)\times(-1,0].$ (62)

Let $(r,\theta)$ be the polar coordinates at the re-entrant corner $(0,0)$. We
are interested in the approximation of the exact solution

$u_{1}(x,y)=u_{1}(r,\theta)=\sin(\pi\,x)\sin(\pi\,y)+r^{\frac{2}{3}}\sin\left(\frac{2}{3}\theta\right).$
(63)

The primal formulation of the problem we are interested in is such that we
have: zero Dirichlet boundary conditions on the edges generating the re-
entrant corner; suitable Dirichlet boundary conditions on all the other edges;
right-hand side computed according to (63). The right-hand side of problem
(11) with exact solution $u_{1}$ is smooth since the singular function is
harmonic.

In Figure 7, we plot the geometry with the reference frame and the adopted
boundary conditions.

$(0,0)$$(1,0)$ Figure 7: Geometry with the reference frame and the adopted
boundary conditions for test case 1, with exact solution $u_{1}$ in (63). In
red and blue, we highlight the edges where homogeneous and nonhomogeneous
Dirichlet boundary conditions are imposed.

##### Test case 2.

The second test case is defined on the slit square domain

$\Omega_{2}:=(-1,1)^{2}\setminus\\{(x,0)\in\mathbb{R}^{2}\mid x\geq 0\\}.$

Let $(r,\theta)$, $\theta\in[0,2\pi)$, be the polar coordinates at $(0,0)$. We
are interested in the approximation of the exact solution

$u_{2}(r,\theta)=r^{\frac{1}{2}}\sin\left(\frac{1}{2}\theta\right).$ (64)

The primal formulation of the problem we are interested in is such that we
have: suitable Dirichlet boundary conditions on all edges; zero right-hand
side, since $u_{2}$ is harmonic.

##### Test case 3.

The third test case is defined on the L-shaped domain (62). Let $(r,\theta)$
be the polar coordinates at the re-entrant corner $(0,0)$. We are interested
in the approximation of the exact solution

$u_{3}(x,y)=u_{3}(r,\theta)=r^{\frac{2}{3}}\sin\left(\frac{2}{3}\theta\right).$
(65)

The primal formulation of the problem we are interested in is such that we
have: zero Dirichlet boundary conditions on the edges generating the re-
entrant corner; suitable Dirichlet boundary conditions on all the other edges;
zero right-hand side.

######  Remark 8.

Solutions $u_{1}$, $u_{2}$, and $u_{3}$ have nonhomogeneous boundary
conditions. In order to cope with them, we refer to Remark 3.

##### Meshes.

In the forthcoming numerical experiments, we employ sequences of uniform
Cartesian meshes; see Figure 8 for examples of such meshes for the two test
cases. As for test case 2, we highlight the slit in colour: the couples of
adjacent squares do _not_ share the same edge.

Figure 8: Meshes that we employ in the numerical experiments. _Left panel_ :
test case 1. _Right panel_ : test case 2.

We depict the first mesh of the sequence of meshes for test case $2$ in Figure
9: it consists of a single heptagonal element with two edges having endpoints
sharing coordinates. Thus, we show that the virtual element method works and
is robust also on degenerate polygons. In some cases, we shall also use
Voronoi meshes for the sake of testing the method on general polygonal meshes;
see Section 5.1.5.

Figure 9: First mesh of the sequence of meshes employed for test case $2$. It
consists of a single heptagonal element with two edges having endpoints
sharing coordinates.

##### Layers.

We consider different distributions of the layers; see Section 3.2. We test
the method assuming that all the elements of the mesh belong to the first
layer defined in (15), i.e., we enrich all the local spaces. In other words,
we pick $\widetilde{\gamma}=+\infty$ in (15). Furthermore, we consider the
partially enriched scheme, i.e., we enrich only the elements in a
neighbourhood of the singular point $\mathbf{A}$. In Figure 10, we depict the
two different layerisations for the same Cartesian mesh for test case 1. We
pick $\widetilde{\gamma}=1/10$ in (15). The number of enriched elements depend
on the size of the mesh: the finer is the mesh, the more elements are
enriched. The partially enriched scheme corresponds to the geometric
enrichement in the XFEM. A similar layerisation is valid for test case 2 and
is therefore not shown. Further, we shall compare the new enriched method with
the nonenriched one [5]. In this case, all the elements belong to the third
layer defined in (16).

$\mathbf{A}$

$\mathbf{A}$

Figure 10: Two different layerisations for the same Cartesian mesh on the
L-shaped domain $\Omega_{1}$ in (62). _Left panel_ : all the elements belong
to the enriched layer $\mathscr{T}_{n}^{1}$. _Right panel_ : there are three
layers; we pick $\widetilde{\gamma}=1/10$ in (15). We depict: in blue, the
elements in the first element layer $\mathscr{T}_{n}^{1}$; in green, the
elements in the second element layer $\mathscr{T}_{n}^{2}$; in red, the
elements in the third layer $\mathscr{T}_{n}^{3}$.

###  5.1 Numerical experiments on the L-shaped domain

In this section, we present several numerical experiments for test case 1,
with exact solution $u_{1}$ defined in (63), using the fully enriched, the
partially enriched, and the nonenriched [5] schemes. Notably, we are
interested in the performance of the $h$\- and the $p$-versions of the method,
which are the topic of Sections 5.1.1 and 5.1.2. We investigate two additional
computational aspects in Sections 5.1.3 and 5.1.4. First, on a single mesh, we
tune parameter $\widetilde{\gamma}$ in (15) in order to optimize the error of
the method; next, we provide a heuristic motivation as for why the fully
enriched scheme turns out to be less stable than the partially enriched one.
Eventually, in Section 5.1.5, we compare the performance of the method on
Voronoi meshes and use the theoretical and practical stabilizations.

####  5.1.1 The $h$-version

In Figure 11, we present numerical results for the $h$-version of the method.
We consider solution $u_{1}$ in (63) and degrees of accuracy $p=1$, $2$, and
$3$. We use sequences of Cartesian meshes as those in Figure 8, and compare
the performance of the fully enriched (left) and partially enriched methods
(right). Further, we plot the error of the nonenriched nonconforming VEM of
[5]; see Figure 11 (bottom).

Figure 11: $h$-version of various versions of the method in order to
approximate solution $u_{1}$ in (63). The polynomial order of accuracy of the
method is $p=1$, $2$, and $3$. We employ sequences of uniform Cartesian meshes
as those in Figure 8 (left). _Left panel_ : fully enriched method. _Right
panel_ : partially enriched method ($\widetilde{\gamma}=1/10$). Bottom panel:
nonenriched method.

In Figure 11 (left) and (right), we observe an optimal rate of convergence. In
particular, for all the three degrees of accuracy $p$, the two enriched method
converge with order $p$. On the other hand, the nonenriched method converges
suboptimally, which we can expect from the analysis of [5]; see Figure 11
(bottom).

The fully enriched scheme suffers from a loss a convergence for fine meshes.
We investigate the reasons of this phenomenon in Section 5.1.4 below. Further,
we shall provide a remedy for such loss of accuracy in Section 5.3 and
Appendix B below.

####  5.1.2 The $p$-version

In Figure 12, we investigate the behaviour of the $p$-version of the fully and
partially enriched ($\widetilde{\gamma}=1/10$) versions of the method.
Further, we plot the decay of the error employing the nonenriched
nonconforming VEM of [5]. We consider solution $u_{1}$ in (63) and fix a
uniform Cartesian mesh consisting of $48$ elements; see Figure 10. What we
could expect combining the standard theory for the $p$-version of Galerkin
methods, see, e.g., [9], and the analysis developed in this paper, is that the
enriched method converges exponentially in terms of $p$. On the other hand, we
expect an algebraic rate of convergence for the nonenriched method.

Figure 12: $p$-version of various versions of the method in order to
approximate solution $u_{1}$ in (63). We take a uniform Cartesian mesh
consisting of $48$ elements and employ the fully enriched, the partially
enriched ($\widetilde{\gamma}=1/10$), and the nonenriched methods.

From Figure 12, we observe several facts. The nonenriched version of the
nonconforming VEM converges algebraically, whereas the error computed with the
fully enriched method blows up. This is due to the ill-conditioning of the
final system; see Section 5.1.4 below.

More surprisingly, the partially enriched method presents an exponential pre-
asymptotic behaviour, up to $p=3$, but then the convergence turns to be
algebraic. This is due again to the ill-conditioning. When computing error
(58), round-off errors prevent us to have the correct coefficient in the
expansion of $\widetilde{\Pi}^{\nabla,K}_{p}$ so that we are not able to
eliminate the singularity in the solution.

As a positive note, we observe that the $p$-version of the partially enriched
method performs one order of magnitude better than the nonenriched one. We
shall provide a remedy for such loss of accuracy in Section 5.3 below.

####  5.1.3 On the choice of parameter $\widetilde{\gamma}$ in (15)

In this section, we investigate how the choice of parameter
$\widetilde{\gamma}$ appearing in (15) influences the performance of the
method. To this aim, consider solution $u_{1}$ in (63), fix a uniform
Cartesian mesh with $192$ elements, consider $14$ equispaced values of
$\widetilde{\gamma}$ in $(0.1,1.4)$, and depict the error of the method for
each choice of $\widetilde{\gamma}$ in Figure 13.

Figure 13: Performance of the partially enriched method picking different
choices of $\widetilde{\gamma}$ in (15). In particular, we pick $14$
equispaced values of $\widetilde{\gamma}$ in $(0.1,1.4)$ The exact solution is
$u_{1}$ in (63). We employ a uniform Cartesian mesh with 192 elements.

It turns out that the optimal choice of parameter $\widetilde{\gamma}$ lies in
the range $\widetilde{\gamma}\in(0.1,1)$. For larger choices of
$\widetilde{\gamma}$, the error grows as the ill-conditioning of the system
increases.

Needless to say, the above analysis of the best parameter is valid for the
current test case and ought to be performed for every exact solution.

####  5.1.4 On why the fully enriched scheme is less stable than the partially
enriched one

In Sections 5.1.1 and 5.1.3, we observed that the fully enriched scheme is
less stable than the partially enriched one; see Figure 11. In this section,
we give some heuristic motivations as for the reason why this happens.

The motivation behind the growth of the ill-conditioning for fine meshes is
due to the behaviour of the singular enrichment function
$\mathscr{S}_{\mathbf{A}}^{K}$. On the elements that are close to the singular
vertex $\mathbf{A}$, $\mathscr{S}_{\mathbf{A}}^{K}$ differs from all the
scaled monomials, which span the nonenriched polynomial basis. On the other
hand, on the elements that are far from $\mathbf{A}$, the singular function
$\mathscr{S}_{\mathbf{A}}^{K}$ becomes close to a constant function,
especially for small elements. In a sense, the basis functions of the virtual
element space become close to be linearly dependent, as keeping on refining
the mesh.

In Figure 14, we depict the singular function $\mathscr{S}_{\mathbf{A}}^{K}$
with singular behaviour given by $\alpha=2/3$ and a constant function along
the radial component on two intervals of length $0.03$. The first one (left)
is close to $\mathbf{A}$, the second one (right) is slight far from it.

Figure 14: Behaviour of $\mathscr{S}_{\mathbf{A}}$ and a constant function
along the radial component. _Left panel_ : interval $[0.01,0.04]$. _Right
panel_ : interval $[0.01,0.04]$.

From Figure 14, it is apparent that the singular functions get extremely close
to a constant function, thus leading to an ill-conditioned system. A possible
remedy to this situation could be to use bulkwise, see [33], and edgewise
orthonormalization techniques, see Section 5.3 and Appendix B.

####  5.1.5 Different stabilizations and general polygonal meshes

Here, we present a numerical comparison of the performance of the method
employing the theoretical (48) and practical (54) stabilizations. Moreover, we
consider Voronoi meshes, in order to show the robustness of the proposed
method on more general meshes, albeit our method is new on quadrilateral and
triangular meshes as well.

For the test case with exact solution $u_{1}$ in (63), we run an $h$-version
of the partially enriched method with $\widetilde{\gamma}=1/10$ in (15) and
$p=1$, $2$, and $3$; see Figure 15.

Figure 15: $h$-version of the partially enriched method with
$\widetilde{\gamma}=1/10$ in (15) and $p=1$, $2$, and $3$, employing the
theoretical (48) and the practical (54) stabilizations.

The two stabilizations lead to analogous performance of the method. Yet, the
practical one is slighlty more robust for the higher order case, as could have
been expected [33].

###  5.2 Numerical experiments on the slit domain: the extended patch test

In this section, we verify that the method works also on the test case with
exact solution $u_{2}$ in (64). We investigate the performance of the
$h$-version of the method only, employing sequences of uniform Cartesian
meshes as in Figure 8 (right).

Test case $2$ can be regarded as an extended patch test: the exact solution is
equal, up to constants, to the enrichment function
$\mathscr{S}_{\mathbf{A}}^{K}$. Thence, if we consider the fully enriched
version of the method then the error is zero up to machine precision, thanks
to the enriched consistency of the discrete bilinear forms in (34). In Figure
16, we depict the decay of the error of the method for $p=1$, $2$, and $3$. We
consider the fully enriched version of the method.

Figure 16: $h$-version of the fully enriched version of the method in order to
approximate solution $u_{2}$ in (64). The polynomial order of accuracy of the
method is $p=1$, $2$, and $3$. We employ sequences of uniform Cartesian meshes
as those in Figure 8 (right).

From Figure 16, we realize that the fully enriched method returns an error,
which is zero up to machine precision as expected from the enriched
consistency (34). The growth of the error for this enriched patch test is an
excellent indicator for the growth of the ill-conditioning of the system.
Finally, we have evidence that the VEM works also on elements with internal
cracks: the first element of the sequence of meshes is the heptagon with two
overlapping edges depicted in Figure 9.

###  5.3 Orthonormalization of the enriched edge polynomial basis functions

In the foregoing Sections 5.1 and 5.2, we observed that the method suffer of
ill-conditioning, notably employing fine meshes and the high order case.

In this section, we provide numerical evidence that such ill-conditioning can
be drastically reduced by changing the definition of the edge degrees of
freedom. We postpone to Appendix B the design of the new degrees of freedom
and the implementation details, and focus here on the comparison between the
performance of the two methods. We refer to such modification as the
_orthonormal_ enriched method, whereas that investigated so far goes under the
name of _standard_ enriched method. In few words, the former approach is based
on orthonormalization of the basis of enriched edge polynomials to avoid
situations as those described in Section 5.1.4.

For the extended patch test $u_{3}$ in (65), we consider both the $h$\- and
$p$-version of the _orthonormal_ and _standard_ fully enriched methods. In the
former case, we employ sequences of uniform Cartesian meshes, in the latter a
fixed Voronoi mesh. see Figure 17 (left) and (right). Since this is an
extended patch test, the errors should be zero up to machine precision and
their growth is the _real_ indicator of the ill-conditioning of the system.

Figure 17: Behaviour of the error for solution $u_{3}$ in (65). _Left panel:_
$h$-version of the _orthonormal_ and _standard_ fully enriched method on
sequenced of uniform Cartesian meshes. _Right panel:_ $p$-version of the
_orthonormal_ and _standard_ fully enriched method on a fixed Voronoi mesh.

From Figure 17, it is apparent that the _orthonormal_ enriched method
drastically outperforms the _standard_ one. For fine meshes and high degrees
of accuracy, the former produces a relatively small error, which is even ten
orders smaller than that produced by the latter.

For the sake of completeness, we also investigate the behaviour of the
condition number employing the _orthonormal_ and the _standard_ enriched
methods. We report them in Table 1 and 2 for the $h$\- and $p$-versions.

|  $p=1$ \- std. |  $p=1$ \- orth. |  $p=2$ \- std. |  $p=2$ \- orth. |  $p=3$ \- std. |  $p=3$ \- orth.
---|---|---|---|---|---|---
mesh 1 | 8.07e+02 | 2.09e+01 | 3.32e+03 | 6.00e+01 | 7.90e+03 | 1.57e+03
mesh 2 | 5.18e+05 | 2.74e+01 | 3.22e+07 | 8.27e+01 | 6.73e+09 | 2.03e+03
mesh 3 | 1.06e+08 | 5.86e+01 | 5.15e+10 | 2.04e+02 | 6.21e+14 | 4.99e+03
mesh 4 | 2.09e+10 | 2.05e+02 | 5.14e+13 | 7.58e+02 | 9.01e+19 | 1.84e+04
Table 1: Condition number of the system resulting from the $h$-version of the _orthonormal_ and _standard_ fully enriched method on sequenced of uniform Cartesian meshes. | standard | orthonormal
---|---|---
$p=1$ | 1.74e+07 | 1.41e+01
$p=2$ | 6.20e+11 | 1.01e+02
$p=3$ | 4.31e+15 | 9.39e+03
$p=4$ | 1.57e+19 | 7.53e+05
$p=5$ | 2.47e+24 | 4.90e+07
$p=6$ | 1.49e+25 | 7.00e+14
$p=7$ | 1.66e+28 | 2.72e+11
$p=8$ | 3.91e+30 | 7.63e+12
Table 2: Condition number of the system resulting from the $p$-version of the
_orthonormal_ and _standard_ fully enriched method on sequenced of a fixed
Voronoi mesh.

Also in Tables 1 and 2, the _orthonormal_ enriched method results in much
smaller condition numbers. The orthonormalization procedure detailed in
Appendix B below is particularly effective in the nonconforming setting, and
in general in the context of skeletal methods. An analogous procedure in
partition of unity based methods would result in orthonormal basis functions
with increasing support due to the presence of the partition of unity
functions.

Finally, in Figure 18, we compare the $p$-version of the _orthonormal_ and
_standard_ fully enriched methods with the exact solution $u_{1}$ in (63). We
employ sequences of uniform Cartesian meshes, in the latter a fixed Voronoi
mesh.

Figure 18: Behaviour of the error for solution $u_{1}$ in (63). $p$-version of
the _orthonormal_ and _standard_ fully enriched method on a fixed Voronoi
mesh.

The orthonormalization allows us to recover exponential convergence of the
error in terms of the degree of accuracy $p$.

##  6 Conclusions

We introduced a novel enriched nonconforming virtual element method for the
approximation of solutions to the Laplace problem on polygonal domain. This
has been done in the spirit of the extended Galerkin methods, but with a novel
twist. The nonenriched virtual element spaces are endowed with special
singular functions arising from asymptotic singular expansions at the corners
of the domain. We analyzed the method and presented several numerical results,
including the high-order version of the method, which validate the theoretical
predictions. In Appendix B below, we discuss the implementation details.
Importantly, both the theoretical and practical aspects are extensions of what
is done in the nonenriched nonconforming virtual element methods.

In future works, we plan to investigate some generalizations of this approach:

  * •

multiple singularities;

  * •

the 3D version of the method;

  * •

enriched virtual elements for more general elliptic operators.

Notably, our approach seems to be applicable in several branches of
computational mechanics including modeling of cracking phenomena,
discontinuous media, and highly nonlinear complex materials behaviours.

##### Acknowledgements

L. Mascotto acknowledges the support of the Austrian Science Fund (FWF)
through the project P $33477$.

##  Appendix A Implementation details

Here, we discuss the implementation details of the method. We employ the same
notation as in [13]. As in nonenriched nonconforming finite and virtual
elements, the global stiffness matrix is obtained by assembling the local
ones. Therefore, we show the computation of the local stiffness matrices only.

We focus on the elements $K$ in the first layer $\mathscr{T}_{n}^{1}$ only:
the local stiffness matrices on the elements $K\in\mathscr{T}_{n}^{3}$ are
computed as in [5], whereas it suffices to combine the tools employed for the
other two layers on the elements $K\in\mathscr{T}_{n}^{2}$.

Recall that we are assuming (10), i.e., we enrich the approximation space with
one singular function only. The implementation details are utterly similar in
the case of multiple singularities.

We fix the following notation:

$\begin{split}&n_{p}^{K}=\dim(\mathbb{P}_{p}(K)),\quad
n_{p-1}^{e}=\dim(\mathbb{P}_{p-1}(e)),\\\
&\widetilde{n}^{K}_{p}=\dim(\widetilde{\mathbb{P}}_{p}(K)),\quad\widetilde{n}^{e}_{p-1}=\dim(\widetilde{\mathbb{P}}_{p-1}(e)),\quad\widetilde{N}^{K}=\dim(V_{n}(K)).\\\
\end{split}$

Moreover, we set

$\widetilde{\mathbb{P}}_{p}(K)=\text{span}_{\alpha=1}^{\widetilde{n}^{K}_{p}}\\{\widetilde{m}_{\alpha}^{K}\\},\quad\quad\widetilde{\mathbb{P}}_{p-1}(e)=\text{span}_{\alpha=1}^{\widetilde{n}^{e}_{p-1}}\\{\widetilde{m}_{\alpha}^{e}\\},$

where the functions $\widetilde{m}_{\alpha}^{K}$ and
$\widetilde{m}_{\alpha}^{e}$ are defined in Section 5.

Following [13], the local matrix on element $K\in\mathscr{T}_{n}$ is given by

$\mathbf{A}_{n}^{K}=\mathbf{\Pi}^{*}{}^{T}\mathbf{G}\mathbf{\Pi}^{*}+(\mathbf{I}-\mathbf{\Pi})^{T}\mathbf{S}(\mathbf{I}-\mathbf{\Pi}).$
(66)

We define the various matrices appearing in (66). We begin with

$\mathbf{G}_{\alpha,\beta}=\begin{cases}(\nabla\widetilde{m}_{\beta}^{K},\nabla\widetilde{m}_{\alpha}^{K})_{0,K}&\forall\alpha,\,\beta=2,\dots,\dim(\widetilde{\mathbb{P}}_{p}(K))\\\
0&\text{otherwise}.\\\ \end{cases}$

Matrix $\mathbf{\Pi}^{*}$ is the matrix representation of the expansion of
projector $\widetilde{\Pi}^{\nabla,K}_{p}$ in terms of the basis functions of
$\widetilde{\mathbb{P}}_{p}(K)$:

$\mathbf{\Pi}^{*}=\widetilde{\mathbf{G}}^{-1}\mathbf{B},$

where

$\widetilde{\mathbf{G}}_{\alpha,\beta}=\begin{cases}(\nabla\widetilde{m}_{\beta}^{K},\nabla\widetilde{m}_{\alpha}^{K})_{0,K}&\alpha,\,\beta=2,\dots,\dim(\widetilde{\mathbb{P}}_{p}(K))\\\
\frac{1}{|\partial K|}\int_{\partial
K}\widetilde{m}_{\alpha}^{K}&\text{otherwise},\\\ \end{cases}$ (67)

and

$\mathbf{B}_{\alpha,i}=\begin{cases}(\nabla\varphi_{i},\nabla\widetilde{m}_{\alpha}^{K})&\forall\alpha=2,\dots,\dim(\widetilde{\mathbb{P}}_{p}(K)),\,\forall
i=1,\dots,\dim(V_{n}(K))\\\ \frac{1}{|\partial K|}\int_{\partial
K}\varphi_{i}&\text{otherwise}.\\\ \end{cases}$ (68)

Matrix $\mathbf{\Pi}$ is the matrix representation of the expansion of
projector $\widetilde{\Pi}^{\nabla,K}_{p}$ in terms of the basis functions of
$V_{n}(K)$:

$\mathbf{\Pi}=\mathbf{D}\mathbf{\Pi}^{*},$

where

$\mathbf{D}_{i,\alpha}=\text{dof}_{i}(\widetilde{m}_{\alpha}^{K})\quad\forall
i=1,\dots,\dim(V_{n}(K)),\,\forall\alpha=1,\dots,\dim(\widetilde{\mathbb{P}}_{p}(K)).$

Finally, $\mathbf{S}$ is the matrix representation of the stabilization, i.e.,

$\mathbf{S}_{i,j}=S^{K}_{P}(\varphi_{j},\varphi_{i})\quad\forall
i,\,j=1,\dots,\dim(V_{n}(K)).$

If we employ the stabilization introduced in (54), then matrix $\mathbf{S}$ is
diagonal, with entries given by the maximum between $1$ and the corresponding
diagonal entries of the consistency matrix
$\mathbf{\Pi}^{*}{}^{T}\mathbf{G}\mathbf{\Pi}^{*}$.

We devote the remainder of this section to show how to compute matrices
$\mathbf{G}$, $\widetilde{\mathbf{G}}$, $\mathbf{B}$, and $\mathbf{D}$.

##### The matrices $\mathbf{G}$ and $\widetilde{\mathbf{G}}$.

It suffices to show how to compute matrix $\widetilde{\mathbf{G}}$ in (67).
Define

$\widetilde{\mathbf{G}}^{A}\in\mathbb{R}^{1\times\widetilde{n}^{K}_{p}},\quad\quad\widetilde{\mathbf{G}}^{B}\in\mathbb{R}^{(\widetilde{n}^{K}_{p}-1)\times\widetilde{n}^{K}_{p}},$

as follows. Begin with $\widetilde{\mathbf{G}}^{A}$:

$\widetilde{\mathbf{G}}^{A}_{1,\beta}=\frac{1}{|\partial K|}\int_{\partial
K}\widetilde{m}_{\beta}^{K}\quad\forall\beta=1,\dots,\widetilde{n}^{K}_{p}.$

As for $\widetilde{\mathbf{G}}^{B}$, we set

$\widetilde{\mathbf{G}}^{B}_{\alpha,\beta}=\int_{K}\nabla\widetilde{m}_{\alpha}^{K}\cdot\nabla\widetilde{m}_{\beta}^{K}\quad\forall\alpha=2,\dots,\widetilde{n}^{K}_{p},\quad\beta=1,\dots,\widetilde{n}^{K}_{p}.$

Each entry of $\widetilde{\mathbf{G}}^{A}$ and $\widetilde{\mathbf{G}}^{B}$
can be approximated at any precision employing a sufficiently accurate
quadrature formula.

Matrix
$\widetilde{\mathbf{G}}\in\mathbb{R}^{\widetilde{n}^{K}_{p}\times\widetilde{n}^{K}_{p}}$
is given by

$\begin{bmatrix}\mbox{ $\widetilde{\mathbf{G}}^{A}$}\\\
\mbox{$\widetilde{\mathbf{G}}^{B}$}\\\ \end{bmatrix}.$

##### Matrix $\mathbf{B}$.

Denote the number of edges and vertices of $K$ by $N_{V}$ and define

$\begin{split}&\mathbf{B}^{A}\in\mathbb{R}^{1\times
N_{V}p},\quad\quad\quad\quad\;\,\mathbf{B}^{B}\in\mathbb{R}^{1\times
N_{V}},\quad\quad\quad\quad\;\,\mathbf{B}^{C}\in\mathbb{R}^{1\times
n_{p-2}^{K}},\\\ &\mathbf{B}^{D}\in\mathbb{R}^{(n_{p}^{K}-1)\times
N_{V}p},\quad\quad\mathbf{B}^{E}\in\mathbb{R}^{(n_{p}^{K}-1)\times
N_{V}},\quad\quad\mathbf{B}^{F}\in\mathbb{R}^{(n_{p}^{K}-1)\times
n_{p-2}^{K}},\\\ &\mathbf{B}^{G}\in\mathbb{R}^{1\times
N_{V}p},\quad\quad\quad\quad\;\,\mathbf{B}^{H}\in\mathbb{R}^{1\times
N_{V}},\quad\quad\quad\quad\;\,\mathbf{B}^{I}\in\mathbb{R}^{1\times
n_{p-2}^{K}}.\\\ \end{split}$

Matrix
$\mathbf{B}\in\mathbb{R}^{\widetilde{n}^{K}_{p}\times\widetilde{N}^{K}}$ is
given by

$\mathbf{B}=\begin{bmatrix}\mbox{ $\mathbf{B}^{A}$}&\mbox{
$\mathbf{B}^{B}$}&\mbox{ $\mathbf{B}^{C}$}\\\ \mbox{ $\mathbf{B}^{D}$}&\mbox{
$\mathbf{B}^{E}$}&\mbox{ $\mathbf{B}^{F}$}\\\ \mbox{ $\mathbf{B}^{G}$}&\mbox{
$\mathbf{B}^{H}$}&\mbox{ $\mathbf{B}^{I}$}\\\ \end{bmatrix}.$

In matrix $\mathbf{B}$, the first column represents the contributions due to
the basis elements associated with the nonenriched edge polynomials; the
second with special boundary functions; the third with the nonenriched bulk
polynomials. On the other hand, the first row represents the zero average
constraint; the second the contributions of the nonenriched bulk polynomials;
the third the contributions of the singular bulk function.

Owing to (68) and the definition of the edge degrees of freedom in (23), we
set

$\begin{split}&\mathbf{B}^{B}_{1,i}=0\quad\forall
i=1,\dots,N_{V},\quad\quad\mathbf{B}^{C}_{1,i}=0\quad\forall
i=1,\dots,n_{p-2}^{K}.\end{split}$

Denote the $j$-th edge in the local ordering of $K$ by $e(j)$. Vector
$\mathbf{B}^{A}$ has the entries equal to $\frac{|e(j)|}{|\partial K|}$ in the
$j\,p$-th column, for $j=1,\dots,N_{V}$. Otherwise, it has zero entries.

Next, observe that an integration by parts yields

$(\nabla\widetilde{m}_{\alpha}^{K},\nabla\varphi_{i})_{0,K}=-(\Delta\widetilde{m}_{\alpha}^{K},\varphi_{i})_{0,K}+\sum_{e\in\mathscr{E}^{K}}(\mathbf{n}_{K}\cdot\nabla\widetilde{m}_{\alpha}^{K},\varphi_{i})_{0,e}.$
(69)

If $\varphi_{i}$ is an edge basis function, then the first term on the right-
hand side of (69) vanishes. For all $e\in\mathscr{E}^{K}$, if
$\widetilde{m}_{\alpha}^{K}$ is a monomial, then we expand
$\mathbf{n}_{K}\cdot\nabla\widetilde{m}_{\alpha}^{K}{}_{|e}$ into a linear
combination of scaled Legendre polynomials defined in (61):

$\mathbf{n}_{K}\cdot\nabla\widetilde{m}_{\alpha}^{K}{}_{|e}=\sum_{\beta=0}^{p-1}\lambda_{\beta}\mathbb{L}^{e}_{\beta}.$
(70)

We identify the coefficients $\lambda_{\alpha}$ in expansion (70) as follows:
test (70) with any scaled Legendre polynomial of degree at most $p-1$ and use
the orthogonality property

$(\mathbb{L}^{e}_{\alpha},\mathbb{L}^{e}_{\beta})_{0,e}=\frac{h_{e}}{2\beta+1}\delta_{\alpha,\beta}\quad\forall\alpha,\,\beta=0,\dots,p-1,$
(71)

to get

$\lambda_{\beta}=\frac{2\beta+1}{h_{e}}(\mathbf{n}_{K}\cdot\nabla\widetilde{m}_{\alpha}^{K},\mathbb{L}^{e}_{\beta})_{0,e}\quad\forall\beta=0,\dots,p-1.$

The integral on the right-hand side is computable exactly. Thus, expansion
(70) becomes

$\mathbf{n}_{K}\cdot\nabla\widetilde{m}_{\alpha}^{K}{}_{|e}=\sum_{\beta=0}^{p-1}\frac{2\beta+1}{h_{e}}(\mathbf{n}_{K}\cdot\nabla\widetilde{m}_{\alpha}^{K},\mathbb{L}^{e}_{\beta})_{0,e}\mathbb{L}^{e}_{\beta}.$

Given $e(i)$ the edge where the edge basis function $\varphi_{i}$ has a
nonzero edge moment, we write

$\begin{split}(\nabla\widetilde{m}_{\alpha}^{K},\nabla\varphi_{i})_{0,K}&=(\mathbf{n}_{K}\cdot\nabla\widetilde{m}_{\alpha}^{K}{}_{|e(i)},\varphi_{i})_{0,e(i)}\\\
&=\sum_{\beta=0}^{p-1}\frac{2\beta+1}{h_{e(i)}}(\mathbf{n}_{K}\cdot\nabla\widetilde{m}_{\alpha}^{K}{}_{|e(i)},\mathbb{L}^{e(i)}_{\beta})_{0,e(i)}(\mathbb{L}^{e(i)}_{\beta},\varphi_{i})_{0,e(i)}.\end{split}$

Let $j_{e(i)}$ denote the numbering of $\varphi_{i}$ as a basis function on
edge $e(i).$ Using the definition of the edge degrees of freedom (23), we set

$\begin{split}&\mathbf{B}^{D}_{\alpha,i}=(2(j_{e(i)}-1)+1)(\mathbf{n}_{K}\cdot\nabla\widetilde{m}_{\alpha}^{K}{}_{|e(i)},\mathbb{L}^{e(i)}_{i})_{0,e(i)}\quad\forall\alpha=2,\dots,n_{p}^{K}-1,\,\forall
i=1,\dots,pN_{V},\\\
&\mathbf{B}^{E}_{\alpha,i}=0\quad\forall\alpha=2,\dots,n_{p}^{K}-1,\,\forall
i=1,\dots,N_{V}.\\\ \end{split}$

Next, consider the case of $\varphi_{i}$ being an edge basis function and
$\widetilde{m}_{\alpha}^{K}$ being the special function
$\mathscr{S}_{\mathbf{A}}^{K}$ defined in (17). From (19) and (69), we deduce

$(\nabla\widetilde{m}_{\alpha}^{K},\nabla\varphi_{i})_{0,K}=(\mathbf{n}_{K}\cdot\nabla\mathscr{S}_{\mathbf{A}}^{K},\varphi_{i})_{0,e(i)}=\mathbf{n}_{K}\cdot\mathbf{n}_{e(i)}\left(\frac{h_{e(i)}}{h_{K}}\right)^{\alpha}(\mathbf{n}_{e(i)}\cdot\nabla\mathscr{S}_{\mathbf{A}}^{e(i)},\varphi_{i})_{0,e(i)}.$

Using the definition of the edge degrees of freedom (23), we set

$\begin{split}&\mathbf{B}^{G}_{1,i}=0\quad\forall
i=1,\dots,N_{V}p,\quad\quad\mathbf{B}^{H}_{1,i}=\mathbf{n}_{K}\cdot\mathbf{n}_{e(i)}\left(\frac{h_{e(i)}}{h_{K}}\right)^{\alpha}\quad\forall
i=1,\dots,N_{V}.\\\ \end{split}$

Finally, focus on the case given by $\varphi_{i}$ being a bulk basis function.
Firstly, assume that $\widetilde{m}_{\alpha}^{K}=m_{\alpha}^{K}$ is a
nonenriched polynomial. Given $(x_{K},y_{K})$ the centroid of $K$, the
following splitting is valid:

$\begin{split}\Delta
m_{\alpha}^{K}&=\Delta\left(\left(\frac{x-x_{K}}{h_{K}}\right)^{\alpha_{1}}\left(\frac{y-y_{K}}{h_{K}}\right)^{\alpha_{2}}\right)\\\
&=\frac{1}{h_{K}^{2}}\left(\alpha_{1}(\alpha_{1}-1)\left(\frac{x-x_{K}}{h_{K}}\right)^{\alpha_{1}-2}\left(\frac{y-y_{K}}{h_{K}}\right)^{\alpha_{2}}+\alpha_{2}(\alpha_{2}-1)\left(\frac{x-x_{K}}{h_{K}}\right)^{\alpha_{1}}\left(\frac{y-y_{K}}{h_{K}}\right)^{\alpha_{2}-2}\right).\end{split}$

Let $\alpha$ be associated to $(\alpha_{1},\alpha_{2})$ via bijection (59).
Whenever it makes sense, set $\widetilde{\alpha}_{1}$ and
$\widetilde{\alpha}_{2}$ the natural numbers associated with
$(\alpha_{1}-2,\alpha_{2})$ and $(\alpha_{1}-2,\alpha_{2})$, via the same
bijection (59).

We have

$\begin{split}&(\nabla m_{\alpha}^{K},\nabla\varphi_{i})_{0,K}=-(\Delta
m_{\alpha}^{K},\varphi_{i})_{0,K}\\\ &=\begin{cases}0&\text{if
}\alpha=1,2,3\\\
-\frac{|K|}{h_{K}^{2}}\left((\alpha_{1}-1)\alpha_{1}\frac{1}{|K|}(m_{\widetilde{\alpha}_{1}}^{K},\varphi_{i})_{0,K}+(\alpha_{2}-1)\alpha_{2}\frac{1}{|K|}(m_{\widetilde{\alpha}_{2}}^{K},\varphi_{i})_{0,K}\right)&\text{otherwise}.\end{cases}\end{split}$

In other words, we get

$\begin{split}&\mathbf{B}^{F}_{\alpha,i}=\begin{cases}0&\text{if
}\alpha=1,2,3,\;\forall i=2,\dots,n_{p-2}^{K}\\\
-\frac{|K|}{h_{K}^{2}}\left((\alpha_{1}-1)\alpha_{1}\delta_{\widetilde{\alpha}_{1},i}+(\alpha_{2}-1)\alpha_{2}\delta_{\widetilde{\alpha}_{2},i}\right)&\forall\alpha=4,\dots,n_{p}^{K},\,\forall
i=1,\dots,n_{p-2}^{K},\\\ \end{cases}\\\ &\mathbf{B}^{I}_{1,i}=0\quad\forall
i=1,\dots,n_{p-2}^{K}.\\\ \end{split}$

##### Matrix $\mathbf{D}$.

We introduce

$\mathbf{D}^{A}\in\mathbb{R}^{pN_{V}\times\widetilde{n}^{K}_{p}},\quad\quad\mathbf{D}^{B}\in\mathbb{R}^{N_{V}\times\widetilde{n}^{K}_{p}},\quad\quad\mathbf{D}^{C}\in\mathbb{R}^{n_{p-2}^{K}\times\widetilde{n}^{K}_{p}},$

so that matrix
$\mathbf{D}\in\mathbb{R}^{\widetilde{N}^{K}\times\widetilde{n}^{K}_{p}}$ is
given by

$\begin{bmatrix}\mbox{ $\mathbf{D}^{A}$}\\\ \mbox{ $\mathbf{D}^{B}$}\\\ \mbox{
$\mathbf{D}^{C}$}\\\ \end{bmatrix}.$

Matrix $\mathbf{D}^{A}$ represents the contributions of the nonenriched edge
basis functions; $\mathbf{D}^{B}$ the contributions of the special edge
functions; $\mathbf{D}^{C}$ the contributions of the bulk functions.

Given $\varphi_{i}$ an edge basis function, let $e(i)$ the edge, where
$\varphi_{i}$ has a nonzero moment. If $\varphi_{i}$ is associated with the
nonenriched polynomial moments, then denote the nonzero order moment by
$\beta(i)$. Recalling the definition of the edge moments (23), we set

$\begin{split}&\mathbf{D}^{A}_{i,\alpha}=\frac{1}{|e(i)|}(\mathbb{L}_{\beta(i)}^{e(i)},\widetilde{m}_{\alpha}^{K}{}_{|e(i)})_{0,e(i)}\quad\forall
i=1,\dots,pN_{V},\,\forall\alpha=1,\dots,\widetilde{n}^{K}_{p},\\\
&\mathbf{D}^{B}_{i,\alpha}=(\mathbf{n}_{e}\cdot\nabla\mathscr{S}_{\mathbf{A}}^{e(i)},\widetilde{m}_{\alpha}^{K}{}_{|e(i)})_{0,e(i)}\quad\forall
i=1,\dots,N_{V},\,\forall\alpha=1,\dots,\widetilde{n}^{K}_{p}.\\\ \end{split}$

As for matrix $\mathbf{D}^{C}$, we simply write

$\mathbf{D}^{C}_{i,\alpha}=\frac{1}{|K|}\int_{K}m^{K}_{i}\widetilde{m}_{\alpha}^{K}\quad\forall
i=1,\dots,n_{p-2}^{K},\,\forall\alpha=1,\dots,\widetilde{n}^{K}_{p}.$ (72)

All the entries of the three matrices above can be computed exactly or
approximated at any precision with a sufficiently accurate quadrature formula.

######  Remark 9.

As for the computation of $\mathbf{D}^{C}$ in (72) in the case
$\widetilde{m}_{\alpha}^{K}=\mathscr{S}_{\mathbf{A}}^{K}$, we suggest to use
the following strategy. Given $m^{K}_{i}\in\mathbb{P}_{p-2}(K)$, it is
possible to write

$m^{K}_{i}=\Delta m^{K}_{p,i}$

for some $m^{K}_{p,i}\in\mathbb{P}_{p}(K)$. We provide an explicit
representation of $m^{K}_{p,i}$ in Appendix C.

Then, the integral in (72) can be rewritten using an integration by parts
twice and the fact that $\Delta\mathscr{S}_{\mathbf{A}}=0$ as

$\frac{1}{|K|}\int_{K}m^{K}_{i}\mathscr{S}_{\mathbf{A}}^{K}=\frac{1}{|K|}\int_{K}\Delta
m^{K}_{p,i}\mathscr{S}_{\mathbf{A}}^{K}=\frac{1}{|K|}\left[\int_{\partial
K}\mathbf{n}\cdot\nabla m^{K}_{p,i}\mathscr{S}_{\mathbf{A}}^{K}-\int_{\partial
K}m^{K}_{p,i}\mathbf{n}\cdot\nabla\mathscr{S}_{\mathbf{A}}^{K}\right].$

######  Remark 10.

In the computation of matrix $\mathbf{G}^{B}$, if at least one of the two
terms, say $\widetilde{m}_{\alpha}^{K}$, is the singular enrichment function
$\mathscr{S}_{\mathbf{A}}^{K}$, we reduce the computation of the bulk integral
to the computation of the boundary integral

$\int_{K}\nabla\widetilde{m}_{\alpha}^{K}\cdot\nabla\widetilde{m}_{\beta}^{K}=\int_{K}\nabla\mathscr{S}_{\mathbf{A}}^{K}\cdot\nabla\widetilde{m}_{\beta}^{K}=\int_{\partial
K}\mathbf{n}\cdot\nabla\mathscr{S}_{\mathbf{A}}^{K}\,\widetilde{m}_{\beta}^{K}.$

In the light of this fact, in the computation of matrices $\mathbf{B}$,
$\mathbf{D}$, $\mathbf{G}$, and the boundary conditions, the integrals
involving singular functions are always boundary integrals. This fact is
extremely relevant. Indeed, in order to compute integrals involving singular
functions, we resort to Gauß-Jacobi quadrature formulas; see, e.g., [39,
Section 4.8-1]. By doing so, the singular integrals can be computed up to
machine precision with relatively few quadrature knots, whereas, in order to
achieve the same precision with the standard Gauß integration rule, we would
need to require a disproportionate number of quadrature knots.

######  Remark 11.

The “$G=BD$” test of [13, Remark 3.3] is valid also in the enriched framework.
This is an excellent test to check the correctness of the implementation of
the method. In order to fulfil this test correct, the integrals must be
computed up to machine precision. Notably, we suggest to use suitable
quadrature formulas; see Remark 10.

######  Remark 12.

In view of possible extensions to linear elasticity, it might be of interest
to discuss the approximation of the gradient of the discrete solution in the
elements and on faces. In the bulk of the elements, we can consider
$\nabla\widetilde{\Pi}^{\nabla,K}_{p}u_{n}$, whereas, on an edge $e$, we can
consider, e.g., the average of the energy projection on the neighbouring
elements $K^{+}$ and $K^{-}$:

$\frac{1}{2}(\nabla\widetilde{\Pi}^{\nabla,K}_{p}u_{n}{}_{|K^{+}}\cdot\mathbf{n}_{K^{+}}+\nabla\widetilde{\Pi}^{\nabla,K}_{p}u_{n}{}_{|K^{-}}\cdot\mathbf{n}_{K^{-}}).$

##### Computation of the right-hand side.

Proceed as in [13]: no enrichment affects the right-hand side.

###  A.1 Nonhomogeneous Dirichlet boundary conditions

As for the treatment of nonhomogenous Dirichlet boundary conditions, we
identify the boundary edge degrees of freedom of the discrete and exact
solutions $u_{n}$ and $u$. In other words, for all $e\in\mathscr{E}_{n}^{B}$
such that $e\subset\Gamma_{D}$, we impose the following condition:

$\int_{e}(u_{n}-u)\widetilde{m}_{\alpha}^{e}=0\quad\quad\forall\alpha=1,\dots,\dim(\widetilde{\mathbb{P}}_{p-1}(e)).$

###  A.2 Nonhomogeneous Neumann boundary conditions

Here, we address the implementation aspects for the computation of the Neumann
boundary conditions term (36). In particular, given a canonical basis function
$\varphi_{i}$ associated with a nonzero moment on a Neumann edge
$e\subset\Gamma_{N}$, we describe how to compute

$\int_{e}g_{N}\widetilde{\Pi}^{0,e}_{p-1}\varphi_{i}.$ (73)

As highlighted in Remark 2, for all $e\in\mathscr{E}_{n}$, projector
$\widetilde{\Pi}^{0,e}_{p-1}$ can be computed only under assumption
$\mathscr{S}_{\mathbf{A}}\in H^{\frac{3}{2}+\varepsilon}(K)$ with
$\varepsilon>0$, where $K\in\mathscr{T}_{n}$ is such that
$e\in\mathscr{E}^{K}$. At the end of this section, we show that, under
suitable assumptions on $g_{N}$, we can indeed compute nonhomogenous Neumann
boundary conditions for $\mathscr{S}_{\mathbf{A}}\not\in
H^{\frac{3}{2}+\varepsilon}(K)$, $\varepsilon>0$, as well, with no need
whatsoever of resorting to projector $\widetilde{\Pi}^{0,e}_{p-1}$.

In order to compute enriched edge projections of the basis functions, consider
the expansion

$\widetilde{\Pi}^{0,e}_{p-1}\varphi_{i}=\sum_{\alpha=0}^{p-1}\lambda^{0}_{\alpha}m_{\alpha}^{e}+\lambda^{0}_{p}\mathbf{n}_{e}\cdot\nabla\mathscr{S}_{\mathbf{A}}^{e}.$
(74)

Once we know the coefficients $\lambda^{0}_{\alpha}$, $\alpha=0,\dots,p$, we
are able to approximate the integral in (73) at any precision.

Define matrix $\mathbf{G}^{0,e}\in\mathbb{R}^{(p+1)\times(p+1)}$ and vector
$\mathbf{b}^{0,e}\in\mathbb{R}^{p+1,1}$ as follows:

$\mathbf{G}^{0,e}_{\alpha,\beta}:=\begin{cases}\frac{h_{e}}{2}\frac{2}{2(\beta-1)+1}=\frac{h_{e}}{2(\beta-1)+1}&\text{if
}\alpha=\beta,\;\beta=1,\dots,p\\\ 0&\text{if }\alpha\text{ is not equal to
}\beta,\;\alpha,\beta=1,\dots,p\\\
(\mathbf{n}_{e}\cdot\nabla\mathscr{S}_{\mathbf{A}}^{e},m_{\beta}^{e})_{0,e}&\text{if
}\alpha=p+1\text{ and }\beta=1,\dots,p;\;\beta=p+1\text{ and
}\alpha=1,\dots,p\\\
(\mathbf{n}_{e}\cdot\nabla\mathscr{S}_{\mathbf{A}}^{e},\mathbf{n}_{e}\cdot\nabla\mathscr{S}_{\mathbf{A}}^{e})_{0,e}&\text{if
}\alpha=\beta=p+1.\end{cases}$

Moreover, define vector $\mathbf{b}^{0,e}\in\mathbb{R}^{p+1,1}$ as the $i$-th
column of the diagonal matrix
$\mathbf{B}^{0,e}\in\mathbb{R}^{(p+1)\times(p+1)}$, which is given by

$\mathbf{B}^{0,e}_{\alpha,i}=\begin{cases}h_{e}&\text{if }\alpha=i\leq p\\\
1&\text{if }\alpha=i=p+1.\\\ \end{cases}$

Matrix $\boldsymbol{\Lambda}^{0,e}$ of the coefficients in (74) for the
expansion of the basis function element $\varphi$ is computed solving the
system

$\mathbf{G}^{0,e}\boldsymbol{\Lambda}^{0,e}=\mathbf{b}^{0,e}.$

In order to see this, it suffices to test (74) with the elements in a basis of
$\widetilde{\mathbb{P}}_{p-1}(e)$ and use the orthogonality property of the
Legendre polynomials (71).

The computation of nonhomogenous Neumann boundary conditions can be simplified
and extended to the case of general singular functions
$\mathscr{S}_{\mathbf{A}}^{K}\not\in H^{\frac{3}{2}+\varepsilon}(\Omega)$ with
$\varepsilon>0$. In particular, assume that, given $c\in\mathbb{R}$,

$g_{N}=\mathbf{n}_{\Omega}\cdot\nabla
u=c\,\mathbf{n}_{\Omega}\cdot\nabla\mathscr{S}_{\mathbf{A}}.$ (75)

More generally, we can assume that
$g_{N}{}_{|e}\in\widetilde{\mathbb{P}}_{p-1}(e)$ for all $e\subset\Gamma_{N}$.
We employ the following discretization of the Neumann datum contribution:

$\sum_{e\in\mathscr{E}_{n}^{B},\,e\subset\Gamma_{N}}\int_{e}g_{N}\varphi_{i}.$

Fix $e\in\mathscr{E}_{n}^{B}$ with $e\subset\Gamma_{N}$ and let $\varphi_{i}$
be a basis function, with a nonzero moment on edge $e$. Denote the power of
the singularity of $\mathscr{S}_{\mathbf{A}}$ by $\alpha$. Thanks to
assumption (75), we can write

$\begin{split}\int_{e}g_{N}\varphi_{i}&=c\int_{e}\mathbf{n}_{\Omega}\cdot\nabla\mathscr{S}_{\mathbf{A}}\,\varphi_{i}=ch_{e}^{\alpha}\int_{e}\mathbf{n}_{\Omega}\cdot\nabla\mathscr{S}_{\mathbf{A}}^{e}\,\varphi_{i}=ch_{e}^{\alpha}\int_{e}\mathbf{n}_{\Omega}\cdot\nabla\mathscr{S}_{\mathbf{A}}^{e}\,\varphi_{i}.\end{split}$

Due to the definition of the enriched edge degrees of freedom (23), this
quantity is equal to zero for all basis functions dual to the scaled Legendre
polynomials, whereas it is equal to $h_{e}^{\alpha}$ if $\varphi_{i}$ is dual
to the singular edge function.

##  Appendix B Design and implementation of a robust variant of (35)

In order to mitigate the ill-conditioning observed in Section 5, we discuss
method (35) changing the definition of the edge degrees of freedom on enriched
edges; see Section B.1. The numerical results with this version of the method
are provided in Section 5.1.5. We also provide some implementation details for
such version of the method; see Section B.2.

Importantly, the orthonormalization procedure is based on a “diagonalization”
process of matrix $\mathbf{G}^{0}$ in Appendix A.2. Thence, we need to assume
that the singular function $\mathscr{S}_{\mathbf{A}}\in
H^{\frac{3}{2}+\varepsilon}(\Omega)$ with $\varepsilon>0$.

To the best of our understanding, this procedure is not possible to use in
other contexts, e.g., in the setting of partition of unity methods. Applying
there any orthonormalization whatsoever, in fact, would result in a dramatic
loss of localization of the basis functions.

###  B.1 Orthonormalization of enriched edge polynomials

Consider method (35) with a modified definition of the degrees of freedom on
the enriched edges $e\in\mathscr{E}_{n}^{1}$. More precisely, recall that
$\\{\widetilde{m}_{\alpha}^{e}\\}_{\alpha=0}^{p-1}$ denote the basis of the
enriched edge polynomial space $\widetilde{\mathbb{P}}_{p-1}(e)$ consisting of
the first $p$ scaled Legendre polynomials and
$\mathbf{n}\cdot\nabla\mathscr{S}_{\mathbf{A}}^{e}$. The degrees of freedom on
enriched edges (23) have been defined with respect to such a basis.

As discussed in Section 5.1.4, for small elements and high polynomial degrees,
those basis elements become close to linear dependent. Therefore, we
$L^{2}(e)$-orthonormalize the elements $\widetilde{m}_{\alpha}^{e}$. For
instance, we can use a stable Gram-Schmidt orthonormalization as that
presented in [10, Section 2]. Denote the new $L^{2}(e)$-orthonormal basis
elements by $\\{\overline{m}_{\alpha}^{e}\\}_{\alpha=0}^{p-1}$.

The modified method is based on the same local and global virtual element
spaces, the same bulk degrees of freedom, the same edge degrees of freedom on
nonenriched edges $e\in\mathscr{E}_{n}^{2}$, and edge degrees of freedom on
enriched edges $e\in\mathscr{E}_{n}^{1}$ with respect to the new basis of
$\widetilde{\mathbb{P}}_{p-1}(e)$. Due to the normalization of the basis
functions, the scaling of the edge degrees of freedom is given by

$\frac{1}{h_{e}^{\frac{1}{2}}}\int_{e}v_{n}\widetilde{m}_{\alpha}^{e}.$

Clearly, the analysis of the method is the same as for the original one.

###  B.2 Implementation details

Here, we provide some implementation details for the new setting discussed in
Section B.1. In particular, we explain how to compute the various matrices
needed in the implementation of the method. We denote the matrices computed
with the new method adding a bar on top of their counterparts in Appendix A.
Note that we only modify the matrices associated with enriched edges and
elements. Moreover, for the sake of conciseness, we avoid to discuss the
details of the matrices associated with the elements $K\in\mathscr{T}_{n}^{2}$
and rather focus on those associated with $K\in\mathscr{T}_{n}^{1}$.

Fix $e\in\mathscr{E}_{n}^{1}$ and let
$\mathbf{GS}^{e}\in\mathbb{R}^{(p+1)\times(p+1)}$ be the lower triangular
matrix containing the coefficients obtained via the orthonormalization such
that

$\overline{m}_{\alpha}^{e}=\sum_{\beta=1}^{\alpha}\mathbf{GS}^{e}_{\alpha,\beta}\widetilde{m}_{\beta}^{e}\quad\quad\forall\alpha=0,\dots,p-1.$

Matrix $\mathbf{GS}^{e}$ can be computed, e.g., as in [10, Section 2].

######  Remark 13.

The structure of the Gram-Schmidt orthonormalizing matrix is partially known a
priori. As already mentioned, $\mathbf{GS}^{e}$ is lower triangular. Moreover,
the upper-left block in $\mathbb{R}^{p\times p}$ is diagonal, with diagonal
entry given by $\|m_{\alpha}^{e}\|_{0,e}^{-1}$ for all $\alpha=1,\dots,p$.
This follows from the fact that the scaled Legendre basis we employ is already
orthogonal, albeit it needs to be normalized. Thus, the only “full” row of
$\mathbf{GS}^{e}$ is the last one, because the singular function has no
orthogonality property whatsoever with respect to the scaled Legendre
polynomials. Knowing a priori the structure of $\mathbf{GS}^{e}$ is of extreme
help in the computation of the local matrices below.

In what follows, we employ the same notation as that in Appendix A.

##### Matrix $\overline{\mathbf{G}}^{0,e}$.

For all $e\in\mathscr{E}_{n}^{1}$, matrix
$\overline{\mathbf{G}}^{0,e}\in\mathbb{R}^{(p+1)\times(p+1)}$ is defined as
the identity matrix. This follows from the orthonormality of the elements in
the new basis of the enriched edge polynomial spaces.

##### Matrix $\overline{\mathbf{B}}^{0,e}$.

For all $e\in\mathscr{E}_{n}^{1}$, matrix
$\overline{\mathbf{B}}^{0,e}\in\mathbb{R}^{(p+1)\times(p+1)}$ is defined as
the diagonal matrix with diagonal entries given by $h_{e}^{\frac{1}{2}}$. This
follows from the fact that the elements of the new basis of the enriched edge
polynomial spaces are normalized in $L^{2}(e)$ and therefore the moments need
the proper scaling.

##### Matrix $\overline{\mathbf{G}}$.

Matrix
$\overline{\mathbf{G}}\in\mathbb{R}^{\widetilde{n}^{K}_{p}\times\widetilde{n}^{K}_{p}}$
is equal to matrix $\mathbf{G}$. In fact, it involves only the product of
basis elements of the enriched bulk polynomial space, which has not been
modified in the new setting.

##### Matrix $\overline{\mathbf{B}}$.

Matrix
$\overline{\mathbf{B}}\in\mathbb{R}^{\widetilde{n}^{K}_{p}\times\widetilde{N}^{K}}$
can be split into the four submatrices

$\begin{split}&\overline{\mathbf{B}}^{A,B}\in\mathbb{R}^{1\times
N_{V}(p+1)},\quad\quad\quad\quad\;\overline{\mathbf{B}}^{C}\in\mathbb{R}^{1\times
n_{p-2}^{K}},\\\
&\overline{\mathbf{B}}^{D,E,G,H}\in\mathbb{R}^{n_{p}^{K}\times
N_{V}(p+1)},\quad\quad\overline{\mathbf{B}}^{F,I}\in\mathbb{R}^{n_{p}^{K}\times
n_{p-2}^{K}}.\\\ \end{split}$

In particular, we write

$\overline{\mathbf{B}}=\begin{bmatrix}\mbox{
$\overline{\mathbf{B}}^{A,B}$}&\mbox{ $\overline{\mathbf{B}}^{C}$}\\\ \mbox{
$\overline{\mathbf{B}}^{D,E,G,H}$}&\mbox{ $\overline{\mathbf{B}}^{F,I}$}\\\
\end{bmatrix}.$

Since there is no modification of the bulk degrees of freedom, we have

$\overline{\mathbf{B}}^{C}=\mathbf{B}^{C},\quad\quad\overline{\mathbf{B}}^{F,I}=\begin{bmatrix}\mbox{
$\mathbf{B}^{F}$}\\\ \mbox{ $\mathbf{B}^{I}$}\\\ \end{bmatrix}.$

Next, focus on matrix $\overline{\mathbf{B}}^{A,B}$. Thanks to Remark 13 and
the fact that $\widetilde{m}_{1}^{e}=1$, we have

$\frac{1}{|\partial K|}\int_{\partial
K}\overline{\varphi}_{i}=\frac{(\mathbf{GS}^{e}_{1,1})^{-1}}{|\partial
K|}\int_{\partial
K}\mathbf{GS}^{e}_{1,1}\overline{\varphi}_{i}=\frac{h_{e}^{\frac{1}{2}}}{\mathbf{GS}^{e}_{1,1}\,|\partial
K|}\,\frac{1}{h_{e}^{\frac{1}{2}}}\int_{\partial
K}\overline{m}^{e}_{1}\,\overline{\varphi}_{i}.$

We deduce that $\overline{\mathbf{B}}^{A,B}$ has the entries equal to
$h_{e}^{\frac{1}{2}}/(\mathbf{GS}^{e}_{1,1}\,|\partial K|)$ in the $j\,p$-th
column for all $j=1,\dots,N_{V}$, where $e(j)$ denotes the $j$-th edge in the
local ordering of $K$.

Eventually, we deal with matrix $\overline{\mathbf{B}}^{D,E,G,H}$, i.e., on
the case of $\widetilde{m}_{\alpha}^{K}$ being a bulk enriched basis function
and $\varphi_{i}$ a basis function of the virtual element space, which is dual
to an edge moment. We write

$\int_{\partial
K}\nabla\widetilde{m}_{\alpha}^{K}\cdot\nabla\overline{\varphi}_{i}=\sum_{e\in\mathscr{E}^{K}}\int_{e}(\mathbf{n}_{K}\cdot\nabla\widetilde{m}_{\alpha}^{K})_{|e}\overline{\varphi}_{i}.$
(76)

For all $e\in\mathscr{E}^{K}$, we need to expand
$(\mathbf{n}_{K}\cdot\nabla\widetilde{m}_{\alpha}^{K})_{|e}$ into the
orthonormalized basis elements:

$(\mathbf{n}_{K}\cdot\nabla\widetilde{m}_{\alpha}^{K})_{|e}=\sum_{\beta=1}^{p+1}\overline{\lambda}_{\beta}\overline{m}_{\beta}^{e}.$

Testing the above identity with $\overline{m}_{\gamma}^{e}$, for all
$\gamma=1,\dots,p+1$, gives

$(\mathbf{n}_{K}\cdot\nabla\widetilde{m}_{\alpha}^{K},\overline{m}_{\gamma}^{e})_{0,e}=\overline{\lambda}_{\gamma}.$

Thus, the following decomposition is valid:

$(\mathbf{n}_{K}\cdot\nabla\widetilde{m}_{\alpha}^{K})_{|e}=\sum_{\beta=1}^{p+1}(\mathbf{n}_{K}\cdot\nabla\widetilde{m}_{\alpha}^{K},\overline{m}_{\beta}^{e})_{0,e}\overline{m}_{\beta}^{e}.$

Inserting this into (76) and using the definition of the new degrees of
freedom yield

$\int_{\partial
K}\nabla\widetilde{m}_{\alpha}^{K}\cdot\nabla\overline{\varphi}_{i}=\sum_{e\in\mathscr{E}^{K}}\sum_{\beta=1}^{p+1}h_{e}^{\frac{1}{2}}(\mathbf{n}_{K}\cdot\nabla\widetilde{m}_{\alpha}^{K},\overline{m}_{\beta}^{e})_{0,e}\frac{1}{h_{e}^{\frac{1}{2}}}\int_{e}\overline{m}_{\beta}^{e}\overline{\varphi}_{i}.$

Let $j_{e(i)}$ denote the numbering of $\varphi_{i}$ as a basis function on
edge $e(i)$. Then, we can write, thanks to the definition of the enriched edge
degrees of freedom (23),

$\begin{split}\overline{\mathbf{B}}^{D,E,G,H}_{\alpha,i}=\int_{\partial
K}\nabla\widetilde{m}_{\alpha}^{K}\cdot\nabla\overline{\varphi}_{i}&=\sum_{\beta=1}^{p+1}h_{e(i)}^{\frac{1}{2}}(\mathbf{n}_{K}\cdot\nabla\widetilde{m}_{\alpha}^{K},\overline{m}_{j_{e(i)}}^{e(i)})_{0,e(i)}\frac{1}{h_{e(i)}^{\frac{1}{2}}}\int_{e(i)}\overline{m}_{j_{e(i)}}^{e(i)}\overline{\varphi}_{i}\\\
&=h_{e(i)}^{\frac{1}{2}}(\mathbf{n}_{K}\cdot\nabla\widetilde{m}_{\alpha}^{K},\overline{m}_{j_{e(i)}}^{e(i)})_{0,e(i)}.\\\
\end{split}$

##### Matrix $\overline{\mathbf{D}}$.

Matrix
$\overline{\mathbf{D}}\in\mathbb{R}^{\widetilde{N}^{K}\times\widetilde{n}^{K}_{p}}$
can be split into two submatrices

$\overline{\mathbf{D}}^{A,B}\in\mathbb{R}^{(p+1)N_{V}\times\widetilde{n}^{K}_{p}},\quad\quad\quad\quad\;\overline{\mathbf{D}}^{C}\in\mathbb{R}^{n_{p-2}^{K}\times\widetilde{n}^{K}_{p}}.$

In particular, we write

$\overline{\mathbf{D}}=\begin{bmatrix}\mbox{ $\overline{\mathbf{D}}^{A,B}$}\\\
\mbox{ $\overline{\mathbf{D}}^{C}$}\\\ \end{bmatrix}.$

Since there is no modification of the bulk degrees of freedom, we have

$\overline{\mathbf{D}}^{C}=\mathbf{D}^{C}.$

Let $j_{e(i)}$ denote the numbering of $\varphi_{i}$ as a basis function on
edge $e(i)$. As for matrix $\overline{\mathbf{D}}^{A,B}$, we apply the
definition of the new enriched edge degrees of freedom and get

$\begin{split}\overline{\mathbf{D}}^{A,B}_{i,\alpha}=\frac{1}{h_{e(i)}^{\frac{1}{2}}}\int_{e(i)}\widetilde{m}_{\alpha}^{K}\overline{m}_{j_{e(i)}}^{e(i)}=\sum_{\ell=1}^{j_{e(i)}}\mathbf{GS}^{e}_{j_{e(i)},\ell}\frac{1}{h_{e(i)}^{\frac{1}{2}}}\int_{e(i)}\widetilde{m}_{\alpha}^{K}m_{\ell}^{e(i)}.\end{split}$

All the integrals appearing on the right-hand side can be computed up to
machine precision as detailed in Appendix A.

##  Appendix C Given a polynomial of degree $p$ in two dimensions, how can we
write it as the Laplacian of a polynomial of degree $p+2$?