# CVEM-BEM coupling with decoupled orders for 2D exterior Poisson problems

Luca Desiderio
Dipartimento di Scienze Matematiche, Fisiche e Informatiche
Università di Parma
Parma, 43124, Italy
<EMAIL_ADDRESS>
&Silvia Falletta
Dipartimento di Scienze Matematiche “G.L. Lagrange”
Politecnico di Torino
Torino, 10129, Italy
<EMAIL_ADDRESS>
&Matteo Ferrari
Dipartimento di Scienze Matematiche “G.L. Lagrange”
Politecnico di Torino
Torino, 10129, Italy
<EMAIL_ADDRESS>
&Letizia Scuderi
Dipartimento di Scienze Matematiche “G.L. Lagrange”
Politecnico di Torino
Torino, 10129, Italy
<EMAIL_ADDRESS>

###### Abstract

For the solution of 2D exterior Dirichlet Poisson problems we propose the
coupling of a Curved Virtual Element Method (CVEM) with a Boundary Element
Method (BEM), by using decoupled approximation orders. We provide optimal
convergence error estimates, in the energy and in the weaker
$\textit{L}^{\text{2}}$-norm, in which the CVEM and BEM contributions to the
error are separated. This allows taking advantage of the high order
flexibility of the CVEM to retrieve an accurate discrete solution by using a
low order BEM. The numerical results confirm the a priori estimates and show
the effectiveness of the proposed approach.

_Keywords_ Exterior Poisson problems, curved virtual element method, boundary
element method, coupling, error estimates.

##  1 Introduction

In this paper we deal with the following 2D problem

$\left\\{\begin{aligned} &-\Delta
u_{e}(\mathbf{x})=f(\mathbf{x})&\mathbf{x}\in\Omega_{e},\\\
&u_{e}(\mathbf{x})=0&\mathbf{x}\in\Gamma_{0},\end{aligned}\right.$ (1)

where $\Omega_{e}:=\mathbf{R}^{2}\setminus\overline{\Omega}_{{0}}$ is an
unbounded domain, exterior to an open bounded one $\Omega_{{0}}$, with
Lipschitz boundary $\Gamma_{{0}}$. It is known (see [24] and the references
therein) that Problem (1) admits a unique solution in the space

$W^{1}(\Omega_{e}):=\left\\{v:\omega v\in L^{2}(\Omega_{e}),\nabla
v\in[L^{2}(\Omega_{e})]^{2}\right\\}$

with
$\omega(\mathbf{x}):=\left(\sqrt{1+\|\mathbf{x}\|^{2}}\left(1+\log\left(\sqrt{1+\|\mathbf{x}\|^{2}}\right)\right)\right)^{-1}$,
satisfying the asymptotic conditions

$u_{e}(\mathbf{x})=\alpha+O\left(\frac{1}{\|\mathbf{x}\|}\right)\quad\text{and}\quad\nabla
u_{e}(\mathbf{x})=O\left(\frac{1}{\|\mathbf{x}\|^{2}}\right)\quad\text{for}\quad\|\mathbf{x}\|\rightarrow\infty.$
(2)

The constant $\alpha$ represents the asymptotic behaviour of $u_{e}$ at
infinity and, here, its value is not fixed in advance.

The above problem is of interest in many engineering and physical
applications, for example when studying electric and thermal plane fields on
infinite domains produced by point sources, or when solving problems of fluid
flows around obstacles. Many and various numerical methods have been proposed
and analysed for its solution, among which we mention the traditional BEM.
This latter is the most natural way to deal with unbounded domains (for a
reference, see [28] and the bibliography therein contained). Another common
approach is the coupling of a classical variational or finite difference
method with a transparent (absorbing or non-reflecting) condition defined on
an artificial boundary $\Gamma$, properly chosen to delimit a finite
computational domain. Among the most commonly used Non-Reflecting Boundary
Conditions (NRBCs), those of integral type are exact (i.e. not approximated)
and allow treating artificial boundaries of arbitrary, even non-convex,
shapes.

The aim of this paper is to propose such a coupling by means of the interior
CVEM and the one-equation BEM. Standard VEMs have been applied to a wide
variety of interior problems (see the pioneering [3] for the Poisson problem
and [1, 26, 2, 6] for more recent applications), but only few papers deal with
exterior problems (see [19, 20, 17, 16] for elliptic equations). Among the
CVEM approaches till now investigated, we mention those proposed in [10] and
[5]. Although the latter deals with local polynomial preserving VEM spaces, we
choose the former since it is well-suited for problems characterized by
computational domains with prescribed curved boundaries, like ours.

The choice of using VEM, or the more general CVEM, is mainly motivated by the
following reasons: it allows us to consider meshes whose elements can be of
general shape, and to use local discrete spaces of arbitrarily high order by
maintaining the simplicity of implementation independent of it. Moreover, the
nature of the VEM allows decoupling the approximation orders and the mesh
grids associated with the domain and boundary methods, without the need of
using special auxiliary variables (like mortar ones) for the coupling. Indeed,
by exploiting the peculiar construction of the VEM, it is possible to add
hanging nodes on the edges of the elements that belong to the artificial
boundary, without significantly modifying the structure of the interior mesh.

For what concerns the one-equation BEM, we recall that it has been proposed in
the well known Johnson & Nédélec Coupling (JNC) (see [14, 24]) and it is based
on a single Boundary Integral Equation (BIE) that involves the integral
operators associated to the fundamental solution (and its normal derivative)
of the Laplace equation.

In the recent work [19], in which a similar problem has been studied, the
authors consider the Costabel & Han Coupling (CHC) (see [15, 22]) combined
with an interior VEM. This approach yields to a symmetric and non-positive
definite scheme but, involving a BIE of hypersingular type, turns out to be
quite onerous from the computational point of view. Even if the CHC has been
applied in several contexts, the JNC turns out to be very appealing from the
engineering point of view, this latter being cheaper and easier to implement.
We remark in addition that, unlikely in [19], we deal with the asymptotic
condition (2) that entails

$\int_{\Gamma}\lambda_{e}({\mathbf{y}})\mathrm{d}\Gamma_{{\mathbf{y}}}=0,$ (3)

where $\lambda_{e}({\mathbf{y}}):=\frac{\partial
u_{e}}{\partial\mathbf{n}}(\mathbf{y})$ denotes the normal derivative of
$u_{e}$ along the artificial boundary $\Gamma$. As a consequence, suitable
spaces satisfying identity (3) have to be considered.

For the discretization of our coupled problem we consider a full Galerkin
approach based on a CVEM in the interior of the computational domain and on a
BEM associated to basis functions chosen in such a way that (3) is satisfied.
We study the proposed approach from the theoretical point of view in a quite
general framework, and we provide optimal error estimates in the energy and in
the weaker $L^{2}$-norm. In particular, since we consider here curved domains,
the use of curvilinear elements instead of polygonal ones, allows us to reach
the optimal convergence rate for degrees of accuracy higher than 2, avoiding
the sub-optimal rate caused by the approximation of the domain.

By a careful study, we show that the source of the approximation error of the
discrete solution, both in the energy and in the $L^{2}$-norm, can be split
into two contributions: a boundary (BEM) and an interior (CVEM) one. In
particular, we show that the boundary contribution behaves like
$h_{\partial}^{k_{\partial}}$ ($h_{\partial}$ denoting the maximum edge length
of the artificial boundary and $k_{\partial}$ representing the BEM polynomial
degree of accuracy), and the interior one like $h_{\circ}^{k_{\circ}}$
($h_{\circ}$ being the element diameter and $k_{\circ}$ the CVEM order
degree). Hence, for $h_{\partial}\ll h_{\circ}$ and by fixing $k_{\partial}$,
it results that the bulk error dominates the boundary one up to a certain CVEM
order, an aspect that allows obtaining a high accuracy of the global scheme
with a low BEM order.

The paper is organized as follows: in the next section we present the model
problem for the Poisson equation and its reformulation in a bounded region,
obtained by introducing the artificial boundary and its associated one
equation Boundary Integral Non Reflecting Boundary Condition (BI-NRBC). In
Section 3 we introduce the variational formulation of the problem restricted
to the finite computational domain. In Section 4 we apply the Galerkin method
and we prove error estimates in the energy and in the $L^{2}$-norm in an
abstract framework, provided that suitable hypotheses are assumed. Then we
show that these latter are satisfied by the CVEM-BEM approximation spaces
introduced in Section 5. Finally, in the last section we detail the choice of
the particular basis functions used for the approximation of the normal
derivative unknown, and we present some numerical test which confirm the
theoretical results.

##  2 The model problem

Let $\Omega_{e}:=\mathbf{R}^{2}\setminus\overline{\Omega}_{{0}}$ be an
unbounded domain, exterior to an open bounded domain
$\Omega_{{0}}\subset\mathbf{R}^{2}$, and denote by
$\Gamma_{{0}}:=\partial{\Omega}_{e}$ its Lipschitz boundary having positive
Haussdorf measure (see Figure 1 (a)). We consider the exterior Dirichlet
Poisson problem (1) in the unknown solution $u_{e}$, where $f\in
L^{2}(\Omega_{e})$ represents a source term having a compact support in
$\Omega_{e}$.

To determine the solution $u_{e}$ of Problem (1) by means of an interior
domain method, we surround the physical obstacle $\Omega_{0}$ by an artificial
boundary $\Gamma$; this 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}$ (see Figure 1 (b)). For the theoretical analysis of the
numerical approach we propose, we need to assume that $\Gamma_{0}$ consists of
a finite number of curves of class $C^{m+1}$, with $m\geq 0$, and that
$\Gamma$ is a contour of class $C^{\infty}$.

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}$ (consequently
$\mathbf{n}_{\infty}=-\mathbf{n}$), we consider the following compatibility
and equilibrium conditions on $\Gamma$:

$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.$
(4)

In the above relations and in the sequel we omit, for simplicity, the use of
the trace operators to indicate the restriction of $H^{1}$ functions to the
boundary $\Gamma$ from the exterior or interior.

Assuming, for simplicity, that $\Gamma$ is chosen such that $\text{supp}(f)$
is a bounded subset of $\Omega$, the following Kirchhoff’s formula

$u_{\infty}(\mathbf{x})=\int_{\Gamma}G(\mathbf{x},\mathbf{y})\frac{\partial
u_{\infty}}{\partial\mathbf{n}_{\infty}}(\mathbf{y})\,\differential\Gamma_{\mathbf{y}}-\int_{\Gamma}\frac{\partial
G}{\partial\mathbf{n}_{\infty,\mathbf{y}}}(\mathbf{x},\mathbf{y})u_{\infty}(\mathbf{y})\,\differential\Gamma_{\mathbf{y}}+\alpha,\quad\mathbf{x}\in\Omega_{\infty}\setminus\Gamma,$
(5)

allows us to represent the solution $u_{\infty}$ in $\Omega_{\infty}$. In (5),
$G$ and $\partial G/\partial\mathbf{n}_{\infty,\mathbf{y}}$ denote,
respectively, the fundamental solution of the 2D Laplace equation and its
normal derivative with respect to the unit vector
$\mathbf{n}_{\infty,\mathbf{y}}$ having initial point in
$\mathbf{y}\in\Gamma$. Their expression is given by

$G(\mathbf{x},\mathbf{y})=-\frac{1}{2\pi}\log
r\quad\text{and}\quad\frac{\partial
G}{\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}},$

where $r=\|\mathbf{r}\|=\|\mathbf{x}-\mathbf{y}\|$. It is known that the trace
of (5) on $\Gamma$ reads

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

where $\text{V}\colon H^{-\nicefrac{{1}}{{2}}}(\Gamma)\to
H^{\nicefrac{{1}}{{2}}}(\Gamma)$ and $\text{K}\colon
H^{\nicefrac{{1}}{{2}}}(\Gamma)\to H^{\nicefrac{{1}}{{2}}}(\Gamma)$ represent,
respectively, the continuous (see [23]) single- and double-layer integral
operators, defined by

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

and

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

To determine the solution of Problem (1) in the finite computational domain
$\Omega$, we impose (6) as BI-NRBC on $\Gamma$. In particular, introducing the
additional unknown $\displaystyle\lambda({\mathbf{y}}):=\frac{\partial
u}{\partial\mathbf{n}}(\mathbf{y})$ and taking into account (4), the new
problem defined in $\Omega$ takes the form:

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

We point out that the asymptotic conditions (2) coupled with (7c) imply that
$\langle\lambda,1\rangle_{\Gamma}=0$, where
$\langle\cdot,\cdot\rangle_{\Gamma}$ denotes the duality pairing between
$H^{-\nicefrac{{1}}{{2}}}(\Gamma)$ and $H^{\nicefrac{{1}}{{2}}}(\Gamma)$. This
justifies the introduction of the space
$H_{0}^{-\nicefrac{{1}}{{2}}}(\Gamma):=\left\\{\lambda\in
H^{-\nicefrac{{1}}{{2}}}(\Gamma)\ :\
\langle\lambda,1\rangle_{\Gamma}=0\right\\}$ in which we will look for the
unknown $\lambda$.

##  3 The variational formulation

Let us introduce the bilinear form $a:H^{1}(\Omega)\times
H^{1}(\Omega)\rightarrow\mathbf{R}$

$a(u,v):=\int\limits_{\Omega}\nabla u(\mathbf{x})\cdot\nabla
v(\mathbf{x})\,\differential\mathbf{x}.$

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

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

where I stands for the identity operator and $(\cdot,\cdot)_{L^{2}(\Omega)}$
denotes the $L^{2}(\Omega)$-inner product. It is worth noting that, since we
test Equation (7c) with $\mu\in H_{0}^{-\nicefrac{{1}}{{2}}}(\Gamma)$,
satisfying by definition $\langle\mu,1\rangle_{\Gamma}=0$, the unknown
constant $\alpha$ does not appear in the variational formulation (8).
Nevertheless, the asymptotic behaviour $\alpha$ is intrinsic to the interior
domain problem, and it can be recovered by the numerical scheme when choosing
$\Gamma$ sufficiently far from the obstacle (see Example 7.2).

To reformulate the above problem in operator form, following [24], we
introduce the Hilbert space $V:=H^{1}_{0,\Gamma_{0}}(\Omega)\times
H_{0}^{-\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).$

Then, we define the bilinear form $\mathcal{A}:V\times V\rightarrow\mathbf{R}$

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

for $\hat{u}=(u,\lambda)$ and $\hat{v}=(v,\mu)$, and the linear continuous
operator $\mathcal{L}_{f}:V\rightarrow\mathbf{R}$

$\mathcal{L}_{f}(\hat{v}):=(f,v)_{L^{2}(\Omega)}.$

Hence, we rewrite Problem (8) as follows: find $\hat{u}\in V$ such that

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

whose well-posedness has been proved in [24] (see Lemma 2).

Finally, for the forthcoming analysis, it is convenient to rewrite
$\mathcal{A}=\mathcal{B}+\mathcal{K}$ where the the bilinear forms
$\mathcal{B},\mathcal{K}:V\times V\rightarrow\mathbf{R}$ are defined as
follows:

$\mathcal{B}(\hat{u},\hat{v}):=a(u,v)-\langle\lambda,v\rangle_{\Gamma}+\langle\mu,u\rangle_{\Gamma}+2\langle\mu,\text{V}\lambda\rangle_{\Gamma},\,\,\,\mathcal{K}(\hat{u},\hat{v}):=-2\langle\mu,\text{K}u\rangle_{\Gamma}.$
(10)

In the following sections, for the solution of Problem (9), we will describe a
numerical approach consisting of a CVEM-BEM coupling. This method and the
corresponding theoretical analysis is based on that proposed for the Helmholtz
problem in [16], to which we refer whenever the theoretical results therein
proved hold in our context as well. It is worth noting that the theoretical
analysis for the exterior Poisson problem cannot be obtained as a particular
sub-case of the Helmholtz one given in [16], by simply choosing the wave
number equal to zero. Indeed, the NRBC associated to the Laplace equation is
different from that of the Helmholtz one, both for what concerns the kernel
functions appearing in the boundary integral operators and the choice of the
discrete function spaces for the approximation of the unknown $\lambda$. In
fact, in this case, since we do not know a priori the asymptotic value
$\alpha$ in (7c), the choice of the space
$H_{0}^{-\nicefrac{{1}}{{2}}}(\Gamma)$ becomes mandatory and, as a
consequence, a proper discrete subspace of it is needed. Moreover, another
important novelty of the theoretical study, with respect to that of [16],
consists in the use of decoupled degrees of approximation for the interior
CVEM and the BEM. This allows in particular the application of the CVEM with
order higher than that of the BEM, a key aspect for the global scheme since
the BEM requires high efforts to efficiently compute the associated system
matrices.

##  4 The numerical method

To describe the numerical approach we propose to solve (9), we start by
introducing a suitable decomposition of the domain $\Omega$, which consists of
generic elements and is not limited to the more commonly used triangles.

Let us denote by $E$ a generic “polygon” having at most one curved edge and by
$h_{E}$ its diameter; similarly we denote by $e$ a generic “edge”, eventually
curved, and by $h_{e}$ its length. We introduce a sequence
$\\{\mathcal{T}_{h_{\circ}}\\}_{h_{\circ}}$ of unstructured meshes
$\mathcal{T}_{h_{\circ}}=\left\\{E\right\\}$, which cover the domain $\Omega$,
where $h_{\circ}:=\max_{E\in\mathcal{T}_{h_{\circ}}}h_{E}$. We denote by
$\mathcal{T}_{h_{\partial}}^{\Gamma}$ the decomposition of the artificial
boundary $\Gamma$ which, according to the regularity assumption required for
$\Gamma$, consists of curvilinear parts. The subscript $h_{\partial}$ denotes
the mesh size defined by
$h_{\partial}=\max_{e\in\mathcal{T}^{\Gamma}_{h_{\partial}}}h_{e}$.

We suppose there exists a constant $\varrho>0$ such that for each $h_{\circ}$
and for each element $E\in\mathcal{T}_{h_{\circ}}$, $E$ is star-shaped with
respect to a ball of radius greater than $\varrho h_{E}$ and the length of any
(eventually curved) edge of $E$ is greater than $\varrho h_{E}$.

For any $k\in\mathbf{N}$, let $P_{k}(E)$ be the space of polynomials of degree
$k$ defined on $E$, and $\Pi_{k}^{\nabla,E}:H^{1}(E)\rightarrow P_{k}(E)$ be
the local polynomial $H^{1}$-projection, defined such that for $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}$

The local projector $\Pi_{k}^{\nabla,E}$ can be naturally extended to the
global one $\Pi_{k}^{\nabla}:H^{1}(\Omega)\rightarrow
P_{k}(\mathcal{T}_{h_{\circ}})$ as follows:

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

$P_{k}(\mathcal{T}_{h_{\circ}})$ being the space of piecewise polynomials with
respect to the decomposition $\mathcal{T}_{h_{\circ}}$ of $\Omega$. Moreover,
let $\Pi_{k}^{0,E}:L^{2}(E)\rightarrow P_{k}(E)$ be the local polynomial
$L^{2}$-projection operator, defined such that for $v\in L^{2}(E)$

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

By introducing the local bilinear form $a^{\text{{E}}}:H^{1}(E)\times
H^{1}(E)\rightarrow\mathbf{R}$ given by

$\displaystyle a^{\text{{E}}}(u,v):=\int\limits_{E}\nabla
u(\mathbf{x})\cdot\nabla v(\mathbf{x})\,\differential\mathbf{x},$ (11)

we can write
$a(u,v)=\sum\limits_{E\in\mathcal{T}_{\bm{h}}}a^{\text{{E}}}(u,v)$.

Finally, we introduce the product space
$H^{1}(\mathcal{T}_{h_{\circ}}):=\prod_{E\in\mathcal{T}_{h_{\circ}}}H^{1}(E)$
and define the associated broken $H^{1}$-norm:

$\|v\|_{H^{1}(\mathcal{T}_{h_{\circ}})}^{2}:=\sum_{E\in\mathcal{T}_{h_{\circ}}}\|v\|^{2}_{H^{1}(E)}.$

To apply the Galerkin method to Problem (9), we introduce the discrete spaces
$Q_{h_{\circ}}^{k_{\circ}}\subset H_{0,\Gamma_{0}}^{1}(\Omega)$ and
$X_{h_{\partial}}^{k_{\partial}}\subset H_{0}^{-\nicefrac{{1}}{{2}}}(\Gamma)$
associated to the meshes $\mathcal{T}_{h_{\circ}}$ and
$\mathcal{T}_{h_{\partial}}^{\Gamma}$, respectively, and the product space
$V_{\bm{h}}^{\bm{k}}:=Q_{h_{\circ}}^{k_{\circ}}\times
X_{h_{\partial}}^{k_{\partial}}$. Then, the Galerkin method consists in
finding $\hat{u}_{\bm{h}}\in V_{\bm{h}}^{\bm{k}}$ such that

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

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

Proceeding analogously as in [16], we introduce sufficient conditions on the
discrete spaces, on the bilinear form $\mathcal{B}_{\bm{h}}$ and on the linear
operator $\mathcal{L}_{f,\bm{h}}$ to guarantee existence and uniqueness of the
solution $\hat{u}_{\bm{h}}\in V_{\bm{h}}^{\bm{k}}$ and to prove convergence
error estimates.

In particular we assume: for any $s\geq 1$

  1. (H1.a)

approximation in $Q_{h_{\circ}}^{k_{\circ}}$: for all $v\in H^{s+1}(\Omega)$
$\underset{v_{h_{\circ}}\in
Q_{h_{\circ}}^{k_{\circ}}}{\text{inf}}\left\|v-v_{h_{\circ}}\right\|_{H^{1}(\Omega)}\apprle
h_{\circ}^{\text{min}(s,k_{\circ})}\left\|v\right\|_{H^{s+1}(\Omega)}$;

  2. (H1.b)

approximation in $X_{h_{\partial}}^{k_{\partial}}$: for all $\mu\in
H^{s-\nicefrac{{1}}{{2}}}(\Gamma)\cap H_{0}^{-\nicefrac{{1}}{{2}}}(\Gamma)$
$\underset{\mu_{h_{\partial}}\in
X_{h_{\partial}}^{k_{\partial}}}{\text{inf}}\left\|\mu-\mu_{h_{\partial}}\right\|_{H^{-\nicefrac{{1}}{{2}}}(\Gamma)}\apprle
h_{\partial}^{\text{min}(s,k_{\partial})}\left\|\mu\right\|_{H^{s-\nicefrac{{1}}{{2}}}(\Gamma)}.$

In the above assumptions the notation $Q_{1}\apprle Q_{2}$ (as well as
$Q_{1}\apprge Q_{2}$ in what follows) 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}$.

According to the definition of the $\|\cdot\|_{V}$ norm, (H1.a) and (H1.b)
ensure the following approximation property for the product space
$V_{\bm{h}}^{\bm{k}}$:

for $s\geq 1$, given $\hat{v}=(v,\mu)\in
H^{s+1}(\Omega)\times(H^{s-\nicefrac{{1}}{{2}}}(\Gamma)\cap
H_{0}^{-\nicefrac{{1}}{{2}}}(\Gamma))$, there exists
$\hat{v}_{\bm{h}}^{I}=(v_{h_{\circ}}^{I},\mu_{h_{\partial}}^{I})\in
V_{\bm{h}}^{\bm{k}}$ such that

$\left\|\hat{v}-\hat{v}_{\bm{h}}^{I}\right\|_{V}\apprle
h_{\circ}^{\text{min}(s,k_{\circ})}\left\|v\right\|_{H^{s+1}(\Omega)}+h_{\partial}^{\text{min}(s,k_{\partial})}\left\|\mu\right\|_{H^{s-\nicefrac{{1}}{{2}}}(\Gamma)}.$
(13)

Recalling that the evaluation of the bilinear form $\mathcal{B}$ on elements
of $V_{\bm{h}}^{\bm{k}}$ is well defined provided that $a(\cdot,\cdot)$ is
split into the sum of the local contributions $a^{\text{{E}}}(\cdot,\cdot)$,
and assuming that the approximated bilinear form $\mathcal{B}_{\bm{h}}$ is
well defined on the space $H^{1}(\mathcal{T}_{h_{\circ}})$, we further assume:

  1. (H2.a)

$k_{\circ}$-consistency: for all $\hat{q}\in
P_{k_{\circ}}(\mathcal{T}_{h_{\circ}})\times X_{h_{\partial}}^{k_{\partial}}$
and $\hat{v}_{\bm{h}}\in V_{\bm{h}}^{\bm{k}}$

$\mathcal{B}_{\bm{h}}(\hat{q},\hat{v}_{\bm{h}})=\mathcal{B}(\hat{q},\hat{v}_{\bm{h}}),\qquad\mathcal{B}_{\bm{h}}(\hat{v}_{\bm{h}},\hat{q})=\mathcal{B}(\hat{v}_{\bm{h}},\hat{q});$

  2. (H2.b)

continuity: for all $\hat{v}_{\bm{h}},\hat{w}_{\bm{h}}\in V_{\bm{h}}^{\bm{k}}$

$\absolutevalue*{\mathcal{B}_{\bm{h}}(\hat{v}_{\bm{h}},\hat{w}_{\bm{h}})}\apprle\left\|\hat{v}_{\bm{h}}\right\|_{V}\left\|\hat{w}_{\bm{h}}\right\|_{V};$

  3. (H2.c)

ellipticity: for all $\hat{w}_{\bm{h}}\in V_{\bm{h}}^{\bm{k}}$

$\mathcal{B}_{\bm{h}}(\hat{w}_{\bm{h}},\hat{w}_{\bm{h}})\apprge\left\|\hat{w}_{\bm{h}}\right\|_{V}^{2}.$

In the following theorem we show the validity of the inf-sup condition for the
discrete bilinear form $\mathcal{A}_{\bm{h}}$.

######  Theorem 4.1.

Assuming (H1.a), (H1.b) and (H2.a)–(H2.c), for $h_{\circ}$ and $h_{\partial}$
small enough, it holds that

$\sup_{\hat{v}_{\bm{h}}\in V_{\bm{h}}^{\bm{k}},\hat{v}_{\bm{h}}\neq
0}\frac{\mathcal{A}_{\bm{h}}(\hat{w}_{\bm{h}},\hat{v}_{\bm{h}})}{\lVert\hat{v}_{\bm{h}}\rVert_{V}}\apprge\lVert\hat{w}_{\bm{h}}\rVert_{V}\quad\forall\,\hat{w}_{\bm{h}}\in
V_{\bm{h}}^{\bm{k}}.$

###### Proof.

Following the proof of Lemma 4 in [24], it is possible to assert that, for any
$\hat{w}_{\bm{h}}\in V_{\bm{h}}^{\bm{k}}$, there exists $\hat{v}_{\bm{h}}\in
V_{\bm{h}}^{\bm{k}}$ such that

$\lVert\hat{v}_{\bm{h}}\rVert_{V}\apprle\lVert\hat{w}_{\bm{h}}\rVert_{V}$ (14)

and

$\mathcal{A}(\hat{w}_{\bm{h}},\hat{v}_{\bm{h}})\apprge(1-h_{\circ}-h_{\partial})\lVert\hat{w}_{\bm{h}}\rVert^{2}_{V}.$
(15)

By exploiting the decoupled assumptions (H1.a) and (H1.b) in Lemma 4.5 of
[16], we obtain that for
$\hat{v}_{\bm{h}}=(v_{h_{\circ}},\mu_{h_{\partial}})\in
V_{\bm{h}}^{\bm{k}}\subset V$ there exists a unique
$\hat{v}_{\bm{h}}^{*}=(v_{h_{\circ}}^{*},\mu_{h_{\partial}}^{*})\in
V_{\bm{h}}^{\bm{k}}$ such that

$\mathcal{B}_{\bm{h}}(\hat{w}_{\bm{h}},\hat{v}_{\bm{h}}^{*})=\mathcal{B}(\hat{w}_{\bm{h}},\hat{v}_{\bm{h}})\quad\forall\,\hat{w}_{\bm{h}}\in
V_{\bm{h}}^{\bm{k}}$ (16)

and

$\lVert\mu_{h_{\partial}}^{*}-\mu_{h_{\partial}}\rVert_{H^{-\nicefrac{{3}}{{2}}}(\Gamma)}\apprle(h_{\circ}+h_{\partial})\lVert\hat{v}_{\bm{h}}\rVert_{V}.$
(17)

Recalling (10) and (12), and using (16), we get:

$\displaystyle\mathcal{A}_{\bm{h}}(\hat{w}_{\bm{h}},\hat{v}_{\bm{h}}^{*})$
$\displaystyle=\mathcal{B}_{\bm{h}}(\hat{w}_{\bm{h}},\hat{v}_{\bm{h}}^{*})+\mathcal{K}(\hat{w}_{\bm{h}},\hat{v}_{\bm{h}}^{*})=\mathcal{B}(\hat{w}_{\bm{h}},\hat{v}_{\bm{h}})+\mathcal{K}(\hat{w}_{\bm{h}},\hat{v}_{\bm{h}}^{*})$
$\displaystyle=\mathcal{A}(\hat{w}_{\bm{h}},\hat{v}_{\bm{h}})+\mathcal{K}(\hat{w}_{\bm{h}},\hat{v}_{\bm{h}}^{*})-\mathcal{K}(\hat{w}_{\bm{h}},\hat{v}_{\bm{h}})$
$\displaystyle=\mathcal{A}(\hat{w}_{\bm{h}},\hat{v}_{\bm{h}})-2\langle\mu_{h_{\partial}}^{*}-\mu_{h_{\partial}},\text{K}w_{\bm{h}}\rangle_{\Gamma}.$

By applying the Hölder’s inequality and using (17) we have

$\displaystyle\absolutevalue*{\langle\mu_{h_{\partial}}^{*}-\mu_{h_{\partial}},\text{K}w_{h_{\circ}}\rangle_{H^{-\nicefrac{{3}}{{2}}}(\Gamma)\times
H^{\nicefrac{{3}}{{2}}}(\Gamma)}}$
$\displaystyle\apprle\lVert\mu_{h_{\partial}}^{*}-\mu_{h_{\partial}}\rVert_{H^{-\nicefrac{{3}}{{2}}}(\Gamma)}\lVert\text{K}w_{h_{\circ}}\rVert_{H^{\nicefrac{{3}}{{2}}}(\Gamma)}$
$\displaystyle\apprle(h_{\circ}+h_{\partial})\lVert\hat{v}_{\bm{h}}\rVert_{V}\lVert\text{K}w_{h_{\circ}}\rVert_{H^{\nicefrac{{3}}{{2}}}(\Gamma)}.$
(18)

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

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

which, together with (18), implies

$\absolutevalue*{\langle\mu_{h_{\partial}}^{*}-\mu_{h_{\partial}},\text{K}w_{h_{\circ}}\rangle_{H^{-\nicefrac{{3}}{{2}}}(\Gamma)\times
H^{\nicefrac{{3}}{{2}}}(\Gamma)}}\apprle(h_{\circ}+h_{\partial})\lVert\hat{v}_{\bm{h}}\rVert_{V}\lVert\hat{w}_{\bm{h}}\rVert_{V}.$
(19)

Finally, combining (15), (19) and (14), we get

$\displaystyle\mathcal{A}_{\bm{h}}(\hat{w}_{\bm{h}},\hat{v}_{\bm{h}}^{*})$
$\displaystyle\apprge(1-h_{\circ}-h_{\partial})\lVert\hat{w}_{\bm{h}}\rVert_{V}^{2}-(h_{\circ}+h_{\partial})\lVert\hat{w}_{\bm{h}}\rVert_{V}\lVert\hat{v}_{\bm{h}}\rVert_{V}$
$\displaystyle\apprge(1-h_{\circ}-h_{\partial})\lVert\hat{w}_{\bm{h}}\rVert_{V}^{2},$

whence, for $h_{\circ}$ and $h_{\partial}$ small enough, the claim follows. ∎

Theorem 4.1 allows us to prove the following convergence error estimate in the
$V$-norm for Problem (12).

######  Theorem 4.2.

Assume there exist $k_{\circ},k_{\partial}$ such that for all $s\geq 1$,
(H1.a),(H1.b), and (H2.a)–(H2.c) hold. Furthermore, assume that, for all
$s\geq 1$, there exists $\sigma_{s}:L^{2}(\Omega)\to\mathbf{R}^{+}$ such that

  1. (H3.a)

for all $\hat{v}_{\bm{h}}\in V_{\bm{h}}^{\bm{k}}$

$\absolutevalue*{\mathcal{L}_{f}(\hat{v}_{\bm{h}})-\mathcal{L}_{f,\bm{h}}(\hat{v}_{\bm{h}})}\apprle\left(h_{\circ}^{\text{min}(s,k_{\circ})}+h_{\partial}^{\text{min}(s,k_{\partial})}\right)\left\|\hat{v}_{\bm{h}}\right\|_{V}\,\sigma_{s}(f).$

Then, for $h_{\circ}$ and $h_{\partial}$ small enough, Problem (12) admits a
unique solution $\hat{u}_{\bm{h}}\in V_{\bm{h}}^{\bm{k}}$ and if
$\hat{u}=(u,\lambda)$, solution of Problem (9), satisfies $\hat{u}\in
H^{s+1}(\Omega)\times H^{s-\nicefrac{{1}}{{2}}}(\Gamma)$, the following
estimate holds

$\lVert\hat{u}-\hat{u}_{\bm{h}}\rVert_{V}\apprle\left(h_{\circ}^{\text{min}(s,k_{\circ})}+h_{\partial}^{\text{min}(s,k_{\partial})}\right)\left(\lVert
u\rVert_{H^{s+1}(\Omega)}+\sigma_{s}(f)\right).$

Since the proof of Theorem 4.2 can be obtained by proceeding analogously as in
Theorem 4.8 of [16], for brevity we omit it here.

In what follows we provide the error estimate in the weaker $W$-norm, where
$W:=L^{2}(\Omega)\times H^{-\nicefrac{{3}}{{2}}}(\Gamma)$. To this aim, we
start by assuming the following property:

  1. (H3.b)

consistency: for all $\hat{q}\in P_{1}(\mathcal{T}_{h_{\circ}})\times
X_{h_{\partial}}^{k_{\partial}}$

$\mathcal{L}_{f,\bm{h}}(\hat{q})=\mathcal{L}_{f}(\hat{q}).$

Such assumption, together with some of those previously introduced, allows us
to prove the following approximation error estimate for the operator
$\mathcal{L}_{f}$.

######  Lemma 4.3.

Let $\hat{v}=(v,\mu)\in H^{2}(\Omega)\times H^{-\nicefrac{{1}}{{2}}}(\Gamma)$,
and let $v_{h_{\circ}}^{I}$ be the best approximation of $v$ in
$Q_{h_{\circ}}^{k_{\circ}}$. Under assumptions (H1.a), (H3.a) and (H3.b), for
$s\geq 1$, it holds:

$\absolutevalue*{\mathcal{L}_{f}((v_{h_{\circ}}^{I},\mu))-{\mathcal{L}_{f,\bm{h}}}((v_{h_{\circ}}^{I},\mu))}\apprle
h_{\circ}\left(h_{\circ}^{\text{min}(s,k_{\circ})}+h_{\partial}^{\text{min}(s,k_{\partial})}\right)\lVert
v\rVert_{H^{2}(\Omega)}\sigma_{s}(f).$

###### Proof.

Using (H3.b) and (H3.a), we can write

$\displaystyle|\mathcal{L}_{f}((v_{h_{\circ}}^{I},\mu))-{\mathcal{L}_{f,\bm{h}}}((v_{h_{\circ}}^{I},\mu))|$
$\displaystyle=|\mathcal{L}_{f}((v_{h_{\circ}}^{I},\mu))-{\mathcal{L}_{f,\bm{h}}}((v_{h_{\circ}}^{I},\mu))-\mathcal{L}_{f}((\Pi_{1}^{\nabla}v,\mu))+\mathcal{L}_{f,\bm{h}}((\Pi_{1}^{\nabla}v,\mu))|$
$\displaystyle=|\mathcal{L}_{f}((v_{h_{\circ}}^{I}-\Pi_{1}^{\nabla}v,0))-{\mathcal{L}_{f,\bm{h}}}((v_{h_{\circ}}^{I}-\Pi_{1}^{\nabla}v,0))|$
$\displaystyle\apprle\left(h_{\circ}^{\text{min}(s,k_{\circ})}+h_{\partial}^{\text{min}(s,k_{\partial})}\right)\|v_{h_{\circ}}^{I}-\Pi^{\nabla}_{1}v\|_{H^{1}(\mathcal{T}_{h_{\circ}})}\,\sigma_{s}(f).$

From (H1.a) and by using standard polynomial approximation estimates (see, for
example, Lemma 4.3.8 in [13]), we get

$\displaystyle\|\Pi_{1}^{\nabla}v-v_{h_{\circ}}^{I}\|_{H^{1}(\mathcal{T}_{h_{\circ}})}\leq\|\Pi_{1}^{\nabla}v-v\|_{H^{1}(\mathcal{T}_{h_{\circ}})}+\|v-v_{h_{\circ}}^{I}\|_{H^{1}(\Omega)}\apprle
h_{\circ}\left\|v\right\|_{H^{2}(\Omega)},$ (20)

which, combined with the previous estimate, leads to the claim. ∎

To prove the error estimate in the weaker $W$-norm, we introduce the dual
space $W^{\prime}:=L^{2}(\Omega)\times H^{\nicefrac{{3}}{{2}}}(\Gamma)$ and
denote by $[\cdot,\cdot]:W\times W^{\prime}\to\mathbf{R}$ the associated
duality pairing. Further, we define the adjoint operator
$\mathcal{A}^{*}:V\rightarrow V^{{}^{\prime}}$ as

$\left(\mathcal{A}^{*}\hat{v}\right)(\hat{u}):=\mathcal{A}(\hat{u},\hat{v}),$

which, by Lemma 3 in [24], is an isomorphism, whose inverse
$\mathcal{A}^{*-1}:H^{1}(\Omega)\times
H^{\nicefrac{{3}}{{2}}}(\Gamma)\rightarrow H^{2}(\Omega)\times
H^{\nicefrac{{1}}{{2}}}(\Gamma)$ is continuous.

######  Theorem 4.4.

Assume that there exist $k_{\circ}$ and $k_{\partial}$ and
$\sigma_{s}:L^{2}(\Omega)\to\mathbf{R}^{+}$ such that, for all $s\geq 1$,
(H1.a),(H1.b), (H2.a)–(H2.c), (H3.a) and (H3.b) hold. Then, for $h_{\circ}$
and $h_{\partial}$ small enough, if $\hat{u}_{\bm{h}}\in V_{\bm{h}}^{\bm{k}}$
is the solution of Problem (12) and $\hat{u}$, solution of Problem (9),
satisfies $\hat{u}\in H^{s+1}(\Omega)\times
H^{s-\nicefrac{{1}}{{2}}}(\Gamma)$, for $s\geq 1$, the following estimate

$\lVert\hat{u}-\hat{u}_{\bm{h}}\rVert_{W}\apprle(h_{\circ}+h_{\partial})\left(h_{\circ}^{\text{min}(s,k_{\circ})}+h_{\partial}^{\text{min}(s,k_{\partial})}\right)\left(\lVert
u\rVert_{H^{s+1}(\Omega)}+\sigma_{s}(f)\right)$

holds.

###### Proof.

Let $\hat{w}\in W^{\prime}$ and $\hat{v}:=\mathcal{A}^{*-1}\hat{w}$ the unique
element that, according to the above mentioned property of
$\mathcal{A}^{*-1}$, satisfies $\hat{v}\in H^{2}(\Omega)\times
H^{\nicefrac{{1}}{{2}}}(\Gamma)\subset V$,

$\mathcal{A}(\hat{z},\hat{v})=\mathcal{A}(\mathcal{A}^{*-1}\hat{w},\hat{z})=[\hat{z},\hat{w}]\quad\forall\hat{z}\in
V$ (21)

and $\lVert\hat{v}\rVert_{H^{2}(\Omega)\times
H^{\nicefrac{{1}}{{2}}}(\Gamma)}\apprle\lVert\hat{w}\rVert_{W^{\prime}}$.

Now, choosing $\hat{z}=\hat{u}-\hat{u}_{\bm{h}}$ in (21), where $\hat{u}$ and
$\hat{u}_{\bm{h}}$ are the solutions of (9) and (12) respectively, and
denoting by
$\hat{v}_{\bm{h}}^{I}=(v^{I}_{h_{\circ}},\mu^{I}_{h_{\partial}})\in
V_{\bm{h}}^{\bm{k}}$ the interpolant of $\hat{v}$, we obtain

$\displaystyle|\,[\hat{u}$
$\displaystyle-\hat{u}_{\bm{h}},\hat{w}]\,|=\absolutevalue*{\mathcal{A}(\hat{u}-\hat{u}_{\bm{h}},\hat{v})}$
$\displaystyle=\absolutevalue*{\mathcal{A}(\hat{u}-\hat{u}_{\bm{h}},\hat{v})+\mathcal{A}_{\bm{h}}(\hat{u}_{\bm{h}},\hat{v}^{I}_{\bm{h}})-\mathcal{A}_{\bm{h}}(\hat{u}_{\bm{h}},\hat{v}_{\bm{h}}^{I})+\mathcal{A}(\hat{u}-\hat{u}_{\bm{h}},\hat{v}_{\bm{h}}^{I})-\mathcal{A}(\hat{u}-\hat{u}_{\bm{h}},\hat{v}_{\bm{h}}^{I})}$
$\displaystyle\leq\absolutevalue*{\mathcal{A}(\hat{u}-\hat{u}_{\bm{h}},\hat{v}-\hat{v}_{\bm{h}}^{I})}+\absolutevalue*{\mathcal{L}_{f}(\hat{v}_{\bm{h}}^{I})-\mathcal{L}_{f,\bm{h}}(\hat{v}_{\bm{h}}^{I})}+\absolutevalue*{\mathcal{A}_{\bm{h}}(\hat{u}_{\bm{h}},\hat{v}_{\bm{h}}^{I})-\mathcal{A}(\hat{u}_{\bm{h}},\hat{v}_{\bm{h}}^{I})}$
$\displaystyle\apprle\lVert\hat{u}-\hat{u}_{\bm{h}}\rVert_{V}\lVert\hat{v}-\hat{v}_{\bm{h}}^{I}\rVert_{V}+h_{\circ}\left(h_{\circ}^{\text{min}(s,k_{\circ})}+h_{\partial}^{\text{min}(s,k_{\partial})}\right)\lVert
v\rVert_{H^{2}(\Omega)}\sigma_{s}(f)$
$\displaystyle+\absolutevalue*{\mathcal{B}_{\bm{h}}(\hat{u}_{\bm{h}},\hat{v}_{\bm{h}}^{I})-\mathcal{B}(\hat{u}_{\bm{h}},\hat{v}_{\bm{h}}^{I})}=:I+II+III,$
(22)

the last inequality following from the continuity of $\mathcal{A}$ and Lemma
4.3. By applying Theorem 4.2 and the interpolation property (13), we estimate

$\displaystyle I$
$\displaystyle\apprle(h_{\circ}+h_{\partial})\left(h_{\circ}^{\text{min}(s,k_{\circ})}+h_{\partial}^{\text{min}(s,k_{\partial})}\right)\left(\lVert
u\rVert_{H^{s+1}(\Omega)}+\sigma_{s}(f)\right)\|\hat{v}\|_{H^{2}(\Omega)\times
H^{\nicefrac{{1}}{{2}}}(\Gamma)}$
$\displaystyle\apprle(h_{\circ}+h_{\partial})\left(h_{\circ}^{\text{min}(s,k_{\circ})}+h_{\partial}^{\text{min}(s,k_{\partial})}\right)\left(\lVert
u\rVert_{H^{s+1}(\Omega)}+\sigma_{s}(f)\right)\|\hat{w}\|_{W^{\prime}},$ (23)

and, since $\lVert v\rVert_{H^{2}(\Omega)}\leq\|\hat{v}\|_{H^{2}(\Omega)\times
H^{\nicefrac{{1}}{{2}}}(\Gamma)}$ we have

$II\apprle
h_{\circ}\left(h_{\circ}^{\text{min}(s,k_{\circ})}+h_{\partial}^{\text{min}(s,k_{\partial})}\right)\|\hat{w}\|_{W^{\prime}}\sigma_{s}(f).$
(24)

To estimate $III$ in (22), we add and subtract the terms
$\mathcal{B}_{\bm{h}}((\Pi^{\nabla}_{k_{\circ}}u,\lambda),\hat{v}_{\bm{h}}^{I})$
and $\mathcal{B}((\Pi^{\nabla}_{k_{\circ}}u,\lambda),\hat{v}_{\bm{h}}^{I})$
which, for (H2.a), are equal. Hence we get

$\displaystyle III$
$\displaystyle=|\mathcal{B}_{\bm{h}}((u_{h_{\circ}}-\Pi^{\nabla}_{k_{\circ}}u,0),\hat{v}_{\bm{h}}^{I})-\mathcal{B}((u_{h_{\circ}}-\Pi^{\nabla}_{k_{\circ}}u,0),\hat{v}_{\bm{h}}^{I})|.$

Similarly, adding and subtracting the two equal terms (see (H3.b))
$\mathcal{B}_{\bm{h}}((u_{h_{\circ}}-\Pi_{k_{\circ}}^{\nabla}u,0),(\Pi^{\nabla}_{k_{\circ}}v,\mu_{h_{\partial}}^{I}))$
and
$\mathcal{B}((u_{h_{\circ}}-\Pi_{k_{\circ}}^{\nabla}u,0),(\Pi^{\nabla}_{k_{\circ}}v,\mu_{h_{\partial}}^{I}))$,
we obtain

$\displaystyle
III\leq|\mathcal{B}_{\bm{h}}((u_{h_{\circ}}-\Pi^{\nabla}_{k_{\circ}}u,0),(v_{h_{\circ}}^{I}-\Pi^{\nabla}_{1}v,0))|+|\mathcal{B}((u_{h_{\circ}}-\Pi^{\nabla}_{k_{\circ}}u,0),(v_{h_{\circ}}^{I}-\Pi^{\nabla}_{1}v,0))|.$

Using the continuity of $\mathcal{B}_{\bm{h}}$ (see (H2.b)) and of
$\mathcal{B}$, we have

$III\apprle\lVert
u_{h_{\circ}}-\Pi^{\nabla}_{k_{\circ}}u\rVert_{H^{1}(\mathcal{T}_{\bm{h}})}\lVert
v_{h_{\circ}}^{I}-\Pi^{\nabla}_{1}v\rVert_{H^{1}(\mathcal{T}_{\bm{h}})}.$

The first factor of the above inequality is estimated, by using Theorem 4.2
and standard polynomial approximations, as follows

$\displaystyle\|u_{h_{\circ}}-\Pi^{\nabla}_{k_{\circ}}u\|_{H^{1}(\mathcal{T}_{\bm{h}})}$
$\displaystyle\leq\|u_{h_{\circ}}-u\|_{H^{1}(\Omega)}+\lVert
u-\Pi^{\nabla}_{k_{\circ}}u\rVert_{H^{1}(\mathcal{T}_{\bm{h}})}$
$\displaystyle\apprle\left(h_{\circ}^{\text{min}(s,k_{\circ})}+h_{\partial}^{\text{min}(s,k_{\partial})}\right)(\lVert
u\rVert_{H^{s+1}(\Omega)}+\sigma_{s}(f)).$

Then, using (20), we obtain

$\displaystyle III$ $\displaystyle\apprle
h_{\circ}\left(h_{\circ}^{\text{min}(s,k_{\circ})}+h_{\partial}^{\text{min}(s,k_{\partial})}\right)\left(\lVert
u\rVert_{H^{s+1}(\Omega)}+\sigma_{s}(f)\right)\|v\|_{H^{2}(\Omega)}$
$\displaystyle\apprle
h_{\circ}\left(h_{\circ}^{\text{min}(s,k_{\circ})}+h_{\partial}^{\text{min}(s,k_{\partial})}\right)\left(\lVert
u\rVert_{H^{s+1}(\Omega)}+\sigma_{s}(f)\right)\|\hat{w}\|_{W^{\prime}}.$ (25)

Finally, the assertion easily follows combining (22) with (4), (24) and (4). ∎

##  5 The CVEM-BEM method

In this section we describe the discrete CVEM-BEM coupling procedure for the
solution of Problem (9). In particular, we will show that all the assumptions,
introduced in Section 4 and used to prove Theorems 4.2 and 4.4, are satisfied.
Referring to the notation introduced in Section 4, and denoting
$P_{-1}(E)=\\{0\\}$, we consider for each $E\in\mathcal{T}_{h_{\circ}}$ the
following local virtual space $Q^{k_{\circ}}_{h_{\circ}}(E)$ defined by

$\displaystyle Q^{k_{\circ}}_{h_{\circ}}(E):=$
$\displaystyle\left\\{v_{h_{\circ}}\in H^{1}(E):\Delta v_{h_{\circ}}\in
P_{k_{\circ}-2}(E),\right.$ $\displaystyle\left.\hskip
11.38092ptv_{h_{\circ}}\,\raisebox{-5.0pt}{$|_{e_{1}}$}\in\widetilde{P}_{k_{\circ}}(e_{1}),v_{h_{\circ}}\,\raisebox{-5.0pt}{$|_{e_{i}}$}\in
P_{k_{\circ}}(e_{i}),\ i=2,\ldots,n_{\text{{E}}}\right\\},$

where $e_{1},\ldots,e_{n_{\text{{E}}}}$ denote the edges of the boundary of
$E$, whose first element $e_{1}$ is assumed to be curved and parametrized by a
local map $\gamma_{E}:I_{E}\rightarrow e_{1}$, and
$\widetilde{P}_{k_{\circ}}(e_{1}):=\left\\{\widetilde{q}=q\circ\gamma_{E}^{-1}\
:\ q\in P_{k_{\circ}}(I_{E})\right\\}$.

We omit here, for brevity, the complete description of such space and we refer
to [3, 10] for a deeper presentation. Further, since we will use some of the
theoretical results proved in [16], we also refer to this latter, in
particular for what concerns the notation and the implementation details.

On the basis of the definition of the local virtual space
$Q_{h_{\circ}}^{k_{\circ}}(E)$, we construct the global virtual space

$Q_{h_{\circ}}^{k_{\circ}}:=\\{v_{h_{\circ}}\in
H^{1}_{0,\Gamma_{0}}:v_{{h_{\circ}}_{|_{E}}}\in Q_{h_{\circ}}^{k_{\circ}}(E),\
E\in\mathcal{T}_{h_{\circ}}\\},$

The validity of Assumption (H1.a) is based on the proof of Lemma 5.2 in [16],
in which the results hold both for the space $Q^{k_{\circ}}_{h_{\circ}}$ and
for a suitable _enhanced_ space associated to it. For each
$E\in\mathcal{T}_{h_{\circ}}$, in the spirit of the virtual element method, we
define the approximation $a^{\text{{E}}}_{h_{\circ}}$ of the bilinear form
$a^{\text{{E}}}$ (see definition (11)), as follows:

$\displaystyle a^{\text{{E}}}_{h_{\circ}}(u_{h_{\circ}},v_{h_{\circ}})$
$\displaystyle:=a^{\text{{E}}}\left(\Pi_{k_{\circ}}^{\nabla,E}u_{h_{\circ}},\Pi_{k_{\circ}}^{\nabla,E}v_{h_{\circ}}\right)+s^{\text{{E}}}\left(\left(I-\Pi_{k_{\circ}}^{\nabla,E}\right)u_{h_{\circ}},\left(I-\Pi_{k_{\circ}}^{\nabla,E}\right)v_{h_{\circ}}\right)$

where $s^{E}$ is the standard “dofi-dofi” stabilization term (see Eq. (3.22)
of [4]). The global approximate bilinear form
$a_{h_{\circ}}:Q^{k_{\circ}}_{h_{\circ}}\times
Q^{k_{\circ}}_{h_{\circ}}\rightarrow\mathbf{R}$ is then defined by summing up
the local contributions

$a_{h_{\circ}}(u_{h_{\circ}},v_{h_{\circ}}):=\sum\limits_{E\in\mathcal{T}_{\bm{h}}}a^{\text{{E}}}_{h_{\circ}}(u_{h_{\circ}},v_{h_{\circ}}).$

The boundary element space $X_{h_{\partial}}^{k_{\partial}}$ associated to the
artificial boundary $\Gamma$ is defined as follows:

$X_{h_{\partial}}^{k_{\partial}}:=\left\\{\lambda\in L^{2}(\Gamma)\
:\lambda_{\mkern 1.0mu\vrule height=6.02777pt\mkern
2.0mue}\in\widetilde{P}_{k_{\partial}-1}(e),\
e\in\mathcal{T}_{h_{\partial}}^{\Gamma}\right\\}\cap\,H_{0}^{-\nicefrac{{1}}{{2}}}(\Gamma),$

and we refer to [25] for the validity of the associated Assumption (H1.b). We
then define the approximate bilinear form
$\mathcal{B}_{\bm{h}}:V^{\bm{k}}_{\bm{h}}\times
V^{\bm{k}}_{\bm{h}}\rightarrow\mathbf{R}$ as:

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

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

For these choices, from [16] (see Section 5.2), Assumptions (H2.a)-(H2.c) are
satisfied .

By approximating the linear operator $\mathcal{L}_{f}$ in a standard VEM way,
in particular as in [12] (see Eq. (3.18)), and assuming $f\in
H^{s-1}(\Omega)$, from Lemma 3.4 in [12], Assumption (H3.a) follows with
$\sigma_{s}(f)=|f|_{H^{s-1}(\Omega)}$. Finally, Assumption (H3.b) is trivially
satisfied.

##  6 Algebraic details and computational issues

In this section we briefly detail the construction of the final linear system
associated with the CVEM-BEM scheme, and we give some implementation issues
concerning the BEM matrices.

We start by re-ordering and splitting the complete index set $\mathcal{S}$ of
the basis functions $\left\\{\Phi_{j}\right\\}_{j\in\mathcal{S}}$ of
$Q^{k_{\circ}}_{h_{\circ}}$ as
$\mathcal{S}=\mathcal{S}^{I}\cup\mathcal{S}^{\Gamma}$, where $\mathcal{S}^{I}$
and $\mathcal{S}^{\Gamma}$ denote the sets of indices related to the internal
degrees of freedom and to those lying on $\Gamma$, respectively. Moreover, we
denote by $\left\\{\varphi_{j}\right\\}_{j\in\mathcal{G}}$ the basis functions
of $X^{k_{\partial}}_{h_{\partial}}$, $\mathcal{G}$ being the corresponding
index set. In order to write the linear system associated with the discrete
problem (12), we expand the unknown function
$\hat{u}_{\bm{h}}=(u_{h_{\circ}},\lambda_{h_{\partial}})\in
Q_{h_{\circ}}^{k_{\circ}}\times X^{k_{\partial}}_{h_{\partial}}$ as

$\displaystyle
u_{h_{\circ}}(\mathbf{x})=:\sum\limits_{j\in\mathcal{S}}u_{h_{\circ}}^{j}\Phi_{j}(\mathbf{x})\quad\text{with}\quad
u_{h_{\circ}}^{j}=\text{dof}_{j}(u_{h_{\circ}}),$ (26)
$\displaystyle\lambda_{h_{\partial}}(\mathbf{x})=:\sum\limits_{j\in\mathcal{G}}\lambda_{h_{\partial}}^{j}\varphi_{j}(\mathbf{x})\quad\text{with}\quad\lambda_{h_{\partial}}^{j}=\text{dof}_{j}(\lambda_{h_{\partial}}).$

Hence, using the basis functions of $Q^{k_{\circ}}_{h_{\circ}}$ to test the
discrete counterpart of equation (8a), we get for $i\in\mathcal{S}$

$\displaystyle\sum\limits_{j\in\mathcal{S}}u_{h_{\circ}}^{j}$
$\displaystyle\sum\limits_{E\in\mathcal{T}_{h_{\circ}}}a^{\text{{E}}}_{h_{\circ}}(\Phi_{j},\Phi_{i})-\sum\limits_{j\in\mathcal{G}}\lambda_{h_{\partial}}^{j}\langle\varphi_{j},\Phi_{i}\rangle_{\Gamma}=\mathcal{L}_{f,\bm{h}}((\Phi_{i},0)).$
(27)

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

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

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

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

For what concerns the discretization of the BI-NRBC, by inserting (26) in (8b)
and testing with the functions $\varphi_{i}$, $i\in\mathcal{G}$, we obtain

$\displaystyle\sum\limits_{j\in\mathcal{S}^{\Gamma}}$ $\displaystyle
u_{h_{\circ}}^{j}\left[\frac{1}{2}\int\limits_{\Gamma}\Phi_{j}(\mathbf{x})\varphi_{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)\varphi_{i}(\mathbf{x})\differential\Gamma_{\mathbf{x}}\right]$
(41)
$\displaystyle\,\,+\sum\limits_{j\in\mathcal{G}}\lambda_{h_{\partial}}^{j}\int\limits_{\Gamma}\left(\int\limits_{\Gamma}G(\mathbf{x},\mathbf{y})\varphi_{j}(\mathbf{y})\,\differential\Gamma_{\mathbf{y}}\right)\varphi_{i}(\mathbf{x})\differential\Gamma_{\mathbf{x}}=0,$

that in matrix form reads

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

where

$\displaystyle\mathbb{V}_{ij}$
$\displaystyle:=\int\limits_{\Gamma}\left(\int\limits_{\Gamma}G(\mathbf{x},\mathbf{y})\varphi_{j}(\mathbf{y})\differential\Gamma_{\mathbf{y}}\right)\varphi_{i}(\mathbf{x})\differential\Gamma_{\mathbf{x}},$
$\displaystyle\mathbb{K}_{ij}$
$\displaystyle:=\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)\varphi_{i}(\mathbf{x})\differential\Gamma_{\mathbf{x}}.$

By combining (40) with (42) we obtain the final linear system

$\displaystyle\left[\begin{array}[]{ccc}\mathbb{A}^{II}&\mathbb{A}^{I\Gamma}&\mathbb{O}\\\
&\\\ \mathbb{A}^{\Gamma I}&\mathbb{A}^{\Gamma\Gamma}&-\mathbb{Q}\\\ &\\\
\mathbb{O}&\frac{1}{2}\mathbb{Q}^{T}-\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].$ (58)

It is worth to point out that, since in the theoretical analysis we have
assumed that the boundary integral operators are not approximated, it is
crucial to compute the integrals defining the BEM entries of $\mathbb{V}$ and
$\mathbb{K}$ with a high accuracy. Hence, to retrieve their approximation
without affecting the overall accuracy of the coupled CVEM-BEM scheme,
suitable efficient quadrature formulas must be considered. In [17] we have
proposed and successfully applied a smoothing technique to weaken the $\log$\-
singularity of the single layer operator and to compute the corresponding
entries with high accuracy by using the Gauss-Legendre product quadrature rule
with few nodes. Such a strategy has been tuned for the standard nodal linear
and quadratic basis functions and could be, in principle, adopted also in this
context for the basis functions satisfying property (3). However, since the
strategy associated to the standard Lagrangian basis is a well-established
task, we take advantage of it by applying a computational trick. To describe
this latter, we introduce the space

$\widehat{X}_{h_{\partial}}^{k_{\partial}}:=\\{\lambda\in
L^{2}(\Gamma):{\lambda}_{|_{e}}\in\widetilde{P}_{k_{\partial}-1}(e),\,\,e\in\mathcal{T}_{h_{\partial}}^{\Gamma}\\}\subset
H^{-\nicefrac{{1}}{{2}}}(\Gamma),$

whose Lagrangian basis functions are denoted by $\widehat{\varphi}_{j}$.

Further, we denote by $\widehat{\mathbb{V}},\widehat{\mathbb{K}}$ and
$\widehat{\mathbb{Q}}$ the matrices associated with the choice of the space
$\widehat{X}_{h_{\partial}}^{k_{\partial}}$, which differ from
${\mathbb{V}},{\mathbb{K}}$ and ${\mathbb{Q}}$ by the presence of the
functions $\widehat{\varphi}_{i}$ instead of $\varphi_{i}$. In the forthcoming
Remark 6.1 we detail the quadrature adopted to efficiently compute
$\widehat{\mathbb{V}}$. Here we describe how to retrieve the matrices
$\mathbb{V},\mathbb{K}$ and $\mathbb{Q}$ from the corresponding
$\widehat{\mathbb{V}},\widehat{\mathbb{K}}$ and $\widehat{\mathbb{Q}}$. To
this aim it is sufficient to define the functions $\varphi_{i}$ as a suitable
linear combination of the standard ones $\widehat{\varphi}_{i}$ and, hence, to
combine accordingly the rows and/or the columns of
$\widehat{\mathbb{V}},\widehat{\mathbb{K}}$ and $\widehat{\mathbb{Q}}$. Such a
combination depends on the order $k_{\partial}$, on the shape of the
artificial boundary $\Gamma$ and on the associated mesh
$\mathcal{T}_{h_{\partial}}^{\Gamma}$. In particular, for $k_{\partial}=2$, we
define $\varphi_{i}:=c_{i}\widehat{\varphi}_{i}+\widehat{\varphi}_{i+1}$,
where $\widehat{\varphi}_{i}$ and $\widehat{\varphi}_{i+1}$ are two
consecutive piece-wise linear nodal basis functions. For $k_{\partial}=3$, we
distinguish the following two cases: a)
$\varphi_{2i}:=c_{2i}\widehat{\varphi}_{2i}+\widehat{\varphi}_{2i+1}$; b)
$\varphi_{2i+1}:=\widehat{\varphi}_{2i+1}+c_{2i+1}\widehat{\varphi}_{2i+2}$.
The coefficients $c_{\ell}$ are chosen such that the relation
$\int_{\Gamma}\varphi_{\ell}=0$ is satisfied and are retrieved by applying a
$\nu$-point Gauss-Legendre quadrature formula, with $\nu$ chosen such that the
integral over $\Gamma$ of the associated $\widehat{\varphi}$ functions is
accurately computed. It is worth noting that
$\text{dim}(X_{h_{\partial}}^{k_{\partial}})=\text{dim}(\widehat{X}_{h_{\partial}}^{k_{\partial}})-1$.

In Figures 2 and 3, we show the basis functions of the spaces
$\widehat{X}_{h_{\partial}}^{k_{\partial}}$ and
$X_{h_{\partial}}^{k_{\partial}}$, with ${k_{\partial}}=2,3$ respectively,
associated with the uniformly partitioned parametrization interval $[0,2\pi)$
of the particular choice of a circumference. For this choice, it is immediate
to get $c_{i}=-1$ for ${k_{\partial}}=2$, and $c_{2i}=c_{2i+1}=-2$ for
${k_{\partial}}=3$.

Figure 2: Basis of $\hat{X}_{h_{\partial}}^{2}$ and $X_{h_{\partial}}^{2}$ in
$[0,2\pi)$

Figure 3: $\hat{X}_{h_{\partial}}^{3}$ and $X_{h_{\partial}}^{3}$ in
$[0,2\pi)$

######  Remark 6.1.

We recall that the numerical integration difficulties in the computation of
the $\widehat{\mathbb{V}}$ entries spring from the $\log$-singularity of
$G(r)$ near the origin, the latter being the kernel of the single layer
operator V. To compute the corresponding integrals with high accuracy by few
nodes, we have used the very simple and efficient polynomial smoothing
technique proposed in [27], referred as the _q-smoothing_ technique. It is
worth noting that such technique is applied only when the distance $r$
approaches zero. This case corresponds to the matrix entries belonging to the
main diagonal and to the 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 [18] and Remark 3 in [17] 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 kernel function $\partial_{\bf n}{G}$, appearing in the double
layer operator K, do not require a particular quadrature strategy, since its
singularity $1/r$ is factored out by the same behaviour of
$\partial_{\mathbf{n}}r$ near the origin. Hence, for the computation of the
entries of the matrix $\widehat{\mathbb{K}}$, we have directly applied a
$9\times 8$-point Gauss-Legendre product quadrature rule. The quadrature
strategy described above guarantees the computation of all the mentioned
integrals with a full precision accuracy (16-digit double precision
arithmetic) for both $k_{\partial}=2$ and $k_{\partial}=3$.

##  7 Numerical results

In this section, we present some numerical test to validate the theoretical
results and to show the effectiveness of the proposed method.

For the generation of the partitioning $\mathcal{T}_{h_{\circ}}$ of the
computational domain $\Omega$, we have used the GMSH software to construct
unstructured conforming meshes consisting of quadrilaterals (see [21]). If a
polygon $E$ has a (straight) edge bordering with the interior boundary or with
the artificial one, we transform it into a curved boundary edge by means of a
suitable parametrization. We remark that, even if in principle it is possible
to fully decouple the interior and boundary meshes, we consider here for
simplicity the boundary mesh inherited by the interior one, for which it turns
out $h_{\partial}\leq h_{\circ}$. Furthermore, we point out that in all the
numerical test we have considered $k_{\partial}=2,3$; larger values than those
considered would require a tailored quadrature technique for the accurate
computation of the BEM matrices that we have not performed yet. Since this
usually is considered the bottleneck of the BEM, the use of decoupled
approximation orders allows us to exploit the flexibility of the CVEM to
retrieve high accuracy by increasing only the approximation order $k_{\circ}$.
In Example 7.1 we will investigate this aspect.

###  7.1 Example 1

Let us consider the unbounded region $\Omega_{\text{e}}$, external to the
unitary disk $\Omega_{0}=\\{\mathbf{x}=(x_{1},x_{2})^{\top}\in\mathbf{R}^{2}\
:\ x_{1}^{2}+x_{2}^{2}\leq 1\\}$. We consider Problem (1) with $f=0$ and
$g(\mathbf{x})=x_{1}+2$ prescribed on the boundary
$\Gamma_{0}=\partial\Omega_{0}$. In this case, the exact solution
$u(\mathbf{x})$ is known and its expression is given by

$u(\mathbf{x})=\frac{x_{1}}{x_{1}^{2}+x_{2}^{2}}+2,\quad\
\mathbf{x}\in\Omega_{\text{e}}.$

We choose as artificial boundary $\Gamma$ the circumference of radius 2, so
that the finite computational domain $\Omega$ is the annulus bounded
internally by $\Gamma_{0}$ and externally by $\Gamma$.

To develop a convergence analysis, we start by considering a 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. 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_{\circ}=k_{\partial}=2,3$ (see Figure 4 for
the meshes corresponding to level 0 (left plot) and level 2 (right plot)).

Figure 4: 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
---|---|---|---|---|---|---
$k_{\circ}=2$ | $368$ | $1,376$ | $5,312$ | $20,864$ | $82,688$ | $329,216$
$k_{\circ}=3$ | $792$ | $3,024$ | $11,808$ | $46,656$ | $185,472$ | $739,584$
Table 1: Number of the degrees of freedom associated to the CVEM space.

To test our numerical approach and to validate the theoretical analysis, the
order $k_{\circ}$ of the approximation spaces is chosen equal to 2 (quadratic)
and 3 (cubic), and $k_{\partial}=k_{\circ}$. Moreover, recalling that the
approximate solution $u_{h_{\circ}}$ is not known inside the polygons, as
suggested in [10] we compute the $H^{1}$-seminorm and $L^{2}$-norm relative
errors, and the corresponding EOC, by means of the following formulas:

  * •

$H^{1}$-seminorm
${\varepsilon^{\nabla,k_{\circ}}_{\text{lev}}:=\sqrt{\frac{\sum\limits_{E\in\mathcal{T}_{h_{\circ}}}\absolutevalue*{u-\Pi_{k_{\circ}}^{\nabla,E}u_{h_{\circ}}}^{2}_{H^{1}(E)}}{\absolutevalue*{u}^{2}_{H^{1}(\Omega)}}}}$;

  * •

$L^{2}$-norm
$\varepsilon^{0,k_{\circ}}_{\text{lev}}:=\sqrt{\frac{\sum\limits_{E\in\mathcal{T}_{h_{\circ}}}\left\|u-\Pi_{k_{\circ}}^{0,E}u_{h_{\circ}}\right\|^{2}_{L^{2}(E)}}{\left\|u\right\|^{2}_{L^{2}(\Omega)}}}$;

  * •

$\text{EOC}:=\log_{2}\left(\frac{\varepsilon^{*,k_{\circ}}_{\text{lev}+1}}{\varepsilon^{*,k_{\circ}}_{\text{lev}}}\right),\quad*\in\\{\nabla,0\\}$.

In the above formulas the superscript $k_{\circ}=2,3$ refers to the
approximation order of $u$, and the subscript lev refers to the refinement
level. For what concerns the evaluation of these errors, to compute the
associated integrals over polygons we have used the $n$-point quadrature
formulas proposed in [29] and [30], which are exact for polynomials of degree
at most $2n$. For curved polygons, we have applied the generalization of these
formulas suggested in [10] (see Section 4.3). In both cases, we have chosen
$n=8$.

In Table 2 we report $\varepsilon^{\nabla,k_{\circ}}_{\text{lev}}$ and
$\varepsilon^{0,k_{\circ}}_{\text{lev}}$ and the corresponding EOC by varying
the refinement level from 0 to 5. As we can see the $H^{1}$-seminorm error and
the $L^{2}$-norm error estimates confirm the expected convergence order of the
method.

|  $L^{2}$-norm |  $H^{1}$-seminorm
---|---|---
lev. | $h_{\circ}$ | $\varepsilon^{0,2}_{\text{lev}}$ | EOC | $\varepsilon^{0,3}_{\text{lev}}$ | EOC | $\varepsilon^{\nabla,2}_{\text{lev}}$ | EOC | $\varepsilon^{\nabla,3}_{\text{lev}}$ | EOC
$0$ | $8.02e-01$ | $4.26e-04$ |  | $6.74e-05$ |  | $4.96e-04$ |  | $1.05e-04$ |
|  |  | $2.9$ |  | $3.9$ |  | $1.9$ |  | $2.8$
$1$ | $4.28e-01$ | $5.56e-05$ |  | $4.58e-06$ |  | $1.36e-04$ |  | $1.51e-05$ |
|  |  | $3.0$ |  | $4.0$ |  | $2.0$ |  | $3.0$
$2$ | $2.22e-01$ | $7.05e-06$ |  | $2.92e-07$ |  | $3.46e-05$ |  | $1.95e-06$ |
|  |  | $3.0$ |  | $4.0$ |  | $2.0$ |  | $3.0$
$3$ | $1.13e-01$ | $8.82e-07$ |  | $1.84e-08$ |  | $8.68e-06$ |  | $2.45e-07$ |
|  |  | $3.0$ |  | $4.0$ |  | $2.0$ |  | $3.0$
$4$ | $5.68e-02$ | $1.10e-07$ |  | $1.14e-09$ |  | $2.17e-06$ |  | $3.07e-08$ |
|  |  | $3.0$ |  | $4.0$ |  | $2.0$ |  | $3.0$
$5$ | $2.85e-02$ | $1.38e-08$ |  | $7.35e-11$ |  | $1.35e-07$ |  | $3.93e-09$ |

Table 2: Example 1. Relative errors and EOC.

For this example, we further investigate the possibility of choosing different
approximation orders. In particular, since the meshes we have considered to
generate Table 2 have the property $h_{\circ}=2h_{\partial}$, we analyse the
convergence of the scheme by fixing $k_{\partial}$ and varying $k_{\circ}$. In
Figures 5 and 6 we report the behaviour of the $H^{1}$-seminorm and
$L^{2}$-norm relative error, respectively. For each of them we fix in the left
plots $k_{\partial}=2$ and in the right ones $k_{\partial}=3$ and we report
the errors associated to the refinement levels 0,1,2 by varying $k_{\circ}$.
As we can see the CVEM convergence order dominates on the BEM one for each
$k_{\circ}\leq 4$, while for larger values the BEM error is no longer
negligible. Further, we observe that for $k_{\partial}=3$ and $k_{\circ}=5$
the CVEM $H^{1}$-convergence order is preserved, contrarily to that of the
$L^{2}$ one. It is worth to point out that, for fixed values of $h_{\circ}$
and $h_{\partial}$ and for fixed $k_{\partial}$, the maximum value of
$k_{\circ}$ such that the CVEM error is larger than that of the BEM is related
to the dependence of the implicit constants of the error estimates on
$k_{\circ}$ and $k_{\partial}$. We are aware of a study on such dependency for
some interior VEM problems (see for example [7], [8] and [9]). This is a task
by no mean trivial and is worth to be investigated. Finally, as we can see
from Table 1, while the increasing behaviour of $k_{\circ}$ for a fixed mesh
is approximately linear, that of lev. for a fixed $k_{\circ}$ is quadratic.
Therefore, it is worth noting that it is more efficient, in terms of
computational cost and memory saving, to use high order CVEM rather than to
refine the mesh, the latter choice being also computationally demanding for
what concerns the efficient computation of the BEM matrices.

Figure 5: Example 1. Behaviour of the $H^{1}$-seminorm relative error for
$k_{\partial}=2$ (left plot) and $k_{\partial}=3$ (right plot) by varying
$k_{\circ}$ (lev. 0,1,2).

Figure 6: Example 1. Behaviour of the $L^{2}$-norm relative error for
$k_{\partial}=2$ (left plot) and $k_{\partial}=3$ (right plot) by varying
$k_{\circ}$ (lev. 0,1,2).

###  7.2 Example 2

We consider the example proposed in [25] (and in [11]), for which $\Gamma_{0}$
is the boundary of the unit disk, centered at the origin of the cartesian
axis, $f=0$ and the datum $g$ on $\Gamma_{0}$ is defined as

$g(\mathbf{x})=\begin{cases}x_{1}^{4}&x_{1}\geq 0,\\\ 0&x_{1}<0.\end{cases}$

Solving the Dirichlet Laplace problem in polar coordinates, and expanding the
solution in terms of the eigenvectors of the associated Sturm Liouville
system, the solution in polar coordinates reads

$u(\rho,\theta)=\frac{3}{16}+\frac{\rho^{-2}}{4}\cos(2\theta)+\frac{\rho^{-4}}{16}\cos(4\theta)+\frac{48}{\pi}\sum_{\underset{n\text{~{}odd}}{n=1}}^{\infty}\frac{(-1)^{\nicefrac{{(n-1)}}{{2}}}\rho^{-n}}{n^{5}-20n^{3}+64n}\cos{(n\theta)},$

from which we deduce that the asymptotic behaviour is characterized by the
constant $\alpha=\nicefrac{{3}}{{16}}=0.1875$. We choose as artificial
boundary the ellipse of semi-axes $50$ and $15$, so that the values of the
numerical solution at the points $(-50,0)$ and $(50,0)$ can be considered good
approximations of $\alpha$. In Figure 7 we compare the behaviour of the exact
and numerical solutions in the intervals $[-50,-1]$ (left plot) and $[1,50]$
(right plot) for a fixed mesh of the computational domain and for different
choices of the approximation orders. Besides noting a very good agreement of
the solutions, we report that the corresponding absolute errors at $(-50,0)$
and $(50,0)$ are approximately equal to $4.0e-04$ for $k_{\circ}=1$ and
$k_{\partial}=2$, $5.0e-05$ for $k_{\circ}=k_{\partial}=2$ and $1.0e-08$ for
$k_{\circ}=k_{\partial}=3$.

Figure 7: Exact and numerical solutions in $[-50,-1]$ (left plot) and $[1,50]$
(right plot) by varying $k_{\circ}$ and $k_{\partial}$.

##  8 Conclusions

We have proposed and analysed the coupling of a Curved Virtual Element Method
with the one-equation Boundary Element Method to solve 2D exterior Poisson
problems. The peculiarity of the scheme consists in the use of decoupled
approximation orders for the interior CVEM and the boundary integral NRBC.
This strategy has allowed us to exploit the well-known flexibility of the CVEM
to retrieve an accurate solution by a low order approximation for the BEM.
Since high order BEMs require non-trivial computational efforts to efficiently
evaluate the matrix entries of the associated integral operators, the
advantage of using a low order BEM turns out to be a key aspect to achieve a
good accuracy and convergence rate weighted against computational costs.

The good performances obtained by applying the proposed scheme to elliptic
problems, encourage us to consider it within other contexts, such as time
dependent exterior problems, for which both the pure BEM and its coupling with
standard interior domain methods could become prohibitive.

## Declarations

This work was performed as part of the GNCS-INdAM 2020 research program
_“Metodologie innovative per problemi di propagazione di onde in domini
illimitati: aspetti teorici e computazionali”_.

The third author was partially supported by MIUR grant _“Dipartimenti di
Eccellenza 2018-2022”_ , CUP E11G18000350001.

## References

  * [1] P. F. Antonietti, G. Manzini, and M. Verani.  The conforming virtual element method for polyharmonic problems.  Comput. Math. Appl., 79(7):2021–2034, 2020.
  * [2] E. Artioli, S. Marfia, and E. Sacco.  VEM-based tracking algorithm for cohesive/frictional 2D fracture.  Comput. Methods Appl. Mech. Engrg., 365:112956, 21, 2020.
  * [3] 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.
  * [4] 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.
  * [5] L. Beirão da Veiga, F. Brezzi, L. D. Marini, and A. Russo.  Polynomial preserving virtual elements with curved edges.  Math. Models Methods Appl. Sci., 30(8):1555–1590, 2020.
  * [6] L. Beirão da Veiga, C. Canuto, R. H. Nochetto, and G. Vacca.  Equilibrium analysis of an immersed rigid leaflet by the virtual element method.  Math. Models Methods Appl. Sci., 31(7):1323–1372, 2021.
  * [7] L. Beirão da Veiga, A. Chernov, L. Mascotto, and A. Russo.  Basic principles of $hp$ virtual elements on quasiuniform meshes.  Math. Models Methods Appl. Sci., 26(8):1567–1598, 2016.
  * [8] L. Beirão da Veiga, A. Chernov, L. Mascotto, and A. Russo.  Exponential convergence of the $hp$ virtual element method in presence of corner singularities.  Numer. Math., 138(3):581–613, 2018.
  * [9] L. Beirão da Veiga, G. Manzini, and L. Mascotto.  A posteriori error estimation and adaptivity in $hp$ virtual elements.  Numer. Math., 143(1):139–175, 2019.
  * [10] 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.
  * [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.C. Brenner, Q. Guan, and L.Y. Sung.  Some estimates for virtual element methods.  Comput. Methods Appl. Math., 17(4):553–574, 2017.
  * [13] 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.
  * [14] F. Brezzi and C. Johnson.  On the coupling of boundary integral and finite element methods.  Calcolo, 16(2):189–201, 1979.
  * [15] M. Costabel.  Symmetric methods for the coupling of finite elements and boundary elements (invited contribution).  Comput. Mech., Southampton, 1987.
  * [16] L. Desiderio, S. Falletta, M. Ferrari, and L. Scuderi.  On the coupling of the Curved Virtual Element Method and the one-equation Boundary Element Method for 2d exterior Helmholtz problems.  arXiv:2107.04794 [math.NA], 2021.
  * [17] 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.
  * [18] 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.
  * [19] G. N. Gatica and S. Meddahi.  On the coupling of VEM and BEM in two and three dimensions.  SIAM J. Numer. Anal., 57(6):2493–2518, 2019.
  * [20] 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.
  * [21] 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.
  * [22] 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.
  * [23] G. C. Hsiao and W. L. Wendland.  Boundary Integral Equations, volume 164 of Applied Mathematical Sciences.  Springer, Berlin, 2008.
  * [24] 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.
  * [25] M. N. Le Roux.  Méthode d’éléments finis pour la résolution numérique de problèmes extérieurs en dimension $2$.  RAIRO Anal. Numér., 11(1):27–60, 112, 1977.
  * [26] L. Mascotto, I. Perugia, and A. Pichler.  A nonconforming Trefftz virtual element method for the Helmholtz problem.  Math. Models Methods Appl. Sci., 29(9):1619–1656, 2019.
  * [27] G. Monegato and L. Scuderi.  Numerical integration of functions with boundary singularities.  J. Comput. Appl. Math., 112(1-2):201–214, 1999.
  * [28] S. A. Sauter and C. Schwab.  Boundary element methods, volume 39 of Springer Series in Computational Mathematics, 2011.
  * [29] A. Sommariva and M. Vianello.  Product Gauss cubature over polygons based on Green’s integration formula.  BIT, 47(2):441–453, 2007.
  * [30] A. Sommariva and M. Vianello.  Gauss-Green cubature and moment computation over arbitrary geometries.  J. Comput. Appl. Math., 231(2):886–896, 2009.