|

# On the coupling of the Curved Virtual Element Method with the one-equation
Boundary Element Method for 2D exterior Helmholtz problems

L. Desiderio<EMAIL_ADDRESS>S. Falletta<EMAIL_ADDRESS>M. Ferrari<EMAIL_ADDRESS>L. Scuderi<EMAIL_ADDRESS>Dipartimento di Scienze Matematiche, Fisiche e Informatiche, Università di
Parma,
Parco Area delle Scienze, 53/A, 43124, Parma, Italia  Dipartimento di Scienze
Matematiche “G.L. Lagrange”, Politecnico di Torino,
Corso Duca degli Abruzzi, 24, 10129, Torino, Italia

###### Abstract

We consider the Helmholtz equation defined in unbounded domains, external to
2D bounded ones, endowed with a Dirichlet condition on the boundary and the
Sommerfeld radiation condition at infinity. To solve it, we reduce the
infinite region, in which the solution is defined, to a bounded computational
one, delimited by a curved smooth artificial boundary and we impose on this
latter a non reflecting condition of boundary integral type. Then, we apply
the curved virtual element method in the finite computational domain, combined
with the one-equation boundary element method on the artificial boundary. We
present the theoretical analysis of the proposed approach and we provide an
optimal convergence error estimate in the energy norm. The numerical tests
confirm the theoretical results and show the effectiveness of the new proposed
approach.

###### keywords:

Exterior Helmholtz problems, Curved Virtual Element Method, Boundary Element
Method, Non Reflecting Boundary Condition.

##  1 Introduction

Frequency-domain wave propagation problems defined in unbounded regions,
external to bounded obstacles, turn out to be a difficult physical and
numerical task due to the issue of determining the solution in an infinite
domain. One of the typical techniques to solve such problems is the Boundary
Integral Equation (BIE) method, which allows to reduce by one the dimension of
the problem, requiring only the discretization of the obstacle boundary. Once
the boundary distribution is retrieved by means of a Boundary Element Method
(BEM) [36], the solution of the original problem at each point of the exterior
domain is obtained by computing a boundary integral. However, this procedure
may result not efficient, especially when the solution has to be evaluated at
many points of the infinite domain.

During the last decades much effort has been concentrated on developing
alternative approaches. Among these we mention those based on the coupling of
domain methods, such as Finite Difference Method (FDM), Finite Element Method
(FEM) and the recent Virtual Element Method (VEM), with the BEM. These are
obtained by reducing the unbounded domain to a bounded computational one,
delimited by an artificial boundary, on which a suitable Boundary Integral-Non
Reflecting Boundary Condition (BI-NRBC) is imposed. This latter guarantees
that the artificial boundary is transparent and that no spurious reflections
arise from the resolution of the original problem by means of the interior
domain method applied in the finite computational domain.
The most popular approaches for such a coupling, associated to the use of the
FEM in the interior domain, involve the Green representation of the solution
and are often referred to as the Johnson & Nédélec Coupling (JNC) [29] or the
Costabel & Han Coupling (CHC) [19, 27]. Since the JNC is based on a single
BIE, involving both the single and the double layer integral operators
associated with the fundamental solution, it is known as the _one equation
BEM-FEM coupling_ and it gives rise to a non-symmetric final linear system.

On the contrary, the CHC is based on a couple of BIEs, one of which involves
the second order normal derivative of the fundamental solution (hence a
hypersingular integral operator), and it yields to a symmetric scheme. Despite
the fact that an integration by parts strategy can be applied to weaken the
hypersingularity, the approach turns out to be quite onerous from the
computational point of view, especially in the case of frequency-domain wave
problems for which the accuracy of the BEM is strictly connected to the
frequency parameter and to the density of discretization points per
wavelength. Even if the CHC has been applied in several contexts, among which
we mention the recent paper [25], where the theoretical analysis has been
derived for the solution of the Helmholtz problem by means of a VEM, from the
engineering point of view the JNC is the most natural and appealing way to
deal with unbounded domain problems (for very recent real-life applications
see, for example, [2] and [23]).
In this paper we propose a new approach based on the JNC between the Galerkin
BEM and the Curved Virtual Element Method (CVEM) in the interior of the
computational domain. This choice is based on the fact that the VEM allows to
broaden the classical family of the FEM for the discretization of partial
differential equations for what concerns both the decomposition of domain with
complex geometry and the definition of local high order discrete spaces. In
the standard VEM formulation the discrete spaces, built on meshes made of
polygonal or polyhedral elements, are similar to the usual finite element
spaces with the addition of suitable non-polynomial functions. The novelty of
the VEM consists in defining discrete spaces and degrees of freedom in such a
way that the elementary stiffness and mass matrices can be computed using only
the degrees of freedom, without the need of explicitly knowing the non-
polynomial functions (from which the “virtual” word descends), with a
consequent easiness of implementation even for high approximation orders.
Originally developed as a variational reformulation of the nodal Mimetic
Finite Difference (MFD) method [6, 15, 31], the VEM has been applied to a wide
variety of interior problems (among the most recent papers we refer the reader
to [3, 9, 10, 16]). On the contrary, only few papers deal with VEM applied to
exterior problems, among which the already mentioned [25] and [23]. In this
latter the JNC between the collocation BEM and the VEM has been numerically
investigated for the approximation of the solution of Dirichlet boundary value
problems defined by the 2D Helmholtz equation.

The very satisfactory results we have obtained in [23] have stimulated us to
further investigate on the application of the VEM to the solution of exterior
problems. For this reason, we propose here a novel approach in the CVEM-
Galerkin context that we have studied both from the theoretical and the
numerical point of view. In particular, the choice of the CVEM instead of the
standard (polygonal) VEM relies on the fact that the use of curvilinear
elements allows to avoid the sub-optimal rate of convergence for orders of
accuracy higher than 2, when curvilinear obstacles are considered. Further,
due to the arbitrariness of the choice of the artificial boundary, and dealing
with curved virtual elements, we choose the latter of curvilinear type. It is
worth mentioning that the choice of the CVEM space refers in particular to
that proposed in [8]. For the discretization of the BI-NRBC, we consider a
classical BEM associated to Lagrangian nodal basis functions. As already
remarked, the main challenge in the theoretical analysis is the lack of
ellipticity of the associated bilinear form. However, using the Fredholm
theory for integral operators, it is possible to prove the well-posedness of
the problem in case of computational domains with smooth artificial
boundaries. Moreover, the analysis of the Helmholtz problem and of the
proposed numerical method for its solution, is carried out by interpreting the
new main operators as perturbations of the Laplace ones. We present the
theoretical analysis of the method in a quite general framework and we provide
an optimal error estimate in the energy norm. Since the analysis is based on
the pioneering paper by Jonhson and Nédélec, the smoothness properties of the
artificial boundary represent a key requirement. We remark that, for the
classical Galerkin approach, the breakthrough in the theoretical analysis that
validates the stability of the JNC also in case of non-smooth boundaries, was
proved by Sayas in [37]. However, since we deal with a generalized Galerkin
method, the same analysis can not be straightforwardly applied and needs
further investigations.

The paper is organized as follows: in the next section we present the model
problem for the Helmholtz equation and its reformulation in a bounded region,
by the introduction of the artificial boundary and the associated one equation
BI-NRBC. In Section 3 we introduce the variational formulation of the problem
restricted to the finite computational domain, recalling the corresponding
main theoretical issues, among which existence and uniqueness of the solution.
In Section 4 we apply the Galerkin method providing an error estimate in the
energy norm, for a quite generic class of approximation spaces. Then, in
Section 5 we describe the choice of the CVEM-BEM approximation spaces and we
prove the validity of the error analysis in this specific context.
Additionally, we detail the algebraic formulation of the coupled global
scheme. Finally, in Section 6 we present some numerical results highlighting
the effectiveness of the proposed approach and the validation of the
theoretical results. Furthermore, in the last example we show that the optimal
convergence order of the scheme is guaranteed also when polygonal
computational domains are considered. Finally, some conclusions are drawn in
Section 7.

##  2 The model problem

In a fixed Cartesian coordinates system $\mathbf{x}=\left(x_{1},x_{2}\right)$,
we consider an open bounded domain $\Omega_{\tiny{0}}\subset\mathbf{R}^{2}$
with a Lipschitz boundary $\Gamma_{\tiny{0}}$ having positive Lebesgue
measure. We denote by
$\Omega_{e}:=\mathbf{R}^{2}\setminus\overline{\Omega}_{\tiny{0}}$ the exterior
unbounded domain (see Figure 1 (a)) and we consider the following frequency-
domain wave propagation problem:

$\displaystyle\Delta
u_{e}(\mathbf{x})+\kappa^{2}u_{e}(\mathbf{x})=-f(\mathbf{x})$
$\displaystyle\mathbf{x}\in\Omega_{e},$ (2.1a) $\displaystyle
u_{e}(\mathbf{x})=g(\mathbf{x})$ $\displaystyle\mathbf{x}\in\Gamma_{0},$
(2.1b)
$\displaystyle\lim\limits_{\|\mathbf{x}\|\rightarrow\infty}\|\mathbf{x}\|^{\frac{1}{2}}\left(\nabla
u_{e}(\mathbf{x})\cdot\frac{\mathbf{x}}{\|\mathbf{x}\|}-\imath\kappa
u_{e}(\mathbf{x})\right)=0.$ (2.1c)

In the above problem, Equation (2.1a) is known as the Helmholtz equation, with
source term $f\in L^{2}(\Omega_{e})$, Equation (2.1b) represents a boundary
condition of Dirichlet type with datum $g$, and Equation (2.1c) is the
Sommerfeld radiation condition, that ensures the appropriate behaviour of the
complex-valued unknown function $u_{e}$ at infinity. Furthermore, $\nabla$ and
$\Delta$ denote the nabla and Laplace operators, respectively, and $\imath$
stands for the imaginary unit.
We recall that the wave number $\kappa$ is often real and constant, and it is
complex if the propagation medium is energy absorbing, or a function of the
space if the medium is inhomogeneous. Here, we suppose that $\kappa$ is real,
positive and constant.
In the sequel we assume that $g\in H^{\nicefrac{{1}}{{2}}}(\Gamma_{0})$, to
guarantee existence and uniqueness of the solution $u_{e}$ of Problem (2.1) in
the Sobolev space $H^{1}_{\text{loc}}(\Omega_{e})$ (see [18]).

As many practical situations require, we aim at determining the solution
$u_{e}$ of Problem (2.1) in a bounded subregion of $\Omega_{e}$ surrounding
$\Omega_{0}$. To this end, we introduce an artificial boundary $\Gamma$ which
allows decomposing $\Omega_{e}$ into a finite computational domain $\Omega$,
bounded internally by $\Gamma_{0}$ and externally by $\Gamma$, and an infinite
residual one, denoted by $\Omega_{\infty}$, as depicted in Figure 1 (b). We
choose $\Gamma$ such that $\text{supp}(f)$ is a bounded subset of $\Omega$. We
assume that $\Gamma_{0}$ and $\Gamma$ are made up of a finite number of curves
of class $C^{m+1}$, with $m\geq 0$, so that $\Omega$ is a domain with piece-
wise smooth boundaries. Moreover, we assume that $\Gamma$ is a Lyapunov
regular contour, i.e. the gradient of any local parametrization is Hölder
continuous.

Figure 1: Model problem setting.

Denoting by $u$ and $u_{\infty}$ the restrictions of the solution $u_{e}$ to
$\Omega$ and $\Omega_{\infty}$ respectively, and by $\mathbf{n}$ and
$\mathbf{n}_{\infty}$ the unit normal vectors on $\Gamma$ pointing outside
$\Omega$ and $\Omega_{\infty}$, we impose the following compatibility and
equilibrium conditions on $\Gamma$ (recall that
$\mathbf{n}_{\infty}=-\mathbf{n}$):

$u(\mathbf{x})=u_{\infty}(\mathbf{x}),\qquad\frac{\partial
u}{\partial{\mathbf{n}}}(\mathbf{x})=-\frac{\partial
u_{\infty}}{\partial{\mathbf{n}_{\infty}}}(\mathbf{x}),\qquad\mathbf{x}\in\Gamma.$
(2.2)

In the above relations and in the sequel we omit, for simplicity, the use of
the trace operator to indicate the restriction of $H^{1}$ functions to the
boundary $\Gamma$ from the exterior or interior. In order to obtain a well
posed problem in $\Omega$, we need to impose a proper boundary condition on
$\Gamma$. It is known that the solution $u_{\infty}$ in $\Omega_{\infty}$ can
be represented by the following Kirchhoff formula:

$u_{\infty}(\mathbf{x})=\int_{\Gamma}G_{\kappa}(\mathbf{x},\mathbf{y})\frac{\partial
u_{\infty}}{\partial\mathbf{n}_{\infty}}(\mathbf{y})\,\differential\Gamma_{\mathbf{y}}-\int_{\Gamma}\frac{\partial
G_{\kappa}}{\partial\mathbf{n}_{\infty,\mathbf{y}}}(\mathbf{x},\mathbf{y})u_{\infty}(\mathbf{y})\,\differential\Gamma_{\mathbf{y}}\qquad\mathbf{x}\in\Omega_{\infty}\setminus\Gamma,$
(2.3)

in which $G_{\kappa}$ is the fundamental solution of the 2D Helmholtz problem
and $\mathbf{n}_{\infty,\mathbf{y}}$ denotes the normal unit vector with
initial point in $\mathbf{y}\in\Gamma$. The expression of $G_{\kappa}$ and of
its normal derivative in (2.3) are given by

$G_{\kappa}(\mathbf{x},\mathbf{y}):=\frac{\imath}{4}H_{0}^{(1)}(\kappa
r)\quad\text{and}\quad\frac{\partial
G_{\kappa}}{\partial\mathbf{n}_{\infty,\mathbf{y}}}(\mathbf{x},\mathbf{y})=\frac{\imath\kappa}{4}\frac{\mathbf{r}\cdot\mathbf{n}_{\infty,\mathbf{y}}}{r}H_{1}^{(1)}(\kappa
r),$

where $r=|\mathbf{r}|=|\mathbf{x}-\mathbf{y}|$ represents the distance between
the source point $\mathbf{x}$ and the field point $\mathbf{y}$, and
$H_{m}^{(1)}$ denotes the $m$-th order Hankel function of the first kind.

We introduce the single-layer integral operator $\text{V}_{\kappa}\colon
H^{-\nicefrac{{1}}{{2}}}(\Gamma)\to H^{\nicefrac{{1}}{{2}}}(\Gamma)$

$\text{V}_{\kappa}\psi(\mathbf{x}):=\int_{\Gamma}G_{\kappa}(\mathbf{x},\mathbf{y})\psi(\mathbf{y})\,\differential\Gamma_{\mathbf{y}},\qquad\mathbf{x}\in\Gamma$

and the double-layer integral operator $\text{K}_{\kappa}\colon
H^{\nicefrac{{1}}{{2}}}(\Gamma)\to H^{\nicefrac{{1}}{{2}}}(\Gamma)$

$\text{K}_{\kappa}\varphi(\mathbf{x}):=-\int_{\Gamma}\frac{\partial
G_{\kappa}}{\partial\mathbf{n}_{\infty,\mathbf{y}}}(\mathbf{x},\mathbf{y})\varphi(\mathbf{y})\,\differential\Gamma_{\mathbf{y}},\qquad\mathbf{x}\in\Gamma,$

which are continuous for all $\kappa>0$ (see [28]). Then, the trace of (2.3)
on $\Gamma$ reads (see [18])

$\frac{1}{2}u_{\infty}(\mathbf{x})-\text{V}_{\kappa}\frac{\partial
u_{\infty}}{\partial\mathbf{n}_{\infty}}({\mathbf{x}})-\text{K}_{\kappa}u_{\infty}(\mathbf{x})=0,\qquad\mathbf{x}\in\Gamma.$
(2.4)

Equation (2.4), which expresses the natural relation that $u_{\infty}$ and its
normal derivative have to satisfy at each point of the artificial boundary, is
imposed on $\Gamma$ as an exact (non local) BI-NRBC to solve Problem (2.1) in
the finite computational domain. Thus, taking into account the compatibility
and equilibrium conditions (2.2), and introducing the notation
$\displaystyle\lambda:=\frac{\partial u}{\partial\mathbf{n}}$, the new problem
defined in the domain of interest $\Omega$ takes the form:

$\displaystyle\Delta u(\mathbf{x})+\kappa^{2}u(\mathbf{x})=-f(\mathbf{x})$
$\displaystyle\mathbf{x}\in\Omega$ (2.5a) $\displaystyle
u(\mathbf{x})=g(\mathbf{x})$ $\displaystyle\mathbf{x}\in\Gamma_{0}$ (2.5b)
$\displaystyle\frac{1}{2}u(\mathbf{x})+\text{V}_{\kappa}\lambda({\mathbf{x}})-\text{K}_{\kappa}u(\mathbf{x})=0$
$\displaystyle\mathbf{x}\in\Gamma.$ (2.5c)

We point out that $\lambda$, which is defined on the boundary $\Gamma$ in
general by means of a trace operator (see [34]), is an additional unknown
function.

For the theoretical analysis we will present in the forthcoming section, we
further need to introduce the fundamental solution $G_{0}$ of the Laplace
equation and its normal derivative:

$G_{0}(\mathbf{x},\mathbf{y}):=-\frac{1}{2\pi}\log
r\quad\text{and}\quad\frac{\partial
G_{0}}{\partial\mathbf{n}_{\infty,\mathbf{y}}}(\mathbf{x},\mathbf{y})=\frac{1}{2\pi}\frac{\mathbf{r}\cdot\mathbf{n}_{\infty,\mathbf{y}}}{r^{2}}.$

Denoting by $\text{V}_{0}$ and $\text{K}_{0}$ the associated single and double
layer operators, the following regularity property of the operators
$\text{V}_{\kappa}-\text{V}_{0}$ and $\text{K}_{\kappa}-\text{K}_{0}$ holds.

######  Lemma 2.1.

The operators
$\text{V}_{\kappa}-\text{V}_{0}:H^{-\nicefrac{{1}}{{2}}}(\Gamma)\rightarrow
H^{\nicefrac{{5}}{{2}}}(\Gamma)$ and
$\text{K}_{\kappa}-\text{K}_{0}:H^{\nicefrac{{1}}{{2}}}(\Gamma)\rightarrow
H^{\nicefrac{{3}}{{2}}}(\Gamma)$ are continuous.

###### Proof.

We preliminary recall that the Hankel functions $H_{m}^{(1)}$, with $m=0,1$,
have the following asymptotic behaviour when $r\rightarrow 0$ (see formulae
(2.14) and (2.15) in [35]):

$\displaystyle H_{0}^{(1)}(r)=\frac{\imath 2}{\pi}\log{(r)}+1+\frac{\imath
2}{\pi}\left(\gamma-\log{(2)}\right)+O(r^{2}),$ (2.6a) $\displaystyle
H_{1}^{(1)}(r)=-\frac{\imath 2}{\pi r}+O(1)$ (2.6b)

where $\gamma\simeq 0.577216$ is the Euler constant. Then it easily follows
that, when $r\rightarrow 0$

$\displaystyle
G_{\kappa}(\mathbf{x},\mathbf{y})-G_{0}(\mathbf{x},\mathbf{y})=\frac{\imath}{4}-\frac{1}{2\pi}\left(\gamma-\log{\left(\frac{\kappa}{2}\right)}\right)+O(r^{2}),$
$\displaystyle\frac{\partial
G_{\kappa}}{\partial\mathbf{n}_{\infty,\mathbf{y}}}(\mathbf{x},\mathbf{y})-\frac{\partial
G_{0}}{\partial\mathbf{n}_{\infty,\mathbf{y}}}(\mathbf{x},\mathbf{y})=O(1).$

Following [28] (see Section 7.1), we can therefore deduce that
$G_{\kappa}-G_{0}$ and $\displaystyle\frac{\partial
G_{\kappa}}{\partial\mathbf{n}_{\infty,\mathbf{y}}}-\frac{\partial
G_{0}}{\partial\mathbf{n}_{\infty,\mathbf{y}}}$ are kernel functions with
pseudo-homogeneous expansions of degree 2 and 0, respectively. From these
properties, and proceeding as in [36] (see Remark 3.13), the thesis is proved.
∎

######  Remark 1.

Similarly, since $G_{0}$ is a kernel function with a pseudo-homogeneous
expansion of degree 0, we can deduce that the operator
$\text{V}_{0}:H^{s}(\Gamma)\rightarrow H^{s+1}(\Gamma)$ is continuous for all
$s\in\mathbf{R}$.

##  3 The weak formulation of the model problem

We start by noting that, as usual, we can reduce the non homogeneous boundary
condition on $\Gamma_{0}$ in (2.1) to a homogeneous one by splitting $u_{e}$
as the sum of a suitable fixed function in
$H^{1}_{g,\Gamma_{0}}(\Omega):=\\{u\in H^{1}(\Omega):u=g\,\text{ on
}\Gamma_{0}\\}$ satisfying the Sommerfeld radiation condition and of an
unknown function belonging to the space $H^{1}_{0,\Gamma_{0}}(\Omega)$.
Therefore, from now on, we consider Problem (2.5) with $g=0$.

In order to derive the weak form of Problem (2.5), we introduce the bilinear
forms $a:H^{1}(\Omega)\times H^{1}(\Omega)\rightarrow\mathbf{C}$ and
$m:L^{2}(\Omega)\times L^{2}(\Omega)\rightarrow\mathbf{C}$ given by

$a(u,v):=\int\limits_{\Omega}\nabla u(\mathbf{x})\cdot\nabla
v(\mathbf{x})\,\differential\mathbf{x}\qquad\text{and}\qquad
m(u,v):=\int\limits_{\Omega}u(\mathbf{x})v(\mathbf{x})\,\differential\mathbf{x},$
(3.1)

and the $L^{2}(\Gamma)$-inner product
$(\cdot,\cdot)_{\Gamma}:L^{2}(\Gamma)\times
L^{2}(\Gamma)\rightarrow\mathbf{C}$

$(\lambda,v)_{\Gamma}=\int\limits_{\Gamma}\lambda(\mathbf{x})v(\mathbf{x})\differential\Gamma_{\mathbf{x}},$

extended to the duality pairing $\langle\cdot,\cdot\rangle_{\Gamma}$ on
$H^{-\nicefrac{{1}}{{2}}}(\Gamma)\times H^{\nicefrac{{1}}{{2}}}(\Gamma)$.

The variational formulation of Problem (2.5) consists in finding $u\in
H_{0,\Gamma_{0}}^{1}(\Omega)$ and $\lambda\in
H^{-\nicefrac{{1}}{{2}}}(\Gamma)$ such that

$\displaystyle
a(u,v)-\kappa^{2}m(u,v)-\langle\lambda,v\rangle_{\Gamma}=m(f,v)$
$\displaystyle\forall\,v\in H^{1}_{0,\Gamma_{0}}(\Omega),$ (3.2a)
$\displaystyle\langle\mu,\left(\frac{1}{2}I-\text{K}_{\kappa}\right)u\rangle_{\Gamma}+\langle\mu,\text{V}_{\kappa}\lambda\rangle_{\Gamma}=0$
$\displaystyle\forall\,\mu\in H^{-\nicefrac{{1}}{{2}}}(\Gamma),$ (3.2b)

where $I$ stands for the identity operator. In order to reformulate the above
problem as an equation in operator form, following [29], we consider the
Hilbert space $V:=H^{1}_{0,\Gamma_{0}}(\Omega)\times
H^{-\nicefrac{{1}}{{2}}}(\Gamma)$, equipped with the norm

$\left\|\hat{u}\right\|_{V}^{2}:=\left\|u\right\|_{H^{1}(\Omega)}^{2}+\left\|\lambda\right\|_{H^{-\nicefrac{{1}}{{2}}}(\Gamma)}^{2},\quad\text{for}\
\hat{u}=(u,\lambda),$

induced by the scalar product

$\left(\hat{u},\hat{v}\right)_{V}:=\left(u,v\right)_{H^{1}(\Omega)}+\left(\lambda,\mu\right)_{H^{-\nicefrac{{1}}{{2}}}(\Gamma)},\quad\text{for}\
\hat{u}=(u,\lambda),\,\hat{v}=(v,\mu).$

We introduce the bilinear form $\mathcal{A}_{\kappa}:V\times
V\rightarrow\mathbf{C}$ defined, for $\hat{u}=(u,\lambda)$ and
$\hat{v}=(v,\mu)$, by

$\mathcal{A}_{\kappa}(\hat{u},\hat{v}):=a(u,v)-\kappa^{2}m(u,v)-\langle\lambda,v\rangle_{\Gamma}+\langle\mu,u\rangle_{\Gamma}+2\langle\mu,\text{V}_{\kappa}\lambda\rangle_{\Gamma}-2\langle\mu,\text{K}_{\kappa}u\rangle_{\Gamma},$
(3.3)

and the linear continuous operator $\mathcal{L}_{f}:V\rightarrow\mathbf{C}$

$\mathcal{L}_{f}(\hat{v}):=m(f,v),\quad\hat{v}=(v,\mu).$

Thus, Problem (3.2) can be rewritten as follows: find $\hat{u}\in V$ such that

$\mathcal{A}_{\kappa}(\hat{u},\hat{v})=\mathcal{L}_{f}(\hat{v})\qquad\forall\,\hat{v}\in
V.$ (3.4)

The well-posedness of the above problem has been proved in [32] (see Theorem
3.2), provided that $\kappa^{2}$ is not an eigenvalue of the Dirichlet-Laplace
problem in $\Omega$.

For what follows, it will be useful to rewrite $\mathcal{A}_{\kappa}$ by means
of the bilinear forms $\mathcal{B}_{\kappa},\mathcal{K}_{\kappa}:V\times
V\rightarrow\mathbf{C}$, defined as:

$\displaystyle\mathcal{A}_{\kappa}(\hat{u},\hat{v}):=\mathcal{B}_{\kappa}(\hat{u},\hat{v})+\mathcal{K}_{\kappa}(\hat{u},\hat{v})$
(3.5a)
$\displaystyle\mathcal{B}_{\kappa}(\hat{u},\hat{v}):=a(u,v)-\kappa^{2}m(u,v)-\langle\lambda,v\rangle_{\Gamma}+\langle\mu,u\rangle_{\Gamma}+2\langle\mu,\text{V}_{\kappa}\lambda\rangle_{\Gamma}$
(3.5b)
$\displaystyle\mathcal{K}_{\kappa}(\hat{u},\hat{v}):=-2\langle\mu,\text{K}_{\kappa}u\rangle_{\Gamma}$
(3.5c)

for $\hat{u}=(u,\lambda),\hat{v}=(v,\mu)\in V$. Due to the continuity property
of both the operators $\text{V}_{\kappa}$ and $\text{K}_{\kappa}$ when
$\kappa$ is real and non-negative (see [25]), by using the trace theorem and
the Cauchy-Schwarz inequality, it is easy to prove that the corresponding
linear mappings
$\mathcal{A}_{\kappa},\mathcal{B}_{\kappa},\mathcal{K}_{\kappa}:V\rightarrow
V^{{}^{\prime}}$, defined by

$\left(\mathcal{A}_{\kappa}\hat{u}\right)(\hat{v}):=\mathcal{A}_{\kappa}(\hat{u},\hat{v}),\qquad\left(\mathcal{B}_{\kappa}\hat{u}\right)(\hat{v}):=\mathcal{B}_{\kappa}(\hat{u},\hat{v}),\qquad\left(\mathcal{K}_{\kappa}\hat{u}\right)(\hat{v}):=\mathcal{K}_{\kappa}(\hat{u},\hat{v}),$

are continuous from $V$ to its dual $V^{{}^{\prime}}$. Finally, we introduce
the adjoint operators
$\mathcal{A}^{*}_{\kappa},\mathcal{B}^{*}_{\kappa}:V\rightarrow
V^{{}^{\prime}}$ defined by:

$\displaystyle\left(\mathcal{A}^{*}_{\kappa}\hat{v}\right)(\hat{u}):=\left(\mathcal{A}_{\kappa}\hat{u}\right)(\hat{v})=\mathcal{A}_{\kappa}(\hat{u},\hat{v})$
$\displaystyle\left(\mathcal{B}^{*}_{\kappa}\hat{v}\right)(\hat{u}):=\left(\mathcal{B}_{\kappa}\hat{u}\right)(\hat{v})=\mathcal{B}_{\kappa}(\hat{u},\hat{v}).$

In the following remarks we recall classical results about the afore
introduced maps.

######  Remark 2.

Theorem 3.2 in [32] and the closed graph theorem ensure that, if $\kappa^{2}$
is not an eigenvalue of the Dirichlet-Laplace problem in $\Omega$, the inverse
linear mappings
$\mathcal{A}_{\kappa}^{-1},\mathcal{A}_{\kappa}^{*-1}:V^{{}^{\prime}}\rightarrow
V$ are continuous.

######  Remark 3.

Denoting by $H_{0}^{-\nicefrac{{1}}{{2}}}(\Gamma):=\left\\{\lambda\in
H^{-\nicefrac{{1}}{{2}}}(\Gamma)\ :\
\langle\lambda,1\rangle_{\Gamma}=0\right\\}$, we set
$\widetilde{V}:=H^{1}_{0,\Gamma_{0}}(\Omega)\times
H^{-\nicefrac{{1}}{{2}}}_{0}(\Gamma)$. It has been proved in [29] (see Lemmas
1, 2 and 3) that the mappings
$\mathcal{A}_{0},\mathcal{A}^{*}_{0},\mathcal{B}_{0},\mathcal{B}^{*}_{0}:\widetilde{V}\rightarrow\widetilde{V}^{{}^{\prime}}$
are isomorphisms. Moreover, for $s\geq 0$, the mappings
$\mathcal{A}_{0}^{-1},\mathcal{A}_{0}^{*-1},\mathcal{B}_{0}^{-1},\mathcal{B}_{0}^{*-1}:H^{s-1}(\Omega)\times
H^{s-\nicefrac{{1}}{{2}}}(\Gamma)\times H^{s+\nicefrac{{1}}{{2}}}(\Gamma)\to
H^{s+1}(\Omega)\times H^{s-\nicefrac{{1}}{{2}}}(\Gamma)$ are continuous.
Finally, we recall that $\mathcal{B}_{0}$ is coercive in the
$\widetilde{V}$-norm.

##  4 The Galerkin method

In what follows, the notation $Q_{1}\apprle Q_{2}$ (resp. $Q_{1}\apprge
Q_{2}$) means that the quantity $Q_{1}$ is bounded from above (resp. from
below) by $c\,Q_{2}$, where $c$ is a positive constant that, unless explicitly
stated, does not depend on any relevant parameter involved in the definition
of $Q_{1}$ and $Q_{2}$.

In order to describe the Galerkin approach applied to (3.4), we introduce a
sequence of unstructured meshes $\\{\mathcal{T}_{h}\\}_{h>0}$, that represent
coverages of the domain $\Omega$ with a finite number of elements $E$, having
diameter $h_{E}$. The mesh width $h>0$, related to the spacing of the grid, is
defined as $h:=\underset{E\in\mathcal{T}_{h}}{\max}h_{E}$. Moreover, we denote
by $\mathcal{T}_{h}^{\Gamma}$ the decomposition of the artificial boundary
$\Gamma$, inherited from $\mathcal{T}_{h}$, which consists of curvilinear
parts joined with continuity. We suppose that for each $h$ and for each
element $E\in\mathcal{T}_{h}$ there exists a constant $\varrho>0$ such that
the following assumptions are fulfilled:

  1. ($A.{1}$)

$E$ is star-shaped with respect to a ball of radius greater than $\varrho
h_{E}$;

  2. ($A.{2}$)

the length of any (eventually curved) edge of $E$ is greater than $\varrho
h_{E}$.

We introduce the splitting of the bilinear forms $a$ and $m$ defined in (3.1)
into a sum of local bilinear forms
$a^{\text{\tiny{E}}},m^{\text{\tiny{E}}}:H^{1}(E)\times
H^{1}(E)\rightarrow\mathbf{C}$, associated to the elements $E$ of the
decomposition of $\Omega$:

$\displaystyle a(u,v)$
$\displaystyle=\sum\limits_{E\in\mathcal{T}_{h}}a^{\text{\tiny{E}}}(u,v):=\sum\limits_{E\in\mathcal{T}_{h}}\int\limits_{E}\nabla
u(\mathbf{x})\cdot\nabla v(\mathbf{x})\,\differential\mathbf{x},$
$\displaystyle m(u,v)$
$\displaystyle=\sum\limits_{E\in\mathcal{T}_{h}}m^{\text{\tiny{E}}}(u,v):=\sum\limits_{E\in\mathcal{T}_{h}}\int\limits_{E}u(\mathbf{x})v(\mathbf{x})\,\differential\mathbf{x}.$

Then, for any $k\in\mathbf{N}$, denoting by $P_{k}(E)$ the space of
polynomials of degree $k$ defined on $E$, we introduce the local polynomial
$H^{1}$-projection $\Pi_{k}^{\nabla,E}:H^{1}(E)\rightarrow P_{k}(E)$, defined
such that for every $v\in H^{1}(E)$:

$\begin{cases}\displaystyle\int_{E}\nabla\Pi_{k}^{\nabla,E}v\cdot\nabla
q\,\differential E=\displaystyle\int_{E}\nabla v\cdot\nabla q\,\differential
E\qquad\forall\,q\in P_{k}(E),\\\\[10.0pt] \displaystyle\int_{\partial
E}\Pi_{k}^{\nabla,E}v\,\differential s=\displaystyle\int_{\partial
E}v\,\differential s\end{cases}$ (4.1)

and the local polynomial $L^{2}$-projection operator
$\Pi_{k}^{0,E}:L^{2}(E)\rightarrow P_{k}(E)$, defined for all $v\in L^{2}(E)$
such that

$\int_{E}\Pi_{k}^{0,E}v\,q\,\differential E=\int_{E}v\,q\,\differential
E\qquad\forall\,q\in P_{k}(E).$ (4.2)

From the definition of $\Pi_{k}^{\nabla,E}$ and of $\Pi_{k}^{0,E}$, it
follows:

$\displaystyle
a^{\text{\tiny{E}}}\left(\Pi_{k}^{\nabla,E}v,q\right)=a^{\text{\tiny{E}}}\left(v,q\right)$
$\displaystyle\forall\,q\in P_{k}(E),$ $\displaystyle
m^{\text{\tiny{E}}}\left(\Pi_{k}^{0,E}v,q\right)=m^{\text{\tiny{E}}}\left(v,q\right)$
$\displaystyle\forall\,q\in P_{k}(E).$

Moreover, since $\Omega$ is the union of star-shaped domains
$E\in\mathcal{T}_{h}$, the local polynomial projectors $\Pi_{k}^{\nabla,E}$
and $\Pi_{k}^{0,E}$ can be extended to the global projectors
$\Pi_{k}^{\nabla}:H^{1}(\Omega)\rightarrow P_{k}(\mathcal{T}_{h})$ and
$\Pi_{k}^{0}:L^{2}(\Omega)\rightarrow P_{k}(\mathcal{T}_{h})$ as follows:

$\displaystyle\left(\Pi_{k}^{\nabla}v\right)_{|_{E}}:=\Pi_{k}^{\nabla,E}v_{|_{E}}$
$\displaystyle\forall\,v\in H^{1}(\Omega)$
$\displaystyle\left(\Pi_{k}^{0}v\right)_{|_{E}}:=\Pi_{k}^{0,E}v_{|_{E}}$
$\displaystyle\forall\,v\in L^{2}(\Omega),$

$P_{k}(\mathcal{T}_{h})$ being the space of piecewise polynomials with respect
to the decomposition $\mathcal{T}_{h}$ of $\Omega$.

In the following lemma we prove a polynomial approximation property of the
above defined projectors. To this aim, since we shall deal with functions
belonging to the space
$H^{1}(\mathcal{T}_{h}):=\underset{E\in\mathcal{T}_{h}}{\prod}H^{1}(E)$, we
need to introduce the following broken $H^{1}$ norm

$\|v\|_{H^{1}(\mathcal{T}_{h})}:=\left(\sum_{E\in\mathcal{T}_{h}}\|v\|^{2}_{H^{1}(E)}\right)^{1/2}.$

######  Lemma 4.1.

Assuming ($A.{1}$), for all $v\in H^{s+1}(\Omega)$ with $0\leq s\leq k$, it
holds:

$\left\|v-\Pi_{k}^{0}v\right\|_{L^{2}(\Omega)}\apprle
h^{s+1}\left\|v\right\|_{H^{s+1}(\Omega)}.$ (4.4)

Moreover, for all $v\in H^{s+1}(\Omega)$ with $1\leq s\leq k$, it holds:

$\left\|v-\Pi_{k}^{\nabla}v\right\|_{H^{1}(\mathcal{T}_{h})}\apprle
h^{s}\left\|v\right\|_{H^{s+1}(\Omega)}.$ (4.5)

###### Proof.

Let us denote by $R_{E}>0$ the radius of the ball in $E\in\mathcal{T}_{h}$
satisfying ($A.{1}$). For any $v\in L^{2}(E)$ and $q\in P_{k}(E)$ we can write

$\left\|v-\Pi_{k}^{0,E}v\right\|_{L^{2}(E)}^{2}=\left(v-\Pi_{k}^{0,E}v,v-q\right)_{L^{2}(E)}+\left(v-\Pi_{k}^{0,E}v,q-\Pi_{k}^{0,E}v\right)_{L^{2}(E)}=\left(v-\Pi_{k}^{0,E}v,v-q\right)_{L^{2}(E)},$

where we have used (4.2) together with $q-\Pi_{k}^{0,E}v\in P_{k}(E)$. Then,
by applying the Cauchy-Schwarz inequality, we easily get

$\left\|v-\Pi_{k}^{0,E}v\right\|_{L^{2}(E)}\leq\left\|v-q\right\|_{L^{2}(E)}\qquad\forall\,q\in
P_{k}(E).$

Now, assuming $v\in H^{s+1}(E)$ with $0\leq s\leq k$ and using the Bramble-
Hilbert Lemma (see Lemma 4.3.8 in [14]), we have

$\left\|v-\Pi_{k}^{0,E}v\right\|_{L^{2}(E)}\leq\underset{q\in
P_{k}(E)}{\text{inf}}\left\|v-q\right\|_{L^{2}(E)}\apprle
C\hskip-2.84544pt\left(\frac{h_{E}}{R_{E}}\right)h_{E}^{s+1}\left\|v\right\|_{H^{s+1}(E)},$

where the implicit constant depends only on $k$ and
$C:\mathbf{R}^{+}\rightarrow\mathbf{R}^{+}$ is an increasing function. Since,
by virtue of Assumption ($A.{1}$), the function
$C\hskip-2.27626pt\left(\nicefrac{{h_{E}}}{{R_{E}}}\right)$ is uniformly
bounded, we can easily get (4.4). Finally, inequality (4.5) can be proved
similarly. ∎

###  4.1 The discrete variational formulation

We present here a class of Galerkin type discretizations of Problem (3.4),
which includes, but is not limited to, VEMs. In Section 5 we will give an
example of CVEM that falls in the framework considered.

Let $Q_{h}^{k}\subset H_{0,\Gamma_{0}}^{1}(\Omega)$ and $X_{h}^{k}\subset
H^{-\nicefrac{{1}}{{2}}}(\Gamma)$ denote two finite dimensional spaces
associated to the meshes $\mathcal{T}_{h}$ and $\mathcal{T}_{h}^{\Gamma}$,
respectively. Introducing the discrete space $V_{h}^{k}:=Q_{h}^{k}\times
X_{h}^{k}$, the Galerkin method applied to Problem (3.2) reads: find
$\hat{u}_{h}\in V_{h}^{k}$ such that

$\mathcal{A}_{\kappa,h}(\hat{u}_{h},\hat{v}_{h}):=\mathcal{B}_{\kappa,h}(\hat{u}_{h},\hat{v}_{h})+\mathcal{K}_{\kappa}(\hat{u}_{h},\hat{v}_{h})=\mathcal{L}_{f,h}(\hat{v}_{h})\quad\forall\,\hat{v}_{h}\in
V_{h}^{k},$ (4.6)

where $\mathcal{A}_{\kappa,h},\mathcal{B}_{\kappa,h}:V_{h}^{k}\times
V_{h}^{k}\rightarrow\mathbf{C}$ and
$\mathcal{L}_{f,h}:V_{h}^{k}\rightarrow\mathbf{C}$ are suitable approximations
of $\mathcal{A}_{\kappa}$, $\mathcal{B}_{\kappa}$ and $\mathcal{L}_{f}$,
respectively.

In order to prove existence and uniqueness of the solution $\hat{u}_{h}\in
V_{h}^{k}$ and to derive a priori error estimates, we preliminary introduce
some assumptions on the discrete spaces, on the bilinear form
$\mathcal{B}_{\kappa,h}$ and on the linear operator $\mathcal{L}_{f,h}$.

We assume that the following properties for $Q_{h}^{k}$, $X_{h}^{k}$ and
$\widetilde{X}_{h}^{k}:=X_{h}^{k}\cap H_{0}^{-\nicefrac{{1}}{{2}}}(\Gamma)$
hold: for $1\leq s\leq k$

  1. ($H1.{a}$)

$\underset{v_{h}\in
Q_{h}^{k}}{\text{inf}}\left\|v-v_{h}\right\|_{H^{1}(\Omega)}\apprle
h^{s}\left\|v\right\|_{H^{s+1}(\Omega)}\hskip 5.69046pt\quad\quad\forall\,v\in
H^{s+1}(\Omega)$;

  2. ($H1.{b}$)

$\underset{\mu_{h}\in
X_{h}^{k}}{\text{inf}}\left\|\mu-\mu_{h}\right\|_{H^{-\nicefrac{{1}}{{2}}}(\Gamma)}\apprle
h^{s}\left\|\mu\right\|_{H^{s-\nicefrac{{1}}{{2}}}(\Gamma)}\hskip
5.69046pt\quad\forall\,\mu\in H^{s-\nicefrac{{1}}{{2}}}(\Gamma)$;

  3. ($H1.{c}$)

$\underset{\mu_{\scriptscriptstyle
0h}\in\tilde{X}_{h}^{k}}{\text{inf}}\left\|\mu_{\scriptscriptstyle
0}-\mu_{\scriptscriptstyle
0h}\right\|_{H^{-\nicefrac{{1}}{{2}}}(\Gamma)}\apprle
h^{s}\left\|\mu_{\scriptscriptstyle
0}\right\|_{H^{s-\nicefrac{{1}}{{2}}}(\Gamma)}\hskip
5.69046pt\quad\forall\,\mu_{\scriptscriptstyle 0}\in
H^{s-\nicefrac{{1}}{{2}}}(\Gamma)\cap H_{0}^{-\nicefrac{{1}}{{2}}}(\Gamma)$.

According to the definition of the $\|\cdot\|_{V}$ norm, the above assumptions
ensure the following approximation property for the product spaces $V_{h}^{k}$
and $\tilde{V}_{h}^{k}:=Q_{h}^{k}\times\tilde{X}_{h}^{k}$:

  * •

given $\hat{v}=(v,\mu)\in H^{s+1}(\Omega)\times
H^{s-\nicefrac{{1}}{{2}}}(\Gamma)$, there exists
$\hat{v}_{h}^{I}=(v_{h}^{I},\mu_{h}^{I})\in V_{h}^{k}$ such that

$\left\|\hat{v}-\hat{v}_{h}^{I}\right\|_{V}\apprle
h^{s}\left(\left\|v\right\|_{H^{s+1}(\Omega)}+\left\|\mu\right\|_{H^{s-\nicefrac{{1}}{{2}}}(\Gamma)}\right);$
(4.7)

  * •

given $\hat{v}_{\scriptscriptstyle 0}=(v,\mu_{\scriptscriptstyle 0})\in
H^{s+1}(\Omega)\times(H^{s-\nicefrac{{1}}{{2}}}(\Gamma)\cap
H_{0}^{-\nicefrac{{1}}{{2}}}(\Gamma))$, there exists
$\hat{v}_{\scriptscriptstyle 0h}^{I}=(v_{h}^{I},\mu_{\scriptscriptstyle
0h}^{I})\in\tilde{V}_{h}^{k}$ such that

$\left\|\hat{v}_{\scriptscriptstyle 0}-\hat{v}_{\scriptscriptstyle
0h}^{I}\right\|_{V}\apprle
h^{s}\left(\left\|v\right\|_{H^{s+1}(\Omega)}+\left\|\mu_{\scriptscriptstyle
0}\right\|_{H^{s-\nicefrac{{1}}{{2}}}(\Gamma)}\right).$ (4.8)

Concerning the bilinear form $\mathcal{B}_{\kappa,h}$, we assume that for all
$\kappa\geq 0$:

  1. ($H2.{a}$)

$k$-consistency:
$\mathcal{B}_{\kappa,h}(\hat{q},\hat{v}_{h})=\mathcal{B}_{\kappa}(\hat{q},\hat{v}_{h})\quad\forall\,\hat{q}\in
P_{k}(\mathcal{T}_{h})\times X_{h}^{k}\ \text{and}\ \forall\,\hat{v}_{h}\in
V_{h}^{k}$;

  2. ($H2.{b}$)

continuity in $V$-norm:
$\left|\mathcal{B}_{\kappa,h}(\hat{v}_{h},\hat{w}_{h})\right|\apprle\left\|\hat{v}_{h}\right\|_{V}\left\|\hat{w}_{h}\right\|_{V}\quad\forall\hat{v}_{h},\hat{w}_{h}\in
V_{h}^{k}$.

######  Remark 4.

It is worth noting that, in Assumption ($H2.{a}$), the evaluation of the
bilinear form $\mathcal{B}_{\kappa}$ is well defined provided that the
computation of the bilinear form $a(\cdot,\cdot)$ is split into the sum of the
local contributions associated to the elements $E$ of $\mathcal{T}_{h}$. For
simplicity of notation, here and in what follows, we assume that such
splitting is considered whenever necessary. Moreover, we assume that the
approximated bilinear form $\mathcal{B}_{\kappa,h}$ is well defined on the
space $H^{1}(\mathcal{T}_{h})$.

Assumptions ($H2.{a}$) and ($H2.{b}$) allow to prove the following consistency
result.

######  Lemma 4.2.

Let $\hat{v}:=\left(v,\mu\right)\in H^{s+1}(\Omega)\times
H^{s-\nicefrac{{1}}{{2}}}(\Gamma)$, $1\leq s\leq k$, and
$\hat{v}_{h}^{I}:=\left(v_{h}^{I},\mu_{h}^{I}\right)\in V_{h}^{k}$ the
interpolant of $\hat{v}$ in $V_{h}^{k}$ such that relation (4.7) holds. Then

$\left|\mathcal{B}_{\kappa}(\hat{v}_{h}^{I},\hat{w}_{h})-\mathcal{B}_{\kappa,h}(\hat{v}_{h}^{I},\hat{w}_{h})\right|\apprle
h^{s}\left\|v\right\|_{H^{s+1}(\Omega)}\left\|\hat{w}_{h}\right\|_{V}\qquad\forall\,\hat{w}_{h}\in
V_{h}^{k}.$

###### Proof.

By abuse of notation, let
$\Pi^{\nabla}_{k}\hat{v}:=(\Pi^{\nabla}_{k}v,\mu_{h}^{I})$. We start from the
inequality

$\left|\mathcal{B}_{\kappa}(\hat{v}_{h}^{I},\hat{w}_{h})-\mathcal{B}_{\kappa,h}(\hat{v}_{h}^{I},\hat{w}_{h})\right|\leq\left|\mathcal{B}_{\kappa}(\hat{v}_{h}^{I},\hat{w}_{h})-\mathcal{B}_{\kappa,h}(\Pi_{k}^{\nabla}\hat{v},\hat{w}_{h})\right|+\left|\mathcal{B}_{\kappa,h}(\Pi_{k}^{\nabla}\hat{v},\hat{w}_{h})-\mathcal{B}_{\kappa,h}(\hat{v}_{h}^{I},\hat{w}_{h})\right|=:I+II.$
(4.9)

By using Assumption ($H2.{a}$) and the continuity of $\mathcal{B}_{\kappa}$ in
the $V$-norm, we get

$I=\left|\mathcal{B}_{\kappa}(\Pi_{k}^{\nabla}\hat{v}-\hat{v}_{h}^{I},\hat{w}_{h})\right|\apprle\left\|\Pi_{k}^{\nabla}\hat{v}-\hat{v}_{h}^{I}\right\|_{H^{1}(\mathcal{T}_{h})\times
H^{-\nicefrac{{1}}{{2}}}(\Gamma)}\left\|\hat{w}_{h}\right\|_{V}.$ (4.10)

Concerning the term $II$, from ($H2.{b}$) we obtain

$II=\left|\mathcal{B}_{\kappa,h}(\Pi_{k}^{\nabla}\hat{v}-\hat{v}_{h}^{I},\hat{w}_{h})\right|\apprle\left\|\Pi_{k}^{\nabla}\hat{v}-\hat{v}_{h}^{I}\right\|_{H^{1}(\mathcal{T}_{h})\times
H^{-\nicefrac{{1}}{{2}}}(\Gamma)}\left\|\hat{w}_{h}\right\|_{V}.$ (4.11)

Finally, by definition of $\Pi_{k}^{\nabla}\hat{v}$, using (4.5) and
Assumption ($H1.{a}$) we can write

$\left\|\Pi_{k}^{\nabla}\hat{v}-\hat{v}_{h}^{I}\right\|_{H^{1}(\mathcal{T}_{h})\times
H^{-\nicefrac{{1}}{{2}}}(\Gamma)}=\left\|\Pi_{k}^{\nabla}v-v_{h}^{I}\right\|_{H^{1}(\mathcal{T}_{h})}\leq\left\|\Pi_{k}^{\nabla}v-v\right\|_{H^{1}(\mathcal{T}_{h})}+\left\|v-v_{h}^{I}\right\|_{H^{1}(\Omega)}\apprle
h^{s}\left\|v\right\|_{H^{s+1}(\Omega)},$

from which, combining (4.10) and (4.11) with (4.9), the thesis follows. ∎

Similarly, the following lemma can be proved.

######  Lemma 4.3.

Let $\hat{v}_{\scriptscriptstyle 0}:=\left(v,\mu_{\scriptscriptstyle
0}\right)\in H^{s+1}(\Omega)\times(H^{s-\nicefrac{{1}}{{2}}}(\Gamma)\cap
H_{0}^{-\nicefrac{{1}}{{2}}}(\Gamma))$, $1\leq s\leq k$, and
$\hat{v}_{\scriptscriptstyle 0h}^{I}:=\left(v_{h}^{I},\mu_{\scriptscriptstyle
0h}^{I}\right)\in\tilde{V}_{h}^{k}$ the interpolant of
$\hat{v}_{\scriptscriptstyle 0}$ in $\tilde{V}_{h}^{k}$ such that relation
(4.8) holds. Then

$\left|\mathcal{B}_{\scriptscriptstyle 0}(\hat{v}_{\scriptscriptstyle
0h}^{I},\hat{w}_{\scriptscriptstyle 0h})-\mathcal{B}_{\scriptscriptstyle
0,h}(\hat{v}_{\scriptscriptstyle 0h}^{I},\hat{w}_{\scriptscriptstyle
0h})\right|\apprle
h^{s}\left\|v\right\|_{H^{s+1}(\Omega)}\left\|\hat{w}_{\scriptscriptstyle
0h}\right\|_{V}\qquad\forall\,\hat{w}_{\scriptscriptstyle
0h}\in\tilde{V}_{h}^{k}.$

In order to prove the main results of our theoretical analysis, we need to
introduce further assumptions on the approximate bilinear form
$\mathcal{B}_{\kappa,h}$. Denoting by
$\mathcal{D}_{\kappa,h}:=\mathcal{B}_{\kappa,h}-\mathcal{B}_{0,h}$, we
require:

  1. ($H3.{a}$)

$\mathcal{D}_{\kappa,h}$ is continuous in the weaker $W$-norm, with
$W:=L^{2}(\Omega)\times H^{-\nicefrac{{1}}{{2}}}(\Gamma)$:

$\left|\mathcal{D}_{\kappa,h}(\hat{v}_{h},\hat{w}_{h})\right|\apprle\left\|\hat{v}_{h}\right\|_{W}\left\|\hat{w}_{h}\right\|_{W}\qquad\forall\,\hat{v}_{h},\hat{w}_{h}\in
V_{h}^{k};$

  2. ($H3.{b}$)

$\mathcal{B}_{0,h}$ is $\tilde{V}_{h}^{k}$-elliptic:

$\mathcal{B}_{0,h}(\hat{w}_{\scriptscriptstyle 0h},\hat{w}_{\scriptscriptstyle
0h})\apprge\left\|\hat{w}_{\scriptscriptstyle
0h}\right\|_{V}^{2}\qquad\forall\,\hat{w}_{\scriptscriptstyle
0h}\in\tilde{V}_{h}^{k};$

  3. ($H3.{c}$)

$k$-consistency in the second term of $\mathcal{B}_{0,h}$:

$\mathcal{B}_{0,h}(\hat{w}_{h},\hat{q})=\mathcal{B}_{0}(\hat{w}_{h},\hat{q})\quad\forall\,\hat{q}\in
P_{k}(\mathcal{T}_{h})\times X_{h}^{k}\ \text{and}\ \forall\,\hat{w}_{h}\in
V_{h}^{k}.$

######  Remark 5.

We remark that Assumption ($H3.{a}$) is the discrete counterpart of the
continuity property of the bilinear form
$\mathcal{D}_{\kappa}:=\mathcal{B}_{\kappa}-\mathcal{B}_{0}$. Indeed,
according to the continuity of $\text{V}_{\kappa}-\text{V}_{0}$ and using the
Cauchy-Schwarz inequality, we obtain: for $\hat{v}=(v,\mu)\in V$ and
$\hat{w}=(w,\nu)\in V$

$\left|\mathcal{D}_{\kappa}(\hat{v},\hat{w})\right|=\left|\kappa^{2}m(v,w)-2\langle\nu,(\text{V}_{\kappa}-\text{V}_{0})\mu\rangle\right|\apprle\left\|\hat{v}\right\|_{W}\left\|\hat{w}\right\|_{W}.$
(4.12)

Assumptions ($H3.{a}$)–($H3.{c}$) are used to prove the following Lemmas 4.4,
4.5 and 4.6, which are then crucial to obtain the Ladyzhenskaya-
Babu$\check{\text{s}}$ka-Brezzi condition for the discrete bilinear form
$\mathcal{A}_{\kappa,h}$.

######  Lemma 4.4.

Let $\hat{v}_{h}=(v_{h},\mu_{h})\in\left(H^{s+1}(\Omega)\times
H^{-\nicefrac{{1}}{{2}}}(\Gamma)\right)\cap V_{h}^{k}$ with $0\leq s\leq k$.
Then

$\left|\mathcal{D}_{\kappa}(\hat{v}_{h},\hat{w}_{h})-\mathcal{D}_{\kappa,h}(\hat{v}_{h},\hat{w}_{h})\right|\apprle
h^{s+1}\left\|v_{h}\right\|_{H^{s+1}(\Omega)}\left\|\hat{w}_{h}\right\|_{W},\quad\forall\,\hat{w}_{h}\in
V_{h}^{k}.$ (4.13)

###### Proof.

Let us denote, by abuse of notation,
$\Pi_{k}^{0}\hat{v}_{h}:=\left(\Pi_{k}^{0}v_{h},\mu_{h}\right)$. By adding and
subtracting the term
$\mathcal{D}_{\kappa,h}(\Pi_{k}^{0}\hat{v}_{h},\hat{w}_{h})$ and using
Assumption ($H2.{a}$), we get

$\displaystyle\left|\mathcal{D}_{\kappa}(\hat{v}_{h},\hat{w}_{h})-\mathcal{D}_{\kappa,h}(\hat{v}_{h},\hat{w}_{h})\right|$
$\displaystyle\leq\left|\mathcal{D}_{\kappa}(\hat{v}_{h}-\Pi_{k}^{0}\hat{v}_{h},\hat{w}_{h})\right|+\left|\mathcal{D}_{\kappa,h}(\hat{v}_{h}-\Pi_{k}^{0}\hat{v}_{h},\hat{w}_{h})\right|$
$\displaystyle\apprle\left\|\hat{v}_{h}-\Pi_{k}^{0}\hat{v}_{h}\right\|_{W}\left\|\hat{w}_{h}\right\|_{W},$

the last inequality directly following from (4.12) and ($H3.{a}$). Finally, by
definition of $\Pi_{k}^{0}\hat{v}_{h}$ and using (4.4), we can write

$\left\|\hat{v}_{h}-\Pi_{k}^{0}\hat{v}_{h}\right\|_{W}=\left\|v_{h}-\Pi_{k}^{0}v_{h}\right\|_{L^{2}(\Omega)}\apprle
h^{s+1}\left\|v_{h}\right\|_{H^{s+1}(\Omega)},$

from which the thesis follows. ∎

######  Lemma 4.5.

Let $\hat{v}_{\scriptscriptstyle 0}=(v,\mu_{\scriptscriptstyle
0})\in\tilde{V}$. There exists one and only one $\hat{v}_{\scriptscriptstyle
0h}=(v_{h},\mu_{\scriptscriptstyle 0h})\in\tilde{V}_{h}^{k}$ such that

$\mathcal{B}_{0,h}(\hat{w}_{\scriptscriptstyle 0h},\hat{v}_{\scriptscriptstyle
0h})=\mathcal{B}_{0}(\hat{w}_{\scriptscriptstyle
0h},\hat{v}_{\scriptscriptstyle 0})\quad\forall\,\hat{w}_{\scriptscriptstyle
0h}=(w_{h},\nu_{\scriptscriptstyle 0h})\in\tilde{V}_{h}^{k}.$ (4.14)

Moreover, it holds:

$\displaystyle\left\|{\hat{v}}_{\scriptscriptstyle
0h}\right\|_{V}\apprle\left\|{\hat{v}_{\scriptscriptstyle 0}}\right\|_{V},$
(4.15a) $\displaystyle\left\|\mu_{\scriptscriptstyle
0h}-\mu_{\scriptscriptstyle
0}\right\|_{H^{-\nicefrac{{3}}{{2}}}(\Gamma)}\apprle
h\left\|\hat{v}_{\scriptscriptstyle 0}\right\|_{V},$ (4.15b)
$\displaystyle\left\|v_{h}-v\right\|_{L^{2}(\Omega)}\apprle
h\left\|\hat{v}_{\scriptscriptstyle 0}\right\|_{V}.$ (4.15c)

###### Proof.

Existence and uniqueness of $\hat{v}_{\scriptscriptstyle
0h}\in\tilde{V}_{h}^{k}$, solution of (4.14), follow from Assumptions
($H2.{b}$) and ($H3.{b}$). Moreover, (4.15a) holds according to ($H3.{b}$) and
the continuity of the bilinear form $\mathcal{B}_{0}$ (see Remark 3).

In order to prove (4.15b), by using a duality argument, it is sufficient to
show that:

$\left|\langle\mu_{\scriptscriptstyle 0h}-\mu_{\scriptscriptstyle
0},\eta\rangle_{H^{-\nicefrac{{3}}{{2}}}(\Gamma)\times
H^{\nicefrac{{3}}{{2}}}(\Gamma)}\right|\apprle
h\left\|\hat{v}_{\scriptscriptstyle
0h}\right\|_{V}\left\|\eta\right\|_{H^{\nicefrac{{3}}{{2}}}(\Gamma)}\quad\forall\eta\in
H^{\nicefrac{{3}}{{2}}}(\Gamma).$ (4.16)

Starting from $\eta\in H^{\nicefrac{{3}}{{2}}}(\Gamma)$, we consider
$\widetilde{w}:=\left(0,0,\eta\right)\in L^{2}(\Omega)\times
H^{\nicefrac{{1}}{{2}}}(\Gamma)\times H^{\nicefrac{{3}}{{2}}}(\Gamma)$ and we
set $\hat{w}_{\scriptscriptstyle 0}:=\mathcal{B}_{0}^{-1}\widetilde{w}$. Then,
we have:

$\mathcal{B}_{0}(\hat{w}_{\scriptscriptstyle 0},\hat{z}_{\scriptscriptstyle
0})=\mathcal{B}_{0}(\mathcal{B}_{0}^{-1}\widetilde{w},\hat{z}_{\scriptscriptstyle
0})=(\mathcal{B}_{0}\mathcal{B}_{0}^{-1}\widetilde{w})(\hat{z}_{\scriptscriptstyle
0})=\langle\zeta_{\scriptscriptstyle
0},\eta\rangle_{H^{-\nicefrac{{3}}{{2}}}(\Gamma)\times
H^{\nicefrac{{3}}{{2}}}(\Gamma)}\quad\forall\,\hat{z}_{\scriptscriptstyle
0}=(z,\zeta_{\scriptscriptstyle 0})\in\tilde{V}.$ (4.17)

From the continuity of $\mathcal{B}_{0}^{-1}:L^{2}(\Omega)\times
H^{\nicefrac{{1}}{{2}}}(\Gamma)\times
H^{\nicefrac{{3}}{{2}}}(\Gamma)\rightarrow H^{2}(\Omega)\times
H^{\nicefrac{{1}}{{2}}}(\Gamma)$ (see Remark 3), it follows that:

$\left\|\hat{w}_{\scriptscriptstyle 0}\right\|_{H^{2}(\Omega)\times
H^{\nicefrac{{1}}{{2}}}(\Gamma)}\apprle\left\|\eta\right\|_{H^{\nicefrac{{3}}{{2}}}(\Gamma)}.$
(4.18)

Therefore, by choosing $\hat{z}_{\scriptscriptstyle
0}=\hat{v}_{\scriptscriptstyle 0h}-\hat{v}_{\scriptscriptstyle 0}$ in (4.17),
we can write

$\langle\mu_{\scriptscriptstyle 0h}-\mu_{\scriptscriptstyle
0},\eta\rangle_{H^{-\nicefrac{{3}}{{2}}}(\Gamma)\times
H^{\nicefrac{{3}}{{2}}}(\Gamma)}=\mathcal{B}_{0}(\hat{w}_{\scriptscriptstyle
0},\hat{v}_{\scriptscriptstyle 0h}-\hat{v}_{\scriptscriptstyle
0})=\mathcal{B}_{0}(\hat{w}_{\scriptscriptstyle 0}-\hat{w}_{\scriptscriptstyle
0h},\hat{v}_{\scriptscriptstyle 0h}-\hat{v}_{\scriptscriptstyle
0})+\mathcal{B}_{0}(\hat{w}_{\scriptscriptstyle
0h},\hat{v}_{\scriptscriptstyle
0h})-\mathcal{B}_{0}(\hat{w}_{\scriptscriptstyle
0h},\hat{v}_{\scriptscriptstyle 0}).$

Since $\hat{v}_{\scriptscriptstyle 0h}\in\tilde{V}_{h}^{k}$ is the solution of
(4.14), we rewrite the previous identity as follows:

$\displaystyle\left|\langle\mu_{\scriptscriptstyle 0h}-\mu_{\scriptscriptstyle
0},\eta\rangle_{H^{-\nicefrac{{3}}{{2}}}(\Gamma)\times
H^{\nicefrac{{3}}{{2}}}(\Gamma)}\right|=\left|\mathcal{B}_{0}(\hat{w}_{\scriptscriptstyle
0}-\hat{w}_{\scriptscriptstyle 0h},\hat{v}_{\scriptscriptstyle
0h}-\hat{v}_{\scriptscriptstyle
0})+\mathcal{B}_{0}(\hat{w}_{\scriptscriptstyle
0h},\hat{v}_{\scriptscriptstyle
0h})-\mathcal{B}_{0,h}(\hat{w}_{\scriptscriptstyle
0h},\hat{v}_{\scriptscriptstyle 0h})\right|$
$\displaystyle\leq\left|\mathcal{B}_{0}(\hat{w}_{\scriptscriptstyle
0}-\hat{w}_{\scriptscriptstyle 0h},\hat{v}_{\scriptscriptstyle
0h}-\hat{v}_{\scriptscriptstyle
0})\right|+\left|\mathcal{B}_{0}(\hat{w}_{\scriptscriptstyle
0h},\hat{v}_{\scriptscriptstyle
0h})-\mathcal{B}_{0,h}(\hat{w}_{\scriptscriptstyle
0h},\hat{v}_{\scriptscriptstyle 0h})\right|=:I+II.$ (4.19)

Choosing in (4.1) $\hat{w}_{\scriptscriptstyle 0h}=\hat{w}_{\scriptscriptstyle
0h}^{I}$, the interpolant of $\hat{w}_{\scriptscriptstyle 0}\in\tilde{V}$ in
$\tilde{V}_{h}^{k}$ such that (4.8) holds, due to Lemma 4.3 with $s=1$ and
(4.15a), we can estimate $II$ as follows:

$II\apprle h\left\|\hat{v}_{\scriptscriptstyle
0}\right\|_{V}\left\|\hat{w}_{\scriptscriptstyle
0}\right\|_{H^{2}(\Omega)\times H^{\nicefrac{{1}}{{2}}}(\Gamma)}.$ (4.20)

Combining the continuity of $\mathcal{B}_{0}$, (4.8) with $s=1$ and (4.15a),
we have:

$I\apprle\left\|\hat{v}_{\scriptscriptstyle 0h}-\hat{v}_{\scriptscriptstyle
0}\right\|_{V}\left\|\hat{w}_{\scriptscriptstyle
0}-\hat{w}_{\scriptscriptstyle 0h}^{I}\right\|_{V}\apprle
h\left\|\hat{v}_{\scriptscriptstyle 0h}-\hat{v}_{\scriptscriptstyle
0}\right\|_{V}\left\|\hat{w}_{\scriptscriptstyle
0}\right\|_{H^{2}(\Omega)\times H^{\nicefrac{{1}}{{2}}}(\Gamma)}\apprle
h\left\|\hat{v}_{\scriptscriptstyle
0}\right\|_{V}\left\|\hat{w}_{\scriptscriptstyle
0}\right\|_{H^{2}(\Omega)\times H^{\nicefrac{{1}}{{2}}}(\Gamma)}.$ (4.21)

Finally, from (4.18) and (4.1) together with (4.20) and (4.21), inequality
(4.16) directly follows.

Inequality (4.15c) can be proved similarly to (4.15b). Indeed, if we consider
$\widetilde{w}:=\left(v_{h}-v,0,0\right)\in L^{2}(\Omega)\times
H^{\nicefrac{{1}}{{2}}}(\Gamma)\times H^{\nicefrac{{3}}{{2}}}(\Gamma)$ and
$\hat{w}_{\scriptscriptstyle 0}=\mathcal{B}_{0}^{-1}\widetilde{w}$ in (4.17),
we get

$\mathcal{B}_{0}(\hat{w}_{\scriptscriptstyle 0},\hat{z}_{\scriptscriptstyle
0})=\left(v_{h}-v,z\right)_{L^{2}(\Omega)}\qquad\forall\,\hat{z}_{\scriptscriptstyle
0}=(z,\zeta_{\scriptscriptstyle 0})\in\tilde{V}.$ (4.22)

Then, choosing $\hat{z}_{\scriptscriptstyle 0}=\hat{v}_{\scriptscriptstyle
0h}-\hat{v}_{\scriptscriptstyle 0}$ in (4.22), we have

$\mathcal{B}_{0}(\hat{w}_{\scriptscriptstyle 0},\hat{v}_{\scriptscriptstyle
0h}-\hat{v}_{\scriptscriptstyle
0})=\left(v_{h}-v,v_{h}-v\right)_{L^{2}(\Omega)}=\left\|v_{h}-v\right\|_{L^{2}(\Omega)}^{2}.$

Finally, by taking into account the continuity of $\mathcal{B}_{0}^{-1}$, we
obtain

$\left\|\hat{w}_{\scriptscriptstyle 0}\right\|_{H^{2}(\Omega)\times
H^{\nicefrac{{1}}{{2}}}(\Gamma)}\apprle\left\|v_{h}-v\right\|_{L^{2}(\Omega)}$

and, proceeding as we did to estimate (4.1), we write

$\left\|v_{h}-v\right\|_{L^{2}(\Omega)}^{2}\apprle
h\left\|\hat{w}_{\scriptscriptstyle 0}\right\|_{H^{2}(\Omega)\times
H^{\nicefrac{{1}}{{2}}}(\Gamma)}\left\|\hat{v}_{\scriptscriptstyle
0}\right\|_{V}\apprle
h\left\|v_{h}-v\right\|_{L^{2}(\Omega)}\left\|\hat{v}_{\scriptscriptstyle
0}\right\|_{V},$

from which (4.15c) follows. ∎

######  Lemma 4.6.

Let $\hat{v}=(v,\mu)\in V$. There exists $\hat{v}_{h}=(v_{h},\mu_{h})\in
V_{h}^{k}$ such that

$\mathcal{B}_{0,h}(\hat{w}_{h},\hat{v}_{h})=\mathcal{B}_{0}(\hat{w}_{h},\hat{v})+\mathcal{B}_{0}((0,\bar{\eta}_{h}),\hat{v}_{h}-\hat{v})\quad\forall\hat{w}_{h}\in(w_{h},\eta_{h})\in
V_{h}^{k}$

where $\bar{\eta}_{h}=\nicefrac{{\langle\eta_{h},1\rangle}}{{\langle
1,1\rangle}}$. Moreover, it holds

$\displaystyle\lVert\hat{v}_{h}\rVert_{V}\apprle\left\|\hat{v}\right\|_{V},$
(4.23a)
$\displaystyle\lVert\mu_{h}-\mu\rVert_{H^{-\nicefrac{{3}}{{2}}}(\Gamma)}\apprle
h\lVert\hat{v}\rVert_{V},$ (4.23b)
$\displaystyle\left\|v_{h}-v\right\|_{L^{2}(\Omega)}\apprle
h\lVert\hat{v}\rVert_{V}.$ (4.23c)

###### Proof.

Let consider $\hat{v}_{\scriptscriptstyle 0}=(v,\mu_{0}):=(v,\mu-\bar{\mu})$,
with $\bar{\mu}=\nicefrac{{\langle\mu,1\rangle}}{{\langle 1,1\rangle}}$, and
$\hat{w}_{\scriptscriptstyle
0h}:=(w_{h},\eta_{h}-\bar{\eta}_{h})\in\tilde{V}_{h}^{k}$. Then we have

$\mathcal{B}_{0}(\hat{w}_{h},\hat{v})=\mathcal{B}_{0}(\hat{w}_{\scriptscriptstyle
0h},\hat{v})+\mathcal{B}_{0}((0,\bar{\eta}_{h}),\hat{v})=\mathcal{B}_{0}(\hat{w}_{\scriptscriptstyle
0h},\hat{v}_{\scriptscriptstyle
0})+\mathcal{B}_{0}((0,\bar{\eta}_{h}),\hat{v})+\mathcal{B}_{0}(\hat{w}_{\scriptscriptstyle
0h},(0,\bar{\mu})).$

According to Lemma 4.5 applied to the first term at the right hand side of the
above equality, there exists a unique $\hat{v}_{\scriptscriptstyle
0h}=(v_{h},{\mu_{\scriptscriptstyle 0h}})\in\tilde{V}_{h}^{k}$ such that

$\mathcal{B}_{0}(\hat{w}_{h},\hat{v})=\mathcal{B}_{0,h}(\hat{w}_{\scriptscriptstyle
0h},\hat{v}_{\scriptscriptstyle
0h})+\mathcal{B}_{0}((0,\bar{\eta}_{h}),\hat{v})+\mathcal{B}_{0}(\hat{w}_{\scriptscriptstyle
0h},(0,\bar{\mu})).$

Using ($H3.{c}$) with $\hat{q}=(0,\bar{\mu})$, we can write

$\mathcal{B}_{0}(\hat{w}_{h},\hat{v})=\mathcal{B}_{0,h}(\hat{w}_{\scriptscriptstyle
0h},\hat{v}_{\scriptscriptstyle
0h})+\mathcal{B}_{0}((0,\bar{\eta}_{h}),\hat{v})+\mathcal{B}_{0,h}(\hat{w}_{\scriptscriptstyle
0h},(0,\bar{\mu}))=\mathcal{B}_{0,h}(\hat{w}_{\scriptscriptstyle
0h},\hat{v}_{h})+\mathcal{B}_{0}((0,\bar{\eta}_{h}),\hat{v}),$

where we have set
$\hat{v}_{h}=(v_{h},\mu_{h}):=(v_{h},{\mu_{\scriptscriptstyle
0h}}+\bar{\mu})\in V_{h}^{k}$. Moreover, by adding and subtracting in this
latter the term $\mathcal{B}_{0}((0,\bar{\eta}_{h}),\hat{v}_{h})$ and using
($H2.{a}$), we get

$\mathcal{B}_{0}(\hat{w}_{h},\hat{v})=\mathcal{B}_{0,h}(\hat{w}_{h},\hat{v}_{h})+\mathcal{B}_{0}((0,\bar{\eta}_{h}),\hat{v}-\hat{v}_{h}).$

By applying the Cauchy-Schwarz inequality to estimate the term

$\lVert\bar{\mu}\rVert_{H^{-\nicefrac{{1}}{{2}}}(\Gamma)}\apprle\lVert\mu\rVert_{H^{-\nicefrac{{1}}{{2}}}(\Gamma)},$

and (4.15a), we prove (4.23a) as follows:

$\lVert\hat{v}_{h}\rVert_{V}=\lVert\hat{v}_{\scriptscriptstyle
0h}+(0,\bar{\mu})\rVert_{V}\apprle\lVert\hat{v}_{\scriptscriptstyle
0}\rVert_{V}+\lVert\bar{\mu}\rVert_{H^{-\nicefrac{{1}}{{2}}}(\Gamma)}\apprle\left\|\hat{v}\right\|_{V}+\lVert\mu\rVert_{H^{-\nicefrac{{1}}{{2}}}(\Gamma)}\apprle\left\|\hat{v}\right\|_{V}.$

Finally, by using (4.15b)-(4.15c) we easily prove (4.23b) and (4.23c):

$\displaystyle\lVert\mu_{h}-\mu\rVert_{H^{-\nicefrac{{3}}{{2}}}(\Gamma)}=\left\|\mu_{\scriptscriptstyle
0h}-\mu_{0}\right\|_{H^{-\nicefrac{{3}}{{2}}}(\Gamma)}\apprle
h\left\|\hat{v}_{\scriptscriptstyle 0}\right\|_{V}\apprle
h\lVert\hat{v}\rVert_{V},$
$\displaystyle\left\|v_{h}-v\right\|_{L^{2}(\Omega)}\apprle
h\left\|\hat{v}_{\scriptscriptstyle 0}\right\|_{V}\apprle
h\lVert\hat{v}\rVert_{V}.$

∎

###  4.2 Error estimate in the energy norm

In the present section we show the validity of the inf-sup condition for the
discrete bilinear form $\mathcal{A}_{\kappa,h}$, with $\kappa>0$, and we prove
that the discrete scheme has the optimal approximation order, providing for
the optimal error estimate in the $V$-norm.

######  Theorem 4.7.

Assume that $\kappa^{2}$ is not an eigenvalue of the Laplacian in $\Omega$
endowed with a Dirichlet boundary condition on $\Gamma$. Then, for $h$ small
enough,

$\sup_{\begin{subarray}{c}\hat{v}_{h}\in V_{h}\\\ \hat{v}_{h}\neq
0\end{subarray}}\frac{\mathcal{A}_{\kappa,h}(\hat{w}_{h},\hat{v}_{h})}{\lVert\hat{v}_{h}\rVert_{V}}\apprge\lVert\hat{w}_{h}\rVert_{V}\quad\forall\,\hat{w}_{h}\in
V_{h}.$

###### Proof.

Given $\hat{w}_{h}\in V_{h}^{k}$, let
$\hat{v}:=\mathcal{A}_{\kappa}^{*-1}J\hat{w}_{h}\in V$ where $J:V\to
V^{\prime}$ denotes the canonical continuous map
$(J\hat{w})(\hat{z}):=(\hat{w},\hat{z})_{V}$. Therefore we can write:

$\displaystyle\mathcal{A}_{\kappa}(\hat{z},\hat{v})$
$\displaystyle=\mathcal{A}_{\kappa}(\hat{z},\mathcal{A}_{\kappa}^{*-1}J\hat{w}_{h})=(\mathcal{A}_{\kappa}\hat{z})(\mathcal{A}_{\kappa}^{*-1}J\hat{w}_{h})=(\mathcal{A}_{\kappa}^{*}\mathcal{A}_{\kappa}^{*-1}J\hat{w}_{h})(\hat{z})$
$\displaystyle=(J\hat{w}_{h})(\hat{z})=(\hat{w}_{h},\hat{z})_{V},\qquad\forall\,\hat{z}\in
V.$ (4.24)

Moreover, according to the continuity of $\mathcal{A}_{\kappa}^{*-1}$ (see
Remark 2) and of $J$, we obtain

$\lVert\hat{v}\rVert_{V}\apprle\lVert\hat{w}_{h}\rVert_{V}.$ (4.25)

Now, by virtue of Lemma 4.6, writing $\hat{v}=(v,\mu)\in V$, there exists
$\hat{v}_{h}=(v_{h},\mu_{h})\in V_{h}^{k}$ such that

$\mathcal{B}_{0,h}(\hat{w}_{h},\hat{v}_{h})=\mathcal{B}_{0}(\hat{w}_{h},\hat{v})+\mathcal{B}_{0}((0,\bar{\eta}_{h}),\hat{v}_{h}-\hat{v})\quad\forall\,\hat{w}_{h}=(w_{h},\eta_{h})\in
V_{h}^{k}$ (4.26)

where $\bar{\eta}_{h}=\nicefrac{{\langle\eta_{h},1\rangle}}{{\langle
1,1\rangle}}$ and such that (4.23a)-(4.23c) hold. Proceeding as in Theorem 5.2
of [25], and recalling the definitions of $\mathcal{A}_{0,h}$ and
$\mathcal{A}_{\kappa,h}$ and of $\mathcal{D}_{\kappa,h}$ and
$\mathcal{D}_{\kappa}$ (see assumption ($H3.{a}$) and Remark 5), we rewrite
$\mathcal{A}_{\kappa,h}$ as follows:

$\displaystyle\mathcal{A}_{\kappa,h}$
$\displaystyle=\mathcal{A}_{0,h}+(\mathcal{A}_{k}-\mathcal{A}_{0})+(\mathcal{A}_{0}-\mathcal{A}_{0,h})+(\mathcal{A}_{\kappa,h}-\mathcal{A}_{\kappa})$
$\displaystyle=\mathcal{A}_{0,h}+(\mathcal{A}_{k}-\mathcal{A}_{0})+(\mathcal{B}_{0}-\mathcal{B}_{0,h})+(\mathcal{B}_{\kappa,h}-\mathcal{B}_{\kappa})$
$\displaystyle=\mathcal{A}_{0,h}+(\mathcal{A}_{k}-\mathcal{A}_{0})+(\mathcal{D}_{\kappa,h}-\mathcal{D}_{\kappa}).$
(4.27)

Using Lemma 4.4 with $s=0$, we have

$\left*(\mathcal{D}_{\kappa,h}-\mathcal{D}_{\kappa})(\hat{w}_{h},\hat{v}_{h})\right\lvert\apprle
h\lVert\hat{w}_{h}\rVert_{V}\lVert\hat{v}_{h}\rVert_{W}\apprle
h\lVert\hat{w}_{h}\rVert_{V}\lVert\hat{v}_{h}\rVert_{V}.$ (4.28)

Recalling (3.5a)–(3.5c) and (4.6), and using (4.26), we get:

$\displaystyle\mathcal{A}_{0,h}(\hat{w}_{h},\hat{v}_{h})$
$\displaystyle=\mathcal{B}_{0,h}(\hat{w}_{h},\hat{v}_{h})+\mathcal{K}_{0}(\hat{w}_{h},\hat{v}_{h})=\mathcal{B}_{0}(\hat{w}_{h},\hat{v})+\mathcal{B}_{0}((0,\bar{\eta}_{h}),\hat{v}_{h}-\hat{v})+\mathcal{K}_{0}(\hat{w}_{h},\hat{v}_{h})$
$\displaystyle=\mathcal{B}_{0}(\hat{w}_{h},\hat{v})+\mathcal{K}_{0}(\hat{w}_{h},\hat{v})+\mathcal{B}_{0}((0,\bar{\eta}_{h}),\hat{v}_{h}-\hat{v})+\mathcal{K}_{0}(\hat{w}_{h},\hat{v}_{h})-\mathcal{K}_{0}(\hat{w}_{h},\hat{v})$
$\displaystyle=\mathcal{A}_{0}(\hat{w}_{h},\hat{v})+\mathcal{B}_{0}((0,\bar{\eta}_{h}),\hat{v}_{h}-\hat{v})-2\langle\mu_{h}-\mu,\text{K}_{0}w_{h}\rangle_{\Gamma}.$
(4.29)

By applying the Hölder inequality and (4.23b), we can estimate the last term
in (4.29) as follows

$\displaystyle\left*\langle\mu_{h}-\mu,\text{K}_{0}w_{h}\rangle_{H^{-\nicefrac{{3}}{{2}}}(\Gamma)\times
H^{\nicefrac{{3}}{{2}}}(\Gamma)}\right\lvert$
$\displaystyle\apprle\lVert\mu_{h}-\mu\rVert_{H^{-\nicefrac{{3}}{{2}}}(\Gamma)}\lVert\text{K}_{0}w_{h}\rVert_{H^{\nicefrac{{3}}{{2}}}(\Gamma)}$
$\displaystyle\apprle
h\lVert\hat{v}\rVert_{V}\lVert\text{K}_{0}w_{h}\rVert_{H^{\nicefrac{{3}}{{2}}}(\Gamma)}.$
(4.30)

Then, using the continuity of $\text{K}_{0}:H^{\nicefrac{{1}}{{2}}}(\Gamma)\to
H^{\nicefrac{{3}}{{2}}}(\Gamma)$ (see [29], formula (2.11)) and the trace
theorem, we obtain

$\norm{\text{K}_{0}w_{h}}_{H^{\nicefrac{{3}}{{2}}}(\Gamma)}\apprle\norm{w_{h}}_{H^{\nicefrac{{1}}{{2}}}(\Gamma)}\apprle\lVert
w_{h}\rVert_{H^{1}(\Omega)}\leq\lVert\hat{w}_{h}\rVert_{V},$

and, hence, combining this latter with (4.30), it follows that

$\left*\langle\mu_{h}-\mu,\text{K}_{0}w_{h}\rangle_{H^{-\nicefrac{{3}}{{2}}}(\Gamma)\times
H^{\nicefrac{{3}}{{2}}}(\Gamma)}\right\lvert\apprle
h\lVert\hat{v}\rVert_{V}\lVert\hat{w}_{h}\rVert_{V}.$ (4.31)

Then, from (4.29) and (4.31), we obtain

$\displaystyle\mathcal{A}_{0,h}(\hat{w}_{h},\hat{v}_{h})$
$\displaystyle\apprge\mathcal{A}_{0}(\hat{w}_{h},\hat{v})+\mathcal{B}_{0}((0,\bar{\eta}_{h}),\hat{v}_{h}-\hat{v})-h\lVert\hat{v}\rVert_{V}\lVert\hat{w}_{h}\rVert_{V}.$
(4.32)

By explicitly writing

$\displaystyle\mathcal{B}_{0}((0,\bar{\eta}_{h}),\hat{v}_{h}-\hat{v})$
$\displaystyle=-\langle\bar{\eta}_{h},v_{h}-v\rangle_{\Gamma}+2\langle\mu_{h}-\mu,\text{V}_{0}\bar{\eta}_{h}\rangle_{\Gamma}$
$\displaystyle=-\bar{\eta}_{h}\langle
1,v_{h}-v\rangle_{\Gamma}+2\bar{\eta}_{h}\langle\mu_{h}-\mu,\text{V}_{0}1\rangle_{\Gamma}=:I+II,$

and using the Cauchy-Schwarz inequality to bound
$\left*\bar{\eta}_{h}\right\lvert\apprle\lVert\eta_{h}\rVert_{H^{-\nicefrac{{1}}{{2}}}}$,
we can estimate $II$ by using Hölder inequality, the continuity of
$\text{V}_{0}:H^{\nicefrac{{1}}{{2}}}(\Gamma)\to
H^{\nicefrac{{3}}{{2}}}(\Gamma)$ (see Remark 1) and (4.23b):

$\displaystyle\left*II\right\lvert$
$\displaystyle\apprle\lVert\eta_{h}\rVert_{H^{-\nicefrac{{1}}{{2}}}(\Gamma)}\lVert\mu_{h}-\mu\rVert_{H^{-\nicefrac{{3}}{{2}}}(\Gamma)}\lVert
V_{0}1\rVert_{H^{\nicefrac{{3}}{{2}}}(\Gamma)}\apprle
h\lVert\hat{v}\rVert_{V}\lVert\hat{w}_{h}\rVert_{V}.$

To estimate the term $I$, we use Hölder inequality and the trace theorem (see
e.g. Eq. (2.1) of [20]) and we obtain, for
$0<\varepsilon<\nicefrac{{1}}{{2}}$:

$\displaystyle\left*I\right\lvert$
$\displaystyle\apprle\lVert\eta_{h}\rVert_{H^{-\nicefrac{{1}}{{2}}}(\Gamma)}\left*\langle
1,v_{h}-v\rangle_{\Gamma}\right\lvert\apprle\lVert\hat{w}_{h}\rVert_{V}\left*\langle
1,v_{h}-v\rangle_{\Gamma}\right\lvert\apprle\lVert\hat{w}_{h}\rVert_{V}\lVert
v_{h}-v\rVert_{H^{\varepsilon}(\Gamma)}$
$\displaystyle\apprle\lVert\hat{w}_{h}\rVert_{V}\lVert
v_{h}-v\rVert_{H^{\nicefrac{{1}}{{2}}+\varepsilon}(\Omega)}.$

Then, using the characterization of the fractional Sobolev space
$H^{\nicefrac{{1}}{{2}}+\varepsilon}(\Omega)$ as the real interpolation
between $L^{2}(\Omega)$ and $H^{1}(\Omega)$, by a standard result concerning
the norm of real interpolation spaces (see Prop. 2.3 of [30]), it holds that
$\lVert v_{h}-v\rVert_{H^{\nicefrac{{1}}{{2}}+\varepsilon}(\Omega)}\leq\lVert
v_{h}-v\rVert^{\nicefrac{{1}}{{2}}-\varepsilon}_{L^{2}(\Omega)}\lVert
v_{h}-v\rVert^{\nicefrac{{1}}{{2}}+\varepsilon}_{H^{1}(\Omega)}$. Hence, by
applying (4.23a) and (4.23c), we finally get:

$\displaystyle\left*I\right\lvert$ $\displaystyle\apprle
h^{\nicefrac{{1}}{{2}}-\varepsilon}\lVert\hat{v}\rVert_{V}^{\nicefrac{{1}}{{2}}-\varepsilon}\lVert
v_{h}-v\rVert^{\nicefrac{{1}}{{2}}+\varepsilon}_{H^{1}(\Omega)}\lVert\hat{w}_{h}\rVert_{V}\apprle
h^{\nicefrac{{1}}{{2}}-\varepsilon}\lVert\hat{v}\rVert_{V}^{\nicefrac{{1}}{{2}}-\varepsilon}(\lVert
v\rVert_{H^{1}(\Omega)}+\lVert
v_{h}\rVert_{H^{1}(\Omega)})^{\nicefrac{{1}}{{2}}+\varepsilon}\lVert\hat{w}_{h}\rVert_{V}$
$\displaystyle\apprle
h^{\nicefrac{{1}}{{2}}-\varepsilon}\lVert\hat{v}\rVert_{V}^{\nicefrac{{1}}{{2}}-\varepsilon}(\lVert\hat{v}\rVert_{V}+\lVert\hat{v}_{h}\rVert_{V})^{\nicefrac{{1}}{{2}}+\varepsilon}\lVert\hat{w}_{h}\rVert_{V}\apprle
h^{\nicefrac{{1}}{{2}}-\varepsilon}\lVert\hat{v}\rVert_{V}^{\nicefrac{{1}}{{2}}-\varepsilon}\lVert\hat{v}\rVert_{V}^{\nicefrac{{1}}{{2}}+\varepsilon}\lVert\hat{w}_{h}\rVert_{V}=h^{\nicefrac{{1}}{{2}}-\varepsilon}\lVert\hat{v}\rVert_{V}\lVert\hat{w}_{h}\rVert_{V}.$

Combining (4.32) with the bounds for $I$ and $II$, we can write

$\displaystyle\mathcal{A}_{0,h}(\hat{w}_{h},\hat{v}_{h})$
$\displaystyle\apprge\mathcal{A}_{0}(\hat{w}_{h},\hat{v})-h^{\nicefrac{{1}}{{2}}-\varepsilon}\lVert\hat{v}\rVert_{V}\lVert\hat{w}_{h}\rVert_{V}-h\lVert\hat{v}\rVert_{V}\lVert\hat{w}_{h}\rVert_{V}$
$\displaystyle\apprge\mathcal{A}_{0}(\hat{w}_{h},\hat{v})-h^{\nicefrac{{1}}{{2}}-\varepsilon}\lVert\hat{v}\rVert_{V}\lVert\hat{w}_{h}\rVert_{V}.$
(4.33)

Starting from (4.27), using (4.28) and (4.2), it follows

$\displaystyle\mathcal{A}_{\kappa,h}(\hat{w}_{h},\hat{v}_{h})$
$\displaystyle\apprge\mathcal{A}_{0,h}(\hat{w}_{h},\hat{v}_{h})+(\mathcal{A}_{k}-\mathcal{A}_{0})(\hat{w}_{h},\hat{v}_{h})-h\lVert\hat{w}_{h}\rVert_{V}\lVert\hat{v}_{h}\rVert_{V}$
$\displaystyle\apprge\mathcal{A}_{0}(\hat{w}_{h},\hat{v})-h^{\nicefrac{{1}}{{2}}-\varepsilon}\lVert\hat{w}_{h}\rVert_{V}\lVert\hat{v}\rVert_{V}+(\mathcal{A}_{k}-\mathcal{A}_{0})(\hat{w}_{h},\hat{v}_{h})$
$\displaystyle=\mathcal{A}_{\kappa}(\hat{w}_{h},\hat{v})-h^{\nicefrac{{1}}{{2}}-\varepsilon}\lVert\hat{v}\rVert_{V}\lVert\hat{w}_{h}\rVert_{V}+(\mathcal{A}_{\kappa}-\mathcal{A}_{0})(\hat{w}_{h},\hat{v}_{h})+(\mathcal{A}_{0}-\mathcal{A}_{\kappa})(\hat{w}_{h},\hat{v})$
$\displaystyle=\lVert\hat{w}_{h}\rVert_{V}^{2}-h^{\nicefrac{{1}}{{2}}-\varepsilon}\lVert\hat{v}\rVert_{V}\lVert\hat{w}_{h}\rVert_{V}+(\mathcal{A}_{\kappa}-\mathcal{A}_{0})(\hat{w}_{h},\hat{v}_{h}-\hat{v})$
(4.34)

having used (4.24) in the last equality.

Concerning the last term in (4.2), we explicitly write:

$(\mathcal{A}_{\kappa}-\mathcal{A}_{0})(\hat{w}_{h},\hat{v}_{h}-\hat{v})=-\kappa^{2}m(w_{h},v_{h}-v)+2\langle\mu_{h}-\mu,(\text{V}_{\kappa}-\text{V}_{0})\eta_{h}-(\text{K}_{\kappa}-\text{K}_{0})w_{h}\rangle$

and, by using the continuity of $m$, the Hölder inequality and the continuity
of $\text{V}_{\kappa}-\text{V}_{0}:H^{-\nicefrac{{1}}{{2}}}(\Gamma)\to
H^{\nicefrac{{3}}{{2}}}(\Gamma)$ and of
$\text{K}_{\kappa}-\text{K}_{0}:H^{\nicefrac{{1}}{{2}}}(\Gamma)\to
H^{\nicefrac{{3}}{{2}}}(\Gamma)$ (see Lemma 2.1),

$\displaystyle\left*(\mathcal{A}_{\kappa}-\mathcal{A}_{0})(\hat{w}_{h},\hat{v}_{h}-\hat{v})\right\lvert$
$\displaystyle\apprle\lVert\hat{w}_{h}\rVert_{V}\lVert
v_{h}-v\rVert_{L^{2}(\Omega)}+\lVert\hat{w}_{h}\rVert_{V}\lVert\mu_{h}-\mu\rVert_{H^{-\nicefrac{{3}}{{2}}}(\Gamma)}\apprle
h\lVert\hat{w}_{h}\rVert_{V}\lVert\hat{v}_{h}\rVert_{V},$ (4.35)

having used, in the last bound, (4.23b) and (4.23c). Finally, combining (4.2)
with (4.35) and (4.25) we get

$\mathcal{A}_{\kappa,h}(\hat{w}_{h},\hat{v}_{h})\apprge\lVert\hat{w}_{h}\rVert_{V}^{2}-h^{\nicefrac{{1}}{{2}}-\varepsilon}\lVert\hat{w}_{h}\rVert_{V}\lVert\hat{v}\rVert_{V}\apprge(1-h^{\nicefrac{{1}}{{2}}-\varepsilon})\lVert\hat{w}_{h}\rVert_{V}^{2}$

from which, for $h$ small enough, the claim follows. ∎

We conclude this section by proving the convergence error estimate for Problem
(4.6).

######  Theorem 4.8.

Assume that $\kappa^{2}$ is not an eigenvalue of the Laplacian in $\Omega$
endowed with a Dirichlet boundary condition on $\Gamma$. Furthermore, assume
that there exist $k\in\mathbf{N}$ such that for all $1\leq s\leq k$ and
$\kappa>0$, Assumptions ($H1.{a}$)-($H1.{c}$), ($H2.{a}$), ($H2.{b}$),
($H3.{a}$)-($H3.{c}$) hold, and $\sigma:L^{2}(\Omega)\to\mathbf{R}^{+}$ such
that

  1. ($H4.{a}$)

$\left|\mathcal{L}_{f}(\hat{v}_{h})-\mathcal{L}_{f,h}(\hat{v}_{h})\right|\apprle
h^{s}\left\|\hat{v}_{h}\right\|_{V}\,\sigma(f)\qquad\forall\,\hat{v}_{h}\in
V_{h}^{k}$.

Then, for $h$ small enough, Problem (4.6) admits a unique solution
$\hat{u}_{h}\in V_{h}^{k}$ and if $\hat{u}$, solution of Problem (3.4),
satisfies $\hat{u}\in H^{s+1}(\Omega)\times
H^{s-\nicefrac{{1}}{{2}}}(\Gamma)$, it holds

$\lVert\hat{u}-\hat{u}_{h}\rVert_{V}\apprle h^{s}\left(\lVert
u\rVert_{H^{s+1}(\Omega)}+\sigma(f)\right).$

###### Proof.

Existence and uniqueness of $\hat{u}_{h}$ follows from the discrete inf-sup
condition of Theorem 4.7. Let $\hat{u}_{h}^{I}\in V_{h}^{k}$ be the
interpolant of $\hat{u}$. By virtue of Theorem 4.7 there exists
$\hat{v}_{h}^{*}=(v_{h}^{*},\mu_{h}^{*})\in V_{h}^{k}$ such that

$\lVert\hat{u}_{h}-\hat{u}_{h}^{I}\rVert_{V}\apprle\frac{\mathcal{A}_{\kappa,h}(\hat{u}_{h}-\hat{u}_{h}^{I},\hat{v}_{h}^{*})}{\lVert\hat{v}_{h}^{*}\rVert_{V}}.$

Since $\hat{u}$ and $\hat{u}_{h}$ are solution of (3.4) and (4.6)
respectively, we have

$\displaystyle\lVert\hat{u}_{h}-\hat{u}_{h}^{I}\rVert_{V}\lVert\hat{v}_{h}^{*}\rVert_{V}$
$\displaystyle\apprle\mathcal{A}_{\kappa,h}(\hat{u}_{h}-\hat{u}_{h}^{I},\hat{v}^{*}_{h})=\mathcal{A}_{\kappa,h}(\hat{u}_{h},\hat{v}^{*}_{h})-\mathcal{A}_{\kappa,h}(\hat{u}_{h}^{I},\hat{v}^{*}_{h})$
$\displaystyle=\mathcal{L}_{f,h}(\hat{v}^{*}_{h})-\mathcal{A}_{\kappa,h}(\hat{u}_{h}^{I},\hat{v}^{*}_{h})+[\mathcal{A}_{\kappa}(\hat{u},\hat{v}^{*}_{h})-\mathcal{L}_{f}(\hat{v}^{*}_{h})]$
$\displaystyle=[\mathcal{L}_{f,h}(\hat{v}^{*}_{h})-\mathcal{L}_{f}(\hat{v}^{*}_{h})]+\mathcal{A}_{\kappa}(\hat{u}-\hat{u}_{h}^{I},\hat{v}^{*}_{h})+[\mathcal{A}_{\kappa}(\hat{u}_{h}^{I},\hat{v}^{*}_{h})-\mathcal{A}_{\kappa,h}(\hat{u}_{h}^{I},\hat{v}^{*}_{h})]$
$\displaystyle=[\mathcal{L}_{f,h}(\hat{v}^{*}_{h})-\mathcal{L}_{f}(\hat{v}^{*}_{h})]+\mathcal{A}_{\kappa}(\hat{u}-\hat{u}_{h}^{I},\hat{v}^{*}_{h})+[(\mathcal{B}_{\kappa}-\mathcal{B}_{\kappa,h})(\hat{u}_{h}^{I},\hat{v}^{*}_{h})].$

Then, by using Assumption ($H4.{a}$), the continuity of $\mathcal{A}_{\kappa}$
and Lemma 4.2, we obtain

$\displaystyle\lVert\hat{u}_{h}-\hat{u}_{h}^{I}\rVert_{V}\lVert\hat{v}_{h}^{*}\rVert_{V}$
$\displaystyle\apprle
h^{s}\left\|\hat{v}^{*}_{h}\right\|_{V}\,\sigma(f)+\lVert\hat{u}-\hat{u}_{h}^{I}\rVert_{V}\lVert\hat{v}_{h}^{*}\rVert_{V}+h^{s}\lVert
u\rVert_{H^{s+1}(\Omega)}\lVert\hat{v}_{h}^{*}\rVert_{V},$

from which it easily follows

$\lVert\hat{u}-\hat{u}_{h}\rVert_{V}\leq\lVert\hat{u}-\hat{u}_{h}^{I}\rVert_{V}+\lVert\hat{u}_{h}-\hat{u}_{h}^{I}\rVert_{V}\apprle\lVert\hat{u}-\hat{u}_{h}^{I}\rVert_{V}+h^{s}\lVert
u\rVert_{H^{s+1}(\Omega)}+h^{s}\sigma(f).$ (4.36)

Finally, combining (4.7) and (4.36) we obtain the thesis. ∎

##  5 The discrete scheme

In this section we introduce the discrete CVEM-BEM scheme for the coupling
procedure (3.2). We start by briefly describing the main tools of the VEM; we
refer the interested reader to [1, 4, 8] for a deeper presentation. In what
follows, we denote by $\mathbf{V}_{1},\ldots,\mathbf{V}_{n_{\text{\tiny{E}}}}$
the $n_{\text{\tiny{E}}}$ vertices of an element $E\in\mathcal{T}_{h}$ and by
$e_{1},\ldots,e_{n_{\text{\tiny{E}}}}$ the edges of its boundary $\partial E$.
For simplicity of presentation, we assume that at most one edge is curved
while the remaining ones are straight. We identify the curved edge by
$e_{1}\subset\partial\Omega$, to which we associate a regular invertible
parametrization $\gamma_{E}:I_{E}\rightarrow e_{1}$, where
$I_{E}\subset\mathbf{R}$ is a closed interval. Furthermore, we denote by
$\mathbf{V}_{E}$, $h_{E}$ and $|E|$ the mass center, the diameter and the
Lebesgue measure of $E$, respectively. Additionally, we call
$N_{\text{\tiny{V}}}$ and $N_{\text{{e}}}$ the numbers of total vertices and
edges of $\mathcal{T}_{h}$, respectively.

In what follows we will show that all the assumptions, used to obtain the
theoretical results in Section 4.1, are satisfied.

###  5.1 The discrete spaces $Q_{h}^{k}$, $X_{h}^{k}$ and
$\widetilde{X}_{h}^{k}$: validity of Assumptions ($H1.{a}$)–($H1.{c}$)

In order to describe the discrete space $Q_{h}^{k}$, introduced in Section 4.1
in a generic setting, we preliminarily consider for each $E\in\mathcal{T}_{h}$
the following local finite dimensional _augmented_ virtual space
$\widetilde{Q}^{k}_{h}(E)$ and the local _enhanced_ virtual space
$Q^{k}_{h}(E)$ defined respectively,

$\widetilde{Q}^{k}_{h}(E):=\left\\{v_{h}\in H^{1}(E)\cap C^{0}(E)\ :\ \Delta
v_{h}\in P_{k}(E),\
v_{h}\,\raisebox{-5.0pt}{$|_{e_{1}}$}\in\widetilde{P}_{k}(e_{1}),\
v_{h}\,\raisebox{-5.0pt}{$|_{e_{i}}$}\in P_{k}(e_{i})\
\forall\,i=2,\ldots,n_{\text{\tiny{E}}}\right\\}$

and

$Q^{k}_{h}(E):=\left\\{v_{h}\in\widetilde{Q}^{k}_{h}(E)\ :\
m^{\text{\tiny{E}}}\left(\Pi_{k}^{\nabla,E}v_{h},q\right)=m^{\text{\tiny{E}}}\left(v_{h},q\right)\
\forall\,q\in P_{k}(E)/P_{k-2}(E)\right\\},$

where $\widetilde{P}_{k}(e_{1}):=\left\\{\widetilde{q}:=q\circ\gamma_{E}^{-1}\
:\ q\in P_{k}(I_{E})\right\\}$ and $P_{k}(E)/P_{k-2}(E)$ stands for the space
of all polynomials of degree $k$ on $E$ that are $L^{2}$-orthogonal to all
polynomials of degree $k-2$ on $E$.

For details on such spaces, we refer the reader to [8] (see Remark 2.6) and to
[1] (see Section 3).

It is possible to prove (see Proposition 2 in [1] and Proposition 2.2 in [8])
that the dimension of $Q^{k}_{h}(E)$ is

$n:=\dim(Q^{k}_{h}(E))=kn_{\text{\tiny{E}}}+\frac{k(k-1)}{2}$

and that a generic element $v_{h}$ of $Q^{k}_{h}(E)$ is uniquely determined by
the following $n$ conditions (see [8], Proposition 2.2):

  * •

its values at the $n_{\text{\tiny{E}}}$ vertices of $E$;

  * •

its values at the $k-1$ internal points of the $(k+1)-$point Gauss-Lobatto
quadrature rule on every straight edge $e_{2},\ldots,e_{n_{E}}\in\partial E$;

  * •

its values at the $k-1$ internal points of $e_{1}$ that are the images,
through $\gamma_{E}$, of the $(k+1)-$point Gauss-Lobatto quadrature rule on
$I_{E}$;

  * •

the internal $k(k-1)/2$ moments of $v_{h}$ against a polynomial basis
$\mathcal{M}_{k-2}(E)$ of $P_{k-2}(E)$ defined for $k\geq 2$, as:

$\frac{1}{|E|}\int\limits_{E}v_{h}(\mathbf{x})p(\mathbf{x})\,\differential\mathbf{x}\qquad\forall\,p\in\mathcal{M}_{k-2}(E)\
\text{with}\ \|p\|_{L^{\infty}(E)}\apprle 1.$ (5.1)

According to the definition of $\widetilde{Q}^{k}_{h}(E)$, it is easy to check
that $P_{0}(E)\subset Q^{k}_{h}(E)$ while, in general, $P_{k}(E)\not\subset
Q^{k}_{h}(E)$, for $k>0$. Now, choosing an arbitrary but fixed ordering of the
degrees of freedom such that these are indexed by $i=1,\cdots,n$, we introduce
as in [4] the operator $\text{dof}_{i}:Q^{k}_{h}(E)\longrightarrow\mathbf{C}$,
defined as

$\text{dof}_{i}(v_{h}):=\text{the value of the $i$-th local degree of freedom
of}\,v_{h}.$

The basis functions $\left\\{\Phi_{j}\right\\}_{j=1}^{n}$ chosen for the space
$Q^{k}_{h}(E)$ are the standard Lagrangian ones, such that

$\text{dof}_{i}(\Phi_{j})=\delta_{ij},\qquad i,j=1,\ldots,n,$

$\delta_{ij}$ being the Kronecker delta.

On the basis of the definition of the local enhanced virtual space
$Q^{k}_{h}(E)$, we are allowed to construct the global enhanced virtual space

$Q^{k}_{h}:=\left\\{v_{h}\in H^{1}_{0,\Gamma_{0}}(\Omega)\ :\
v_{h}\,\raisebox{-5.0pt}{$|_{E}$}\in Q^{k}_{h}(E)\quad\forall
E\in\mathcal{T}_{h}\right\\}.$

We remark that the global degrees of freedom for a generic element $v_{h}\in
Q^{k}_{h}$ are:

  * •

its values at each of the $\widetilde{N}_{\text{\tiny{V}}}$ vertices of
$\mathcal{T}_{h}$ that do not belong to $\Gamma_{0}$;

  * •

its values at the $k-1$ internal points of the $(k+1)-$point Gauss-Lobatto
quadrature rule on each of the $\bar{N}_{\text{\tiny{e}}}$ straight edges of
$\mathcal{T}_{h}$;

  * •

its values at the $k-1$ internal points of the
$\widetilde{N}_{\text{\tiny{e}}}$ curved edge of $\mathcal{T}_{h}$, that do
not belong to $\Gamma_{0}$ and that are the images through the parametrization
$\gamma_{E}$ of the the $(k+1)-$point Gauss-Lobatto quadrature rule on
parametric interval $I_{E}$;

  * •

its moments up to order $k-2$ in each of the $N_{\text{\tiny{E}}}$ elements of
$\mathcal{T}_{h}$, for $k\geq 2$:

$\frac{1}{|E|}\int\limits_{E}v_{h}(\mathbf{x})p(\mathbf{x})\,\differential\mathbf{x}\qquad\forall\,p\in\mathcal{M}_{k-2}(E)(E)\
\text{with}\ \|p\|_{L^{\infty}(E)}\apprle 1.$

Consequently, $Q^{k}_{h}$ has dimension

$N:=\dim(Q^{k}_{h})=\widetilde{N}_{\text{\tiny{V}}}+(k-1)(\widetilde{N}_{\text{\tiny{e}}}+\bar{N}_{\text{\tiny{e}}})+\frac{k(k-1)}{2}N_{\text{\tiny{E}}}.$
(5.2)

######  Remark 6.

We remark that the global enhanced virtual space $Q^{k}_{h}$ defined above is
slightly different from that introduced in the pioneering paper on CVEM [8],
the latter being defined for the solution of the Laplace problem. However, as
highlighted in Remark 2.6 in [8], the theoretical analysis therein contained
can be extended to our context by following the ideas of [1].

In the following lemma we prove that Assumption ($H1.{a}$) holds for the space
$Q_{h}^{k}$.

######  Lemma 5.1.

Let $v\in H^{s+1}(\Omega)$ with $\nicefrac{{1}}{{2}}<s\leq k$. Then

$\inf_{v_{h}\in Q_{h}^{k}}\|v-v_{h}\|_{H^{1}(\Omega)}\apprle
h^{s}\|v\|_{H^{s+1}(\Omega)}.$

###### Proof.

Let $E$ be an element of $\mathcal{T}_{h}$. By virtue of Theorem 3.7 in [8]
there exists $v_{h}^{I}$, interpolant of $v$ in $Q_{h}^{k}$, such that

$|v-v_{h}^{I}|_{H^{1}(E)}\apprle h_{E}^{s}\|v\|_{H^{s+1}(E)}.$

Moreover, by using the Poincaré-Friedrichs inequality (see (2.11) in [13]), we
can write

$\|v-v_{h}^{I}\|_{L^{2}(E)}\apprle
h_{E}|v-v_{h}^{I}|_{H^{1}(E)}+\left|\int_{\partial
E}\left[v(\mathbf{x})-v_{h}^{I}(\mathbf{x})\right]ds\right|\apprle
h_{E}^{s+1}\|v\|_{H^{s+1}(E)}+\int_{\partial
E}\left|v(\mathbf{x})-v_{h}^{I}(\mathbf{x})\right|ds.$

Then, by applying the Hölder inequality, Lemma 3.2 and (3.20) in [8], we can
estimate the second term at the right hand side of the above inequality as
follows

$\displaystyle\int_{\partial
E}\left|v(\mathbf{x})-v_{h}^{I}(\mathbf{x})\right|ds$
$\displaystyle=\sum_{e\subset\partial
E}\int_{e}\left|v(\mathbf{x})-v_{h}^{I}(\mathbf{x})\right|ds\leq\sum_{e\subset\partial
E}|e|^{\nicefrac{{1}}{{2}}}\|v-v_{h}^{I}\|_{L^{2}(e)}\leq
h_{E}^{\nicefrac{{1}}{{2}}}\sum_{e\subset\partial
E}\|v-v_{h}^{I}\|_{L^{2}(e)}$ $\displaystyle\apprle
h_{E}^{\nicefrac{{1}}{{2}}}\sum_{e\subset\partial
E}h_{E}^{s+\nicefrac{{1}}{{2}}}\|v\|_{H^{s+\nicefrac{{1}}{{2}}}(e)}\apprle
h_{E}^{s+1}\|v\|_{H^{s+1}(E)}.$

Combining the local bounds for the $L^{2}$-norm and for the $H^{1}$-seminorm
of $v-v_{h}^{I}$ on $E$, we obtain

$\|v-v_{h}^{I}\|_{H^{1}(\Omega)}\apprle h^{s}\|v\|_{H^{s+1}(\Omega)},$

from which the thesis easily follows. ∎

Finally, we introduce the boundary element space $X_{h}^{k}$ associated to the
artificial boundary $\Gamma$

$X_{h}^{k}:=\left\\{\lambda\in L^{2}(\Gamma)\ :\lambda_{\mkern 1.0mu\vrule
height=6.02777pt\mkern 2.0mue}\in\widetilde{P}_{k}(e),\ \forall
e\in\Gamma\right\\}\qquad\text{with}\qquad|e|<h,$

where $|e|$ denotes the length of the edge $e$. By virtue of Theorem 4.3.20 in
[36], we have that the space $X_{h}^{k}$ satisfies the interpolation property
($H1.{b}$). For what concerns the space $\widetilde{X}_{h}^{k}=X_{h}^{k}\cap
H_{0}^{-\nicefrac{{1}}{{2}}}(\Gamma)$ and the corresponding hypothesis
($H1.{c}$), we refer to (3.2b) in [29]. Moreover, a natural basis for the
space $X_{h}^{k}$ consists in the choice of the functions
$\Phi_{j_{{|_{\Gamma}}}}$, which are the restriction of $\Phi_{j}$ on
$\Gamma$.

###  5.2 The discrete bilinear forms $\mathcal{A}_{\kappa,h}$ and
$\mathcal{B}_{\kappa,h}$: validity of Assumptions ($H2.{a}$), ($H2.{b}$) and
($H3.{a}$)–($H3.{c}$)

In order to define computable discrete local bilinear forms
$a^{\text{\tiny{E}}}_{h}:Q^{k}_{h}(E)\times Q^{k}_{h}(E)\rightarrow\mathbf{C}$
and $m^{\text{\tiny{E}}}_{h}:Q^{k}_{h}(E)\times
Q^{k}_{h}(E)\rightarrow\mathbf{C}$, following [4] and by using the definition
of $\Pi_{k}^{\nabla,E}$ and $\Pi_{k}^{0,E}$, we first split
$a^{\text{\tiny{E}}}$ and $m^{\text{\tiny{E}}}$ in a part that can be computed
exactly (up to the machine precision) and in a part that will be suitably
approximated:

$\displaystyle a^{\text{\tiny{E}}}(u_{h},v_{h})$
$\displaystyle=a^{\text{\tiny{E}}}\left(\Pi_{k}^{\nabla,E}u_{h},\Pi_{k}^{\nabla,E}v_{h}\right)+a^{\text{\tiny{E}}}\left(\left(I-\Pi_{k}^{\nabla,E}\right)u_{h},\left(I-\Pi_{k}^{\nabla,E}\right)v_{h}\right)$
(5.3) $\displaystyle m^{\text{\tiny{E}}}(u_{h},v_{h})$
$\displaystyle=m^{\text{\tiny{E}}}\left(\Pi_{k}^{0,E}u_{h},\Pi_{k}^{0,E}v_{h}\right)+m^{\text{\tiny{E}}}\left(\left(I-\Pi_{k}^{0,E}\right)u_{h},\left(I-\Pi_{k}^{0,E}\right)v_{h}\right),$
(5.4)

$I$ being the identity operator. The implementation steps for the computation
of
$a^{\text{\tiny{E}}}\left(\Pi_{k}^{\nabla,E}u_{h},\Pi_{k}^{\nabla,E}v_{h}\right)$
and $m^{\text{\tiny{E}}}\left(\Pi_{k}^{0,E}u_{h},\Pi_{k}^{0,E}v_{h}\right)$
require the choice of a suitable basis of the space $P_{k}(E)$, that allows to
define in practice the projectors $\Pi_{k}^{\nabla,E}$ and $\Pi_{k}^{0,E}$. In
accordance with the standard literature on VEM (see [5], Section 3.1), we have
considered the basis of the scaled monomials, i.e.

$\mathcal{M}_{k}(E):=\left\\{p_{\bm{\alpha}}(\mathbf{x})=\left(\frac{\mathbf{x}-\mathbf{V}_{E}}{h_{E}}\right)^{\bm{\alpha}},\
\bm{\alpha}=(\alpha_{1},\alpha_{2})\ :\
|\bm{\alpha}|=\alpha_{1}+\alpha_{2}\leq k\right\\},$

where, we recall, ${\bf V}_{E}$ and $h_{E}$ denote the mass centre and the
diameter of $E$, respectively.

Following [5], the second term in (5.3) is approximated by the following
bilinear form which represents a stabilization term:

$\displaystyle
s^{\text{\tiny{E}}}\left(\left(I-\Pi_{k}^{\nabla,E}\right)u_{h},\left(I-\Pi_{k}^{\nabla,E}\right)v_{h}\right)$
$\displaystyle:=\sum\limits_{j=1}^{n}\text{dof}_{j}\left(\left(I-\Pi_{k}^{\nabla,E}\right)u_{h}\right)\text{dof}_{j}\left(\left(I-\Pi_{k}^{\nabla,E}\right)v_{h}\right).$
(5.5)

On the contrary, since the analysis of the method requires the ellipticity
property only for the bilinear form $\mathcal{B}_{0,h}$ (see Assumption
($H3.{b}$)), an analogous stabilizing term is not needed in (5.4). Therefore
we define the approximations $a^{\text{\tiny{E}}}_{h}$ and
$m^{\text{\tiny{E}}}_{h}$ of $a^{\text{\tiny{E}}}$ and $m^{\text{\tiny{E}}}$,
respectively, as follows:

$\displaystyle a^{\text{\tiny{E}}}_{h}(u_{h},v_{h})$
$\displaystyle:=a^{\text{\tiny{E}}}\left(\Pi_{k}^{\nabla,E}u_{h},\Pi_{k}^{\nabla,E}v_{h}\right)+s^{\text{\tiny{E}}}\left(\left(I-\Pi_{k}^{\nabla,E}\right)u_{h},\left(I-\Pi_{k}^{\nabla,E}\right)v_{h}\right)$
(5.6) $\displaystyle m^{\text{\tiny{E}}}_{h}(u_{h},v_{h})$
$\displaystyle:=m^{\text{\tiny{E}}}\left(\Pi_{k}^{0,E}u_{h},\Pi_{k}^{0,E}v_{h}\right).$
(5.7)

As shown in [7] and [8], the approximate bilinear forms satisfy the following
properties:

  * •

$k$-consistency: for all $v_{h}\in Q^{k}_{h}(E)$ and for all $q\in P_{k}(E)$:

$a^{\text{\tiny{E}}}_{h}(v_{h},q)=a^{\text{\tiny{E}}}(v_{h},q)\qquad\text{and}\qquad
m^{\text{\tiny{E}}}_{h}(v_{h},q)=m^{\text{\tiny{E}}}(v_{h},q);$ (5.8)

  * •

stability: for all $v_{h}\in Q^{k}_{h}(E)$:

$a^{\text{\tiny{E}}}(v_{h},v_{h})\apprle
a^{\text{\tiny{E}}}_{h}(v_{h},v_{h})\apprle
a^{\text{\tiny{E}}}(v_{h},v_{h})\qquad\text{and}\qquad
m^{\text{\tiny{E}}}_{h}(v_{h},v_{h})\apprle m^{\text{\tiny{E}}}(v_{h},v_{h}).$
(5.9)

The global approximate bilinear forms $a_{h},m_{h}:Q^{k}_{h}\times
Q^{k}_{h}\rightarrow\mathbf{C}$ are then defined by summing up the local
contributions as follows:

$a_{h}(u_{h},v_{h}):=\sum\limits_{E\in\mathcal{T}_{h}}a^{\text{\tiny{E}}}_{h}(u_{h},v_{h})\qquad\text{and}\qquad
m_{h}(u_{h},v_{h}):=\sum\limits_{E\in\mathcal{T}_{h}}m^{\text{\tiny{E}}}_{h}(u_{h},v_{h}).$

From the right hand side of (5.9), it immediately follows that

$m_{h}(v_{h},v_{h})\apprle\left\|v_{h}\right\|_{L^{2}(\Omega)}^{2}\qquad\forall\,v_{h}\in
Q^{k}_{h}$ (5.10)

while, combining (5.9) with the Poincaré-Friedrichs inequality (see (5.3.3) in
[14]), we have:

$\left\|v_{h}\right\|_{H^{1}(\Omega)}^{2}\apprle
a_{h}(v_{h},v_{h})\apprle\left\|v_{h}\right\|_{H^{1}(\Omega)}^{2}\qquad\forall
v_{h}\in Q^{k}_{h}.$ (5.11)

Moreover, the characterization of the virtual element space $Q_{h}^{k}$ and
the boundary element space $X_{h}^{k}$ allows us to formally define the
bilinear form $\mathcal{B}_{\kappa,h}:V^{k}_{h}\times
V^{k}_{h}\rightarrow\mathbf{C}$,

$\mathcal{B}_{\kappa,h}(\hat{u}_{h},\hat{v}_{h}):=a_{h}(u_{h},v_{h})-\kappa^{2}m_{h}(u_{h},v_{h})-\langle\lambda_{h},v_{h}\rangle_{\Gamma}+\langle\mu_{h},u_{h}\rangle_{\Gamma}+2\langle\mu_{h},\text{V}_{\kappa}\lambda_{h}\rangle_{\Gamma}$

for $\hat{u}_{h}=(u_{h},\lambda_{h}),\hat{v}_{h}=(v_{h},\mu_{h})\in
V^{k}_{h}$.

From the $k$-consistency of the discrete bilinear forms
$a^{\text{\tiny{E}}}_{h}$ and $m^{\text{\tiny{E}}}_{h}$ (see (5.8)), it
immediately follows that $\mathcal{B}_{\kappa,h}$ satisfies Assumption
($H2.{a}$). Furthermore, the continuity of the bilinear forms $a_{h}$ and
$m_{h}$, as well as the continuity of the boundary operator
$\text{V}_{\kappa}$, ensure the $V$-norm continuity of
$\mathcal{B}_{\kappa,h}$, i.e. Assumption ($H2.{b}$). Analogously, the
$W$-norm continuity of
$\mathcal{D}_{\kappa,h}=\mathcal{B}_{\kappa,h}-\mathcal{B}_{0,h}$, i.e.
Assumption ($H3.{a}$), is a consequence of the continuity of $m_{h}$ and of
Lemma 2.1.

Now, to prove the $\widetilde{V}_{h}^{k}$-ellipticity, i.e. assumption
($H3.{b}$), we focus on the term

$\mathcal{B}_{0,h}(\hat{v}_{\scriptscriptstyle 0h},\hat{v}_{\scriptscriptstyle
0h})=a_{h}(v_{h},v_{h})+2\langle\mu_{\scriptscriptstyle
0h},\text{V}_{0}\mu_{\scriptscriptstyle 0h}\rangle_{\Gamma}$

for $\hat{v}_{\scriptscriptstyle 0h}=(v_{h},\mu_{\scriptscriptstyle
0h})\in\widetilde{V}^{k}_{h}$. In order to bound the first and the second term
in the above sum, we use (5.11) and Theorem 6.22 in [40], respectively, and we
get

$\mathcal{B}_{0,h}(\hat{v}_{\scriptscriptstyle 0h},\hat{v}_{\scriptscriptstyle
0h})\apprge\left\|v_{h}\right\|_{H^{1}(\Omega)}^{2}+\left\|\mu_{\scriptscriptstyle
0h}\right\|_{H^{-\nicefrac{{1}}{{2}}}(\Gamma)}^{2}=\left\|\hat{v}_{h}\right\|_{V}^{2}.$

Thus, Assumption ($H3.{c}$) is a direct consequence of the $k$-consistency
(5.8).

###  5.3 The discrete linear operator $\mathcal{L}_{f,h}$: validity of
Assumption ($H4.{a}$)

In the present section, we define the discrete linear operator
$\mathcal{L}_{f,h}:V_{h}^{k}\rightarrow\mathbf{C}$ such that

$\mathcal{L}_{f,h}(\hat{v}_{h}):=\begin{cases}\sum\limits_{E\in\mathcal{T}_{h}}m^{E}(f,\Pi^{0,E}_{1}v_{h})&k=1,2\,,\\\
\sum\limits_{E\in\mathcal{T}_{h}}m^{E}(f,\Pi^{0,E}_{k-2}v_{h})&k\geq
3.\end{cases}$

Assuming $f\in H^{k-1}(\Omega)$, in [13] (see Lemma 3.4) it has been proved
that

$\left|\mathcal{L}_{f}(\hat{v}_{h})-\mathcal{L}_{f,h}(\hat{v}_{h})\right|\apprle
h^{k}\left|f\right|_{H^{k-1}(\Omega)}\left\|v_{h}\right\|_{H^{1}(\Omega)}.$

Hence, Assumption ($H4.{a}$) is fulfilled with
$\sigma(f)=\left|f\right|_{H^{k-1}(\Omega)}$.

###  5.4 Algebraic formulation of the discrete problem

For what follows, in order to detail the algebraic form of the coupled CVEM-
BEM method, it is convenient to re-order and split the complete index set
$\mathcal{S}:=\\{j=1,\cdots,N\\}$ of the basis functions
$\left\\{\Phi_{j}\right\\}_{j\in\mathcal{S}}$ of $Q^{k}_{h}$ as

$\mathcal{S}=\mathcal{S}^{I}\cup\mathcal{S}^{\Gamma},$ (5.12)

where $\mathcal{S}^{I}$ and $\mathcal{S}^{\Gamma}$ denote the sets of indices
related to the internal degrees of freedom and to the degrees of freedom lying
on $\Gamma$, respectively. With this choice, we have

$Q_{h}^{k}=\text{span}\left\\{\Phi_{j}\right\\}_{j\in\mathcal{S}^{I}\cup\mathcal{S}^{\Gamma}},\qquad
X^{k}_{h}=\text{span}\left\\{\Phi_{j_{{|_{\Gamma}}}}\right\\}_{j\in\mathcal{S}^{\Gamma}}.$

In order to derive the linear system associated to the discrete problem (4.6),
we expand the unknown function $\hat{u}_{h}=(u_{h},\lambda_{h})\in
Q_{h}^{k}\times X^{k}_{h}$ as

$\displaystyle
u_{h}(\mathbf{x})=:\sum\limits_{j\in\mathcal{S}}u_{h}^{j}\Phi_{j}(\mathbf{x})\qquad\text{with}\quad
u_{h}^{j}=\text{dof}_{j}(u_{h})$ (5.13)
$\displaystyle\lambda_{h}(\mathbf{x})=:\sum\limits_{j\in\mathcal{S}^{\Gamma}}\lambda_{h}^{j}\Phi_{j_{{|_{\Gamma}}}}(\mathbf{x})\qquad\text{with}\quad\lambda_{h}^{j}=\text{dof}_{j}(\lambda_{h}).$

Hence, using the basis functions of $Q^{k}_{h}$ to test the discrete
counterpart of equation (3.2a), we get for
$i\in\mathcal{S}^{I}\cup\mathcal{S}^{\Gamma}$

$\displaystyle\sum\limits_{j\in\mathcal{S}^{I}\cup\mathcal{S}^{\Gamma}}u_{h}^{j}$
$\displaystyle\sum\limits_{E\in\mathcal{T}_{h}}\left[a^{\text{\tiny{E}}}_{h}(\Phi_{j},\Phi_{i})-\kappa^{2}m^{\text{\tiny{E}}}_{h}(\Phi_{j},\Phi_{i})\right]-\sum\limits_{j\in\mathcal{S}^{\Gamma}}\lambda_{h}^{j}\langle\Phi_{j},\Phi_{i}\rangle_{\Gamma}=\mathcal{L}_{f,h}((\Phi_{i},0)).$
(5.14)

To write the matrix form of the above linear system, we introduce the
stiffness and mass matrices $\mathbb{A}$, $\mathbb{M}$ and the matrix
$\mathbb{Q}$ whose entries are respectively defined by

$\mathbb{A}_{ij}:=\sum\limits_{E\in\mathcal{T}_{h}}a^{\text{\tiny{E}}}_{h}(\Phi_{j},\Phi_{i}),\qquad\mathbb{M}_{ij}:=\sum\limits_{E\in\mathcal{T}_{h}}m^{\text{\tiny{E}}}_{h}(\Phi_{j},\Phi_{i}),\qquad\mathbb{Q}_{ij}:=\langle\Phi_{j},\Phi_{i}\rangle_{\Gamma}$

and the column vectors
$\mathbf{u}=\left[u_{h}^{j}\right]_{j\in\mathcal{S}^{I}\cup\mathcal{S}^{\Gamma}}$,
$\bm{\lambda}=\left[\lambda_{h}^{j}\right]_{j\in\mathcal{S}^{\Gamma}}$ and
$\mathbf{f}=\left[\mathcal{L}_{f,h}((\Phi_{i},0))\right]_{i\in\mathcal{S}^{I}\cup\mathcal{S}^{\Gamma}}$.
In accordance with the splitting of the set of the degrees of freedom (5.12),
we consider the block partitioned representation of the above matrices and
vectors (with obvious meaning of the notation), and we rewrite equation (5.14)
as follows:

$\displaystyle\left[\begin{array}[]{ll}\mathbb{A}^{II}-\kappa^{2}\mathbb{M}^{II}&\mathbb{A}^{I\Gamma}-\kappa^{2}\mathbb{M}^{I\Gamma}\\\
&\\\ \mathbb{A}^{\Gamma I}-\kappa^{2}\mathbb{M}^{\Gamma
I}&\mathbb{A}^{\Gamma\Gamma}-\kappa^{2}\mathbb{M}^{\Gamma\Gamma}\\\
\end{array}\right]\left[\begin{array}[]{l}\mathbf{u}^{I}\\\ \\\
\mathbf{u}^{\Gamma}\end{array}\right]-\left[\begin{array}[]{l}\mathbb{Q}^{I\Gamma}\bm{\lambda}\\\
\\\
\mathbb{Q}^{\Gamma\Gamma}\bm{\lambda}\end{array}\right]=\left[\begin{array}[]{l}\mathbf{f}^{I}\\\
\\\ \mathbf{f}^{\Gamma}\end{array}\right]$ (5.27)

noting that $\mathbb{Q}^{I\Gamma}$ is a null matrix since
$\langle\Phi_{j},\Phi_{i}\rangle_{\Gamma}=0$ for $i\in\mathcal{S}^{I}$,
$j\in\mathcal{S}^{\Gamma}$.

For what concerns the discretization of the BI-NRBC, by inserting (5.13) in
(3.2b) and testing with the functions $\Phi_{i}$, $i\in\mathcal{S}^{\Gamma}$,
we obtain

$\displaystyle\sum\limits_{j\in\mathcal{S}^{\Gamma}}$
$\displaystyle\left\\{u_{h}^{j}\left[\frac{1}{2}\int\limits_{\Gamma}\Phi_{j}(\mathbf{x})\Phi_{i}(\mathbf{x})\differential\Gamma_{\mathbf{x}}-\int\limits_{\Gamma}\left(\int\limits_{\Gamma}\frac{\partial
G}{\partial\mathbf{n}_{\mathbf{y}}}(\mathbf{x},\mathbf{y})\Phi_{j}(\mathbf{y})\,\differential\Gamma_{\mathbf{y}}\right)\Phi_{i}(\mathbf{x})\differential\Gamma_{\mathbf{x}}\right]\right.$
(5.28)
$\displaystyle\,\,\left.+\lambda_{h}^{j}\int\limits_{\Gamma}\left(\int\limits_{\Gamma}G(\mathbf{x},\mathbf{y})\Phi_{j}(\mathbf{y})\,\differential\Gamma_{\mathbf{y}}\right)\Phi_{i}(\mathbf{x})\differential\Gamma_{\mathbf{x}}\right\\}=0.$

To detail the computation of the integrals in (5.28), we start by splitting
the integral on the whole $\Gamma$ into the sum of the contributions
associated to each boundary edge $\Gamma_{{}_{\ell}}$,
$\ell=1,\cdots,N^{\Gamma}$

$\displaystyle\sum\limits_{j\in\mathcal{S}^{\Gamma}}$
$\displaystyle\left\\{u_{h}^{j}\left[\frac{1}{2}\sum_{\ell=1}^{N^{\Gamma}}\int\limits_{\Gamma_{{}_{\ell}}}\Phi_{j}(\mathbf{x})\Phi_{i}(\mathbf{x})\differential\Gamma_{\mathbf{x}}-\sum_{\ell=1}^{N^{\Gamma}}\int\limits_{\Gamma_{\ell}}\left(\sum_{r=1}^{N^{\Gamma}}\int\limits_{\Gamma_{r}}\frac{\partial
G}{\partial\mathbf{n}_{\mathbf{y}}}(\mathbf{x},\mathbf{y})\Phi_{j}(\mathbf{y})\,\differential\Gamma_{\mathbf{y}}\right)\Phi_{i}(\mathbf{x})\differential\Gamma_{\mathbf{x}}\right]\right.$
(5.29)
$\displaystyle\,\,\left.+\lambda_{h}^{j}\sum_{\ell=1}^{N^{\Gamma}}\int\limits_{\Gamma_{\ell}}\left(\sum_{r=1}^{N^{\Gamma}}\int\limits_{\Gamma_{r}}G(\mathbf{x},\mathbf{y})\Phi_{j}(\mathbf{y})\,\differential\Gamma_{\mathbf{y}}\right)\Phi_{i}(\mathbf{x})\differential\Gamma_{\mathbf{x}}\right\\}=0.$

Then, denoting by $E_{\ell}$, $\ell=1,\cdots,N^{\Gamma}$, the mesh element of
$\mathcal{T}_{h}$ that has one of its curved edges on $\Gamma$ and by
$\gamma_{E_{\ell}}:I_{E_{\ell}}\rightarrow\Gamma_{\ell}$ the associated
pameterization, we rewrite (5.29) by introducing $\gamma_{E_{\ell}}$ and hence
by reducing the integration over $\Gamma_{\ell}$ to that over the parametric
interval $I_{E_{\ell}}$:

$\displaystyle\sum\limits_{j\in\mathcal{S}^{\Gamma}}$
$\displaystyle\left\\{u_{h}^{j}\left[\frac{1}{2}\sum_{\ell=1}^{N^{\Gamma}}\int\limits_{I_{E_{\ell}}}\Phi_{j}\Big{(}\gamma_{E_{\ell}}(\vartheta)\Big{)}\Phi_{i}\Big{(}\gamma_{E_{\ell}}(\vartheta)\Big{)}\Big{|}\gamma^{\prime}_{E_{\ell}}(\vartheta)\Big{|}\differential\vartheta\right.\right.$
$\displaystyle\,\,\left.\left.-\sum_{\ell=1}^{N^{\Gamma}}\int\limits_{I_{E_{\ell}}}\left(\sum_{r=1}^{N^{\Gamma}}\int\limits_{I_{E_{r}}}\frac{\partial
G}{\partial\mathbf{n}_{\mathbf{y}}}\Big{(}\gamma_{E_{\ell}}(\vartheta),\gamma_{E_{r}}(\sigma)\Big{)}\Phi_{j}\Big{(}\gamma_{E_{r}}(\sigma)\Big{)}\,\Big{|}\gamma^{\prime}_{E_{r}}(\sigma)\Big{|}\differential\sigma\right)\Phi_{i}\Big{(}\gamma_{E_{\ell}}(\vartheta)\Big{)}\Big{|}\gamma^{\prime}_{E_{\ell}}(\vartheta)\Big{|}\differential\vartheta\right]\right.$
$\displaystyle\,\,\left.+\lambda_{h}^{j}\sum_{\ell=1}^{N^{\Gamma}}\int\limits_{I_{E_{\ell}}}\left(\sum_{r=1}^{N^{\Gamma}}\int\limits_{I_{E_{r}}}G\Big{(}\gamma_{E_{\ell}}(\vartheta),\gamma_{E_{r}}(\sigma)\Big{)}\Phi_{j}\Big{(}\gamma_{E_{r}}(\sigma)\Big{)}\,\Big{|}\gamma^{\prime}_{E_{r}}(\sigma)\Big{|}\differential\sigma\right)\Phi_{i}\Big{(}\gamma_{E_{\ell}}(\vartheta)\Big{)}\Big{|}\gamma^{\prime}_{E_{\ell}}(\vartheta)\Big{|}\differential\vartheta\right\\}=0.$

Finally, introducing the BEM matrices $\mathbb{V}$ and $\mathbb{K}$ whose
entries, for $i,j\in\mathcal{S}^{\Gamma}$, are respectively

$\displaystyle\mathbb{V}_{ij}$
$\displaystyle:=\sum_{\ell=1}^{N^{\Gamma}}\int\limits_{I_{E_{\ell}}}\left(\sum_{r=1}^{N^{\Gamma}}\int\limits_{I_{E_{r}}}G\Big{(}\gamma_{E_{\ell}}(\vartheta),\gamma_{E_{r}}(\sigma)\Big{)}\Phi_{j}\Big{(}\gamma_{E_{r}}(\sigma)\Big{)}\,\Big{|}\gamma^{\prime}_{E_{r}}(\sigma)\Big{|}\differential\sigma\right)\Phi_{i}\Big{(}\gamma_{E_{\ell}}(\vartheta)\Big{)}\Big{|}\gamma^{\prime}_{E_{\ell}}(\vartheta)\Big{|}\differential\vartheta$
$\displaystyle\mathbb{K}_{ij}$
$\displaystyle:=\sum_{\ell=1}^{N^{\Gamma}}\int\limits_{I_{E_{\ell}}}\left(\sum_{r=1}^{N^{\Gamma}}\int\limits_{I_{E_{r}}}\frac{\partial
G}{\partial\mathbf{n}_{\mathbf{y}}}\Big{(}\gamma_{E_{\ell}}(\vartheta),\gamma_{E_{r}}(\sigma)\Big{)}\Phi_{j}\Big{(}\gamma_{E_{r}}(\sigma)\Big{)}\,\Big{|}\gamma^{\prime}_{E_{r}}(\sigma)\Big{|}\differential\sigma\right)\Phi_{i}\Big{(}\gamma_{E_{\ell}}(\vartheta)\Big{)}\Big{|}\gamma^{\prime}_{E_{\ell}}(\vartheta)\Big{|}\differential\vartheta,$

the matrix form of the BI-NRBC reads

$\left(\frac{1}{2}\mathbb{Q}^{\Gamma\Gamma}-\mathbb{K}\right)\mathbf{u}^{\Gamma}+\mathbb{V}\bm{\lambda}=\mathbf{0}.$
(5.30)

By combining (5.27) with (5.30) we obtain the final linear system

$\displaystyle\left[\begin{array}[]{ccc}\mathbb{A}^{II}-\kappa^{2}\mathbb{M}^{II}&\mathbb{A}^{I\Gamma}-\kappa^{2}\mathbb{M}^{I\Gamma}&\mathbb{O}\\\
&\\\ \mathbb{A}^{\Gamma I}-\kappa^{2}\mathbb{M}^{\Gamma
I}&\mathbb{A}^{\Gamma\Gamma}-\kappa^{2}\mathbb{M}^{\Gamma\Gamma}&-\mathbb{Q}^{\Gamma\Gamma}\\\
&\\\
\mathbb{O}&\frac{1}{2}\mathbb{Q}^{\Gamma\Gamma}-\mathbb{K}&\mathbb{V}\end{array}\right]\left[\begin{array}[]{l}\mathbf{u}^{I}\\\
\\\ \mathbf{u}^{\Gamma}\\\ \\\
\bm{\lambda}\end{array}\right]=\left[\begin{array}[]{c}\mathbf{f}^{I}\\\ \\\
\mathbf{f}^{\Gamma}\\\ \\\ \mathbf{0}\end{array}\right]$ (5.46)

which represents the matrix form of the coupling of equations (5.14) and
(5.29).

##  6 Numerical results

In this section, we present some numerical test to validate the theoretical
results and to show the effectiveness of the proposed method. To this aim,
some preliminary features are addressed.
For the generation of the partitioning $\mathcal{T}_{h}$ of the computational
domain $\Omega$, we have used the GMSH software (see [26]). In particular, we
have built unstructured conforming meshes consisting of quadrilaterals, by
employing the Mesh.ElementOrder option within the GMSH code. If an element $E$
has a (straight) edge bordering with a curvilinear part of $\partial\Omega$,
we transform it into a curved boundary edge by means of a suitable
parametrization. In order to develop a convergence analysis, we start by
considering the coarse mesh associated to the level of refinement zero (lev.
0) and all the successive refinements are obtained by halving each side of its
elements.
To test our numerical approach, the order $k$ of the approximation spaces is
chosen equal to 1 (linear) and 2 (quadratic) for both spaces $Q_{h}^{k}$ and
$X_{h}^{k}$. Moreover, recalling that the approximate solution $u_{h}$ is not
known inside the polygons, as suggested in [8] we compute the $H^{1}$-seminorm
and $L^{2}$-norm errors, and the corresponding Estimated Order of Convergence
(EOC), by means of the following formulas:

  * •

$H^{1}$-seminorm error
$\varepsilon^{\nabla,k}_{\text{lev}}:=\sqrt{\frac{\sum\limits_{E\in\mathcal{T}_{h}}\left|u-\Pi_{k}^{\nabla,E}u_{h}\right|^{2}_{H^{1}(E)}}{\sum\limits_{E\in\mathcal{T}_{h}}\left|u\right|^{2}_{H^{1}(E)}}}$
and
$\text{EOC}:=\log_{2}\left(\frac{\varepsilon^{\nabla,k}_{\text{lev}+1}}{\varepsilon^{\nabla,k}_{\text{lev}}}\right)$;

  * •

$L^{2}$-norm error
$\varepsilon^{0,k}_{\text{lev}}:=\sqrt{\frac{\sum\limits_{E\in\mathcal{T}_{h}}\left\|u-\Pi_{k}^{0,E}u_{h}\right\|^{2}_{L^{2}(E)}}{\sum\limits_{E\in\mathcal{T}_{h}}\left\|u\right\|^{2}_{L^{2}(E)}}}$
and
$\text{EOC}:=\log_{2}\left(\frac{\varepsilon^{0,k}_{\text{lev}+1}}{\varepsilon^{0,k}_{\text{lev}}}\right)$.

In the above formulas the superscript $k=1,2$ refers to the linear or
quadratic order approximation of $u$, the subscript lev refers to the
refinement level and, we recall, $\Pi_{k}^{\nabla,E}$ and $\Pi_{k}^{0,E}$ are
the local $H^{1}$ and $L^{2}$-projector defined in (4.1) and (4.2),
respectively.

All the numerical computations have been performed on a cluster with two
Intel® XEON® E5-2683v4 CPUs (2.1 GHz clock frequency and 16 cores) by means of
parallel Matlab® codes.

###  6.1 On the computation of the integrals involved in the proposed approach

We start by describing the quadrature techniques adopted for the computation
of the integrals appearing in the local bilinear forms $a_{h}^{E}$ and
$m_{h}^{E}$ in (5.6) and (5.7). To this aim we point out that, if the element
$E$ is a polygon with straight edges, the choice of $\mathcal{M}_{k}(E)$
defined by (5.2) allows for an exact (up to the machine precision) and easy
computation of the first term in the right hand side of (5.6) and of (5.7)
(for the details we refer to formulas (27)–(30) in [23]). On the contrary, if
the element $E$ is a polygon with a curved edge, the corresponding integrals
can be computed by applying the $n$-point Gauss-Lobatto quadrature rule. This
latter, contrarily to the former, is affected by a quadrature error if the
involved parametrization is not of polynomial type. However, in our test, the
choice $n=8$ has guaranteed the optimal convergence order of the global
scheme.

For what concerns the evaluation of the $H^{1}$-seminorm and $L^{2}$-norm
errors, to compute the associated integrals over polygons we have used the
$n$-point quadrature formulas proposed in [38] and [39], which are exact for
polynomials of degree at most $2n$. For curved polygons, we have applied the
generalization of these formulas suggested in [38] (see Remark 1). In this
case too, we have chosen $n=8$.

Finally, for what concerns the computation of the integrals defining the BEM
matrix elements, since in the theoretical analysis we have assumed that the
boundary integral operators are not approximated, it is crucial to compute
them with a high accuracy. We recall that the numerical integration
difficulties spring from the asymptotic behaviour of the Hankel function
$H_{0}^{(1)}(r)$ near the origin (see (2.6a)), the latter being the kernel of
the single layer operator $\text{V}_{\kappa}$. To compute the corresponding
integrals with high accuracy by a small number of nodes, we have used the very
simple and efficient polynomial smoothing technique proposed in [33], referred
as the _q-smoothing_ technique. It is worth noting that such technique is
applied only when the distance $r$ approaches to zero. This case corresponds
to the matrix entries belonging to the main diagonal and to those co-
diagonals, for which the supports of the basis functions overlap or are
contiguous. After having introduced the _q-smoothing_ transformation, with
$q=3$, we have applied the $n$-point Gauss-Legendre quadrature rule with $n=9$
for the outer integrals, and $n=8$ for the inner ones (see [24] and Remark 3
in [23] for further details). For the computation of all the other integrals,
we have applied a $9\times 8$-point Gauss-Legendre product quadrature rule.
Incidentally, we point out that the integrals involving the Bessel function
$H^{(1)}_{1}(r)$, appearing in the kernel of the double layer operator
$\text{K}_{\kappa}$ (see (2.6b)), do not require a particular quadrature
strategy, since its singularity $1/r$ is factored out by the same behaviour of
the Jacobian near the origin. Hence, for the computation of the entries of the
matrix $\mathbb{K}$, we have directly applied a $9\times 8$-point Gauss-
Legendre product quadrature rule.

The above described quadrature strategy has guaranteed the computation of all
the mentioned integrals with a full precision accuracy (16-digit double
precision arithmetic) for both $k=1$ and $k=2$.

###  6.2 Definition of the test problem

We consider Problem (2.1) with source $f=0$ and Dirichlet condition

$g(\mathbf{x})=\frac{\imath}{4}H_{0}^{(1)}\left(\kappa|\mathbf{x}-\mathbf{x}_{0}|\right)\quad\text{with}\quad\mathbf{x}_{0}=(0,0),\
\mathbf{x}\in\Gamma_{0}$ (6.1)

prescribed on the boundary $\Gamma_{0}$, $H_{0}^{(1)}$ being the 0-th order
Hankel function of the first kind. In this case, the exact solution
$u(\mathbf{x})$ is known and it is the field produced by the point source
$\mathbf{x}_{0}$. Its expression is given by (6.1) for every
$\mathbf{x}\in\mathbf{R}^{2}$. In the sequel, we report the numerical results
corresponding to two choices of the wave number $\kappa=1$ and $\kappa=10$.

It is worth noting that the system (5.46) is associated to the discretization
of the model problem considered for the theoretical analysis in which, we
recall, the Dirichlet boundary condition is null. In the forthcoming examples,
dealing with null source $f$ and non vanishing Dirichlet datum $g$, the right
hand side term of (5.46) involves the function $g$ instead of $f$ (see, for
example, [23] for details on the corresponding algebraic linear system).

In Example 1 we test the CVEM-BEM approach for the choice of a computational
domain having both interior and artificial curved boundaries, for which the
theoretical analysis has been given. As we will show, the use of the CVEM
reveals to be crucial to retrieve the optimal convergence order of the global
scheme. Indeed, the standard VEM approach, defined on the polygon that
approximates the curved computational domain, allows retrieving only a sub-
optimal convergence order (see Figure 5).

However, even if we have not provided theoretical results for the CVEM-BEM
coupling in the case of a non sufficiently smooth artificial boundary
$\Gamma$, in Example 2 we show that the proposed method allows to obtain the
optimal convergence order also when a polygonal $\Gamma$ is considered.

###  6.3 Example 1. Computational domain with curved boundaries

Let us consider the unbounded region $\Omega_{\text{e}}$, external to the
unitary disk $\Omega_{0}:=\\{\mathbf{x}=(x_{1},x_{2})\in\mathbf{R}^{2}\ :\
x_{1}^{2}+x_{2}^{2}\leq 1\\}$. The artificial boundary $\Gamma$ is the
circumference of radius 2,
$\Gamma=\\{\mathbf{x}=(x_{1},x_{2})\in\mathbf{R}^{2}\ :\
x_{1}^{2}+x_{2}^{2}=4\\}$, so that the finite computational domain $\Omega$ is
the annulus bounded internally by $\Gamma_{0}=\partial\Omega_{0}$ and
externally by $\Gamma$.

In Table 1, we report the total number of the degrees of freedom associated to
the CVEM space, corresponding to each decomposition level of the computational
domain, and the approximation orders $k=1,2$. To give an idea of the
curvilinear mesh used, in Figure 2 we plot the meshes corresponding to level 0
(left plot) and level 2 (right plot). We remark that the maximum level of
refinement we have considered is lev. 7 for $k=1$, whose number of degrees of
freedom coincides with that of lev. 6 for $k=2$. Therefore, in the following
tables the symbol $\times$ means that the corresponding simulation has not
been performed.

Figure 2: Meshes of $\Omega$ for lev. 0 (left plot) and lev. 2 (right plot). | lev. 0 | lev. 1 | lev. 2 | lev. 3 | lev. 4 | lev. 5 | lev. 6 | lev. 7
---|---|---|---|---|---|---|---|---
$k=1$ | $104$ | $368$ | $1,376$ | $5,312$ | $20,864$ | $82,688$ | $329,216$ | $1,313,792$
$k=2$ | $368$ | $1,376$ | $5,312$ | $20,864$ | $82,688$ | $329,216$ | $1,313,792$ | $\times$
Table 1: Example 1. Total number of degrees of freedom for $k=1,2$ and for
different levels of refinement.

In Figures 3 and 4, we show the real and imaginary parts of the numerical
solution for the wave numbers $\kappa=1$ and $\kappa=10$ respectively,
obtained by the quadratic approximation associated to the minimum refinement
level for which the graphical behaviour is accurate and not wavy; this latter
is lev. 3 for Figure 3 and lev. 5 for Figure 4.

Figure 3: Example 1. Real (left plot) and imaginary (right plot) part of the
numerical solution for $\kappa=1$, lev. 3 and $k=2$. Figure 4: Example 1. Real
(left plot) and imaginary (right plot) part of the numerical solution for
$\kappa=10$, lev. 5 and $k=2$.

In Tables 2 and 3, we report the errors $\varepsilon^{\nabla,k}_{\text{lev}}$
and $\varepsilon^{0,k}_{\text{lev}}$ and the corresponding EOC. As we can see,
for both $\kappa=1$ and $\kappa=10$, the $H^{1}$-seminorm error confirms the
convergence order $k$ of the method. Although we did not provide the
$L^{2}$-norm error estimate, the reported numerical results show the expected
convergence order $k+1$.

|  $L^{2}$-norm |  $H^{1}$-seminorm
---|---|---
lev. | $h$ | $\varepsilon^{0,1}_{\text{lev}}$ | EOC | $\varepsilon^{0,2}_{\text{lev}}$ | EOC | $\varepsilon^{\nabla,1}_{\text{lev}}$ | EOC | $\varepsilon^{\nabla,2}_{\text{lev}}$ | EOC
$0$ | $8.02e-01$ | $1.64e-02$ |  | $5.83e-04$ |  | $5.22e-02$ |  | $6.07e-03$ |
|  |  | $1.9$ |  | $3.0$ |  | $1.0$ |  | $2.0$
$1$ | $4.28e-01$ | $4.52e-03$ |  | $7.23e-05$ |  | $2.59e-02$ |  | $1.54e-03$ |
|  |  | $1.9$ |  | $3.0$ |  | $1.0$ |  | $2.0$
$2$ | $2.22e-01$ | $1.18e-03$ |  | $9.00e-06$ |  | $1.29e-02$ |  | $3.88e-04$ |
|  |  | $2.0$ |  | $3.0$ |  | $1.0$ |  | $2.0$
$3$ | $1.13e-01$ | $3.00e-04$ |  | $1.12e-06$ |  | $6.44e-03$ |  | $9.72e-05$ |
|  |  | $2.0$ |  | $3.0$ |  | $1.0$ |  | $2.0$
$4$ | $5.68e-02$ | $7.56e-05$ |  | $1.40e-07$ |  | $3.22e-03$ |  | $2.42e-05$ |
|  |  | $2.0$ |  | $3.0$ |  | $1.0$ |  | $2.0$
$5$ | $2.85e-02$ | $1.90e-05$ |  | $1.75e-08$ |  | $1.61e-03$ |  | $6.07e-06$ |
|  |  | $2.0$ |  | $3.0$ |  | $1.0$ |  | $2.0$
$6$ | $1.43e-02$ | $4.75e-06$ |  | $2.20e-09$ |  | $8.04e-04$ |  | $1.52e-06$ |
|  |  | $2.0$ |  | $\times$ |  | $1.0$ |  | $\times$
$7$ | $7.14e-03$ | $1.19e-06$ |  | $\times$ |  | $4.02e-04$ |  | $\times$ |

Table 2: Example 1. $L^{2}$-norm and $H^{1}$-seminorm relative errors and
corresponding EOC, for $\kappa=1$ and $k=1,2$.

|  $L^{2}$-norm |  $H^{1}$-seminorm
---|---|---
lev. | $h$ | $\varepsilon^{0,1}_{\text{lev}}$ | EOC | $\varepsilon^{0,2}_{\text{lev}}$ | EOC | $\varepsilon^{\nabla,1}_{\text{lev}}$ | EOC | $\varepsilon^{\nabla,2}_{\text{lev}}$ | EOC
$0$ | $8.02e-01$ | $6.03e-01$ |  | $2.57e-01$ |  | $5.77e-01$ |  | $3.07e-01$ |
|  |  | $0.8$ |  | $2.7$ |  | $0.6$ |  | $1.8$
$1$ | $4.28e-01$ | $3.52e-01$ |  | $4.00e-02$ |  | $3.92e-01$ |  | $8.59e-02$ |
|  |  | $1.5$ |  | $3.2$ |  | $1.1$ |  | $2.0$
$2$ | $2.22e-01$ | $1.33e-01$ |  | $4.37e-03$ |  | $1.84e-01$ |  | $2.18e-02$ |
|  |  | $1.8$ |  | $3.2$ |  | $1.2$ |  | $2.0$
$3$ | $1.13e-01$ | $3.76e-02$ |  | $4.71e-04$ |  | $7.88e-02$ |  | $5.49e-03$ |
|  |  | $1.9$ |  | $3.1$ |  | $1.1$ |  | $2.0$
$4$ | $5.68e-02$ | $9.74e-03$ |  | $5.51e-05$ |  | $3.65e-02$ |  | $1.38e-03$ |
|  |  | $2.0$ |  | $3.0$ |  | $1.0$ |  | $2.0$
$5$ | $2.85e-02$ | $2.46e-03$ |  | $6.75e-06$ |  | $1.78e-02$ |  | $3.44e-04$ |
|  |  | $2.0$ |  | $3.0$ |  | $1.0$ |  | $2.0$
$6$ | $1.43e-02$ | $6.16e-04$ |  | $8.39e-07$ |  | $8.86e-03$ |  | $8.61e-05$ |
|  |  | $2.0$ |  | $\times$ |  | $1.0$ |  | $\times$
$7$ | $7.14e-03$ | $1.54e-04$ |  | $\times$ |  | $4.42e-03$ |  | $\times$ |

Table 3: Example 1. $L^{2}$-norm and $H^{1}$-seminorm relative errors and
corresponding EOC, for $\kappa=10$ and $k=1,2$.

For completeness, to highlight the importance of the use of CVEM for curved
domains, in Figure 5 we report the behaviour of the $H^{1}$-seminorm and
$L^{2}$-norm errors, obtained by applying the standard VEM defined on
polygonal approximations of the computational domain. As expected, the optimal
rate of convergence is confirmed for the $H^{1}$-seminorm, while the
approximation of the geometry affects the $L^{2}$-norm error only in the case
$k=2$.

Figure 5: Example 1. $H^{1}$-seminorm (left plot) and $L^{2}$-norm (right
plot) error for a sequence of “straight” meshes, for wave number $\kappa=1$

###  6.4 Example 2. Computational domain with piece-wise linear boundaries

Let now $\Omega_{0}:=[-1,1]^{2}$ and $\Gamma$ the contour of the square
$[-2,2]^{2}$ (see Figure 6). Since the domain of interest is a polygon, we can
apply the standard (non curvilinear) VEM on polygonal meshes without
introducing an approximation of the geometry. In Table 4, we report the total
number of degrees of freedom of the VEM space, associated to each
decomposition level of the computational domain, for $k=1,2$. In Figure 6, we
plot the meshes corresponding to lev. 0 (left plot) and lev. 3 (right plot).
We remark that the maximum level of refinement we have considered is lev. 7
for k = 1, whose number of degrees of freedom coincides with that of lev. 6
for k = 2.

Figure 6: Example 2. Meshes of $\Omega$ for lev. 0 (left plot) and lev. 3 (right plot). | lev. 0 | lev. 1 | lev. 2 | lev. 3 | lev. 4 | lev. 5 | lev. 6 | lev. 7
---|---|---|---|---|---|---|---|---
$k=1$ | $120$ | $432$ | $1,632$ | $6,336$ | $24,960$ | $99,072$ | $394,752$ | $1,575,940$
$k=2$ | $432$ | $1,632$ | $6,336$ | $24,960$ | $99,072$ | $394,752$ | $1,575,940$ | $-$
Table 4: Example 2. Total number of degrees of freedom for $k=1,2$ and for
different levels of refinement.

In Figures 7 and 8, we show the real and imaginary parts of the numerical
solution for the wave numbers $\kappa=1$ and $\kappa=10$ respectively,
obtained by the quadratic approximation associated to the minimum refinement
level for which the graphical behaviour is accurate and not wavy; this latter
is lev. 3 for Figure 7 and lev. 5 for Figure 8.

Figure 7: Example 1. Real (left plot) and imaginary (right plot) part of the
numerical solution for $\kappa=1$, lev. 3 and $k=2$. Figure 8: Example 1. Real
(left plot) and imaginary (right plot) part of the numerical solution for
$\kappa=10$, lev. 5 and $k=2$.

As we can see from Tables 5 and 6, the expected convergence order of the VEM-
BEM approach for both $H^{1}$-seminorm and $L^{2}$-norm errors are confirmed,
even if the assumption on the regularity of the artificial boundary $\Gamma$,
required by the theory, is not satisfied.

|  $L^{2}$-norm |  $H^{1}$-seminorm
---|---|---
lev. | $h$ | $\varepsilon^{0,1}_{\text{lev}}$ | EOC | $\varepsilon^{0,2}_{\text{lev}}$ | EOC | $\varepsilon^{\nabla,1}_{\text{lev}}$ | EOC | $\varepsilon^{\nabla,2}_{\text{lev}}$ | EOC
$0$ | $7.60e-01$ | $1.71e-02$ |  | $8.34e-04$ |  | $1.57e-01$ |  | $1.66e-02$ |
|  |  | $2.0$ |  | $3.0$ |  | $1.1$ |  | $2.0$
$1$ | $3.85e-01$ | $4.37e-03$ |  | $1.01e-04$ |  | $7.57e-02$ |  | $4.07e-03$ |
|  |  | $2.0$ |  | $3.0$ |  | $1.0$ |  | $2.0$
$2$ | $1.94e-01$ | $1.10e-03$ |  | $1.26e-05$ |  | $3.78e-02$ |  | $1.02e-03$ |
|  |  | $2.0$ |  | $3.0$ |  | $1.0$ |  | $2.0$
$3$ | $9.73e-02$ | $2.74e-04$ |  | $1.57e-06$ |  | $1.89e-02$ |  | $2.56e-04$ |
|  |  | $2.0$ |  | $3.0$ |  | $1.0$ |  | $2.0$
$4$ | $4.87e-02$ | $6.86e-05$ |  | $1.96e-07$ |  | $9.46e-03$ |  | $6.40e-05$ |
|  |  | $2.0$ |  | $3.0$ |  | $1.0$ |  | $2.0$
$5$ | $2.44e-02$ | $1.71e-05$ |  | $2.46e-08$ |  | $4.73e-03$ |  | $1.60e-05$ |
|  |  | $2.0$ |  | $2.9$ |  | $1.0$ |  | $2.0$
$6$ | $1.22e-02$ | $4.29e-06$ |  | $3.35e-09$ |  | $2.36e-03$ |  | $4.03e-06$ |
|  |  | $2.0$ |  | $\times$ |  | $1.0$ |  | $\times$
$7$ | $6.10e-03$ | $1.07e-06$ |  | $\times$ |  | $1.18e-03$ |  | $\times$ |

Table 5: Example 2. $L^{2}$-norm and $H^{1}$-seminorm relative errors and
corresponding EOC, for $\kappa=1$ and $k=1,2$.

|  $L^{2}$-norm |  $H^{1}$-seminorm
---|---|---
lev. | $h$ | $\varepsilon^{0,1}_{\text{lev}}$ | EOC | $\varepsilon^{0,2}_{\text{lev}}$ | EOC | $\varepsilon^{\nabla,1}_{\text{lev}}$ | EOC | $\varepsilon^{\nabla,2}_{\text{lev}}$ | EOC
$0$ | $7.60e-01$ | $1.02e-00$ |  | $4.21e-01$ |  | $1.05e-00$ |  | $5.54e-01$ |
|  |  | $1.0$ |  | $3.7$ |  | $0.7$ |  | $1.7$
$1$ | $3.85e-01$ | $5.22e-01$ |  | $3.25e-02$ |  | $6.43e-01$ |  | $1.25e-01$ |
|  |  | $1.7$ |  | $3.1$ |  | $1.2$ |  | $2.0$
$2$ | $1.94e-01$ | $1.60e-01$ |  | $3.78e-03$ |  | $2.77e-01$ |  | $3.24e-02$ |
|  |  | $1.9$ |  | $3.1$ |  | $1.2$ |  | $2.0$
$3$ | $9.73e-02$ | $4.22e-02$ |  | $4.55e-04$ |  | $1.23e-01$ |  | $8.16e-03$ |
|  |  | $2.0$ |  | $3.0$ |  | $1.1$ |  | $2.0$
$4$ | $4.87e-02$ | $1.07e-02$ |  | $5.62e-05$ |  | $5.92e-02$ |  | $2.04e-03$ |
|  |  | $2.0$ |  | $3.0$ |  | $1.0$ |  | $2.0$
$5$ | $2.44e-02$ | $2.67e-03$ |  | $7.01e-06$ |  | $2.93e-02$ |  | $5.11e-04$ |
|  |  | $2.0$ |  | $3.0$ |  | $1.0$ |  | $2.0$
$6$ | $1.22e-02$ | $6.68e-04$ |  | $8.85e-07$ |  | $1.46e-02$ |  | $1.28e-04$ |
|  |  | $2.0$ |  | $\times$ |  | $1.0$ |  | $\times$
$7$ | $6.10e-03$ | $1.67e-04$ |  | $\times$ |  | $7.29e-03$ |  | $\times$ |

Table 6: Example 2. $L^{2}$-norm and $H^{1}$-seminorm relative errors and
corresponding EOC, for $\kappa=10$ and $k=1,2$.

##  7 Conclusions and perspectives

In this paper we have proposed a novel numerical approach for the solution of
2D Helmholtz problems defined in unbounded regions, external to bounded
obstacles. It consists in reducing the unbounded domain to a finite
computational one and in the coupling of the CVEM with the one equation BEM,
by means of the Galerkin approach. We have carried out the theoretical
analysis of the method in a quite abstract framework, providing an optimal
error estimate in the energy norm, under assumptions that can, in principle,
include a variety of discretization spaces wider than those we have
considered, for both the BIE and the interior PDE. While the VEM/CVEM has been
extensively and successfully applied to interior problems, its application to
exterior problems is still at an early stage and, to the best of the authors’
knowledge, the CVEM approach has been applied in this paper for the first time
to solve exterior frequency-domain wave propagation problems in the Galerkin
context.

We remark that the above mentioned coupling has been proposed in a conforming
approach context, so that the order of the CVEM and BEM approximation spaces
have been chosen with the same polynomial degree of accuracy, and the grid
used for the BEM discretization is the one inherited by the interior CVEM
decomposition.

It is worth noting that it is possible in principle to decouple the CVEM and
the BEM discretization, both in terms of degree of accuracy and of non-
matching grids. This would lead to a non-conforming coupling approach by
using, for example, a mortar type technique (see for instance [11]). Such an
approach would offer the further advantage of coupling different types of
approximation spaces and of using fast techniques for the discretization of
the BEM (see for example the very recent papers [17, 12, 21, 22]). This will
be the subject of a future investigation.

## Acknowledgments

This research benefits from the HPC (High Performance Computing) facility of
the University of Parma, Italy.

## Funding

This work was performed as part of the GNCS-IDAM 2020 research program
_“Metodologie innovative per problemi di propagazione di onde in domini
illimitati: aspetti teorici e computazionali”_. The second and the third
author were partially supported by MIUR grant _“Dipartimenti di Eccellenza
2018-2022”_ , CUP E11G18000350001.

## References

  * [1] B. Ahmad, A. Alsaedi, F. Brezzi, L. D. Marini, and A. Russo.  Equivalent projectors for virtual element methods.  Comput. Math. Appl., 66(3):376–391, 2013.
  * [2] A. Aimi, L. Desiderio, P. Fedeli, and A. Frangi.  A fast boundary-finite element approach for estimating anchor losses in micro-electro-mechanical system resonators.  Applied Mathematical Modelling, 97:741–753, 2021.
  * [3] P.F. Antonietti, L. Beirão da Veiga, S. Sacchi, and M. Verani.  A C1 virtual element method for the Cahn-Hilliard equation with polygonal meshes.  SIAM J. Numer. Anal., 54(1):34–56, 2016.
  * [4] L. Beirão da Veiga, F. Brezzi, A. Cangiani, G. Manzini, L. D. Marini, and A. Russo.  Basic principles of virtual element methods.  Math. Models Methods Appl. Sci., 23(1):199–214, 2013.
  * [5] L. Beirão da Veiga, F. Brezzi, L. D. Marini, and A. Russo.  The hitchhiker’s guide to the virtual element method.  Math. Models Methods Appl. Sci., 24(8):1541–1573, 2014.
  * [6] L. Beirão da Veiga, K. Lipnikov, and G. Manzini.  Arbitrary order nodal mimetic discretizations of elliptic problems on polygonal meshes.  SIAM J. Numer. Anal., 49(5):1737–1760, 2011.
  * [7] L. Beirão da Veiga, A. Russo, and G. Vacca.  The virtual element method with curved edges.  arXiv:1711.04306 [math.NA], 2018.
  * [8] L. Beirão da Veiga, A. Russo, and G. Vacca.  The virtual element method with curved edges.  ESAIM Math. Model. Numer. Anal., 53(2):375–404, 2019.
  * [9] E. Benvenuti, A. Chiozzi, G. Manzini, and N. Sukumar.  Extended virtual element method for the Laplace problem with singularities and discontinuities.  Comput. Methods Appl. Mech. Engrg., 356:571–597, 2019.
  * [10] S. Berrone, A. Borio, and G. Manzini.  SUPG stabilization for the nonconforming virtual element method for advection-diffusion-reaction equations.  Comput. Methods Appl. Mech. Engrg., 340:500–529, 2018.
  * [11] S. Bertoluzza and S. Falletta.  FEM solution of exterior elliptic problems with weakly enforced integral non reflecting boundary conditions.  J. Sci. Comput., 81(2):1019–1049, 2019.
  * [12] S. Bertoluzza, S. Falletta, and L. Scuderi.  Wavelets and convolution quadrature for the efficient solution of a 2D space-time BIE for the wave equation.  Appl. Math. Comput., 366:124726, 2020.
  * [13] S.C. Brenner, Q. Guan, and L.Y. Sung.  Some estimates for virtual element methods.  Comput. Methods Appl. Math., 17(4):553–574, 2017.
  * [14] S.C. Brenner and L.R. Scott.  The Mathematical Theory of Finite Element Methods, volume 15 of Texts in Applied Mathematics.  Springer, New York, third edition, 2008.
  * [15] F. Brezzi, A. Buffa, and K. Lipnikov.  Mimetic finite differences for elliptic problems.  M2AN Math. Model. Numer. Anal., 43:277–295, 2009.
  * [16] O. Certik, F. Gardini, G. Manzini, L. Mascotto, and G. Vacca.  The p- and hp-versions of the virtual element method for elliptic eigenvalue problems.  Comput. Math. Appl., 79(7):2035–2056, 2020.
  * [17] S. Chaillat, L. Desiderio, and P. Ciarlet.  Theory and implementation of $\mathcal{H}$-matrix based iterative and direct solvers for Helmholtz and elastodynamic oscillatory kernels.  J. Comput. Phys., 351:165–186, 2017.
  * [18] D. Colton and R. Kress.  Inverse Acoustic and Electromagnetic Scattering Theory, volume 93 of Applied Mathematical Sciences.  Springer, Cham, 2019.
  * [19] M. Costabel.  Symmetric methods for the coupling of finite elements and boundary elements (invited contribution).  In Boundary elements IX, Vol. 1 (Stuttgart, 1987), pages 411–420. Comput. Mech., Southampton, 1987.
  * [20] M. Costabel.  Boundary integral operators on Lipschitz domains: elementary results.  SIAM J. Math. Anal., 19(3):613–626, 1988.
  * [21] L. Desiderio.  An $\mathcal{H}$-matrix based direct solver for the Boundary Element Method in 3D elastodynamics.  AIP Conf. Proc., 1978:120005$\\_$1–120005$\\_$4, 2018.
  * [22] L. Desiderio and S. Falletta.  Efficient solution of two-dimensional wave propagation problems by CQ-wavelet BEM: algorithm and applications.  SIAM J. Sci. Comput., 42(4):B894–B920, 2020.
  * [23] L. Desiderio, S. Falletta, and L. Scuderi.  A virtual element method coupled with a boundary integral non reflecting condition for 2D exterior Helmholtz problems.  Comput. Math. Appl., 84:296–313, 2021.
  * [24] S. Falletta, G. Monegato, and L. Scuderi.  A space-time BIE method for wave equation problems: the (two-dimensional) Neumann case.  IMA J. Numer. Anal., 34(1):390–434, 2014.
  * [25] G. N. Gatica and S. Meddahi.  Coupling of virtual element and boundary element methods for the solution of acoustic scattering problems.  J. Numer. Math., 28(4):223–245, 2020.
  * [26] C. Geuzaine and J.F. Remacle.  Gmsh: a three-dimensional finite element mesh generator with built-in pre- and post processing facilities.  Internat. J. Numer. Methods Engrg., (79):1309–1331, 2009.
  * [27] H. D. Han.  A new class of variational formulations for the coupling of finite and boundary element methods.  J. Comput. Math., 8(3):223–232, 1990.
  * [28] G. C. Hsiao and W. L. Wendland.  Boundary Integral Equations, volume 164 of Applied Mathematical Sciences.  Springer-Verlag, Berlin, 2008.
  * [29] C. Johnson and J.-C. Nédélec.  On the coupling of boundary integral and finite element methods.  Math. Comp., 35(152):1063–1079, 1980.
  * [30] J.L. Lions and E. Magenes.  Problèmes aux Limites non Homogènes et Applications. II.  Travaux et Recherches Mathématiques, No. 18. Dunod, Paris, 1968.
  * [31] G. Manzini, K. Lipnikov, J.D. Moulton, and M. Shashkov.  Convergence analysis of the mimetic finite difference method for elliptic problems with staggered discretizations of diffusion coefficients.  SIAM J. Numer. Anal., 55(6):2956–2981, 2017.
  * [32] S. Márquez, A. Meddahi and V. Selgas.  Computing acoustic waves in an inhomogeneous medium of the plane by a coupling of spectral and finite elements.  SIAM J. Numer. Anal., 41(5):1729–1750, 2003.
  * [33] G. Monegato and L. Scuderi.  Numerical integration of functions with boundary singularities.  J. Comput. Appl. Math., 112(1-2):201–214, 1999.
  * [34] A. Quarteroni and A. Valli.  Numerical Approximation of Partial Differential Equations, volume 23 of Springer Series in Computational Mathematics.  Springer-Verlag, Berlin, 1994.
  * [35] J. Saranen and G. Vainikko.  Periodic Integral and Pseudodifferential Equations with Numerical Approximation.  Springer Monographs in Mathematics. Springer-Verlag, Berlin, 2002.
  * [36] S. A. Sauter and C. Schwab.  Boundary Element Methods, volume 39 of Springer Series in Computational Mathematics.  Springer-Verlag, Berlin, 2011.
  * [37] F.J. Sayas.  The validity of Johnson-Nédélec’s BEM-FEM coupling on polygonal interfaces.  SIAM J. Numer. Anal., 47(5):3451–3463, 2009.
  * [38] A. Sommariva and M. Vianello.  Product Gauss cubature over polygons based on Green’s integration formula.  BIT, 47(2):441–453, 2007.
  * [39] A. Sommariva and M. Vianello.  Gauss-Green cubature and moment computation over arbitrary geometries.  J. Comput. Appl. Math., 231(2):886–896, 2009.
  * [40] O. Steinbach.  Numerical Approximation Methods for Elliptic Boundary Value Problems.  Springer, New York, 2008.