# Hybrid coupling of finite element and boundary element methods using
Nitsche’s method and the Calderon projection

Timo Betcke email<EMAIL_ADDRESS>Department of Mathematics, University
College London, 25 Gordon Street, WC1H 0AY, London, United Kingdom Michał Bosy
email<EMAIL_ADDRESS>Department of Mathematics, University College London,
25 Gordon Street, WC1H 0AY, London, United Kingdom Erik Burman email:
<EMAIL_ADDRESS>Department of Mathematics, University College London, 25
Gordon Street, WC1H 0AY, London, United Kingdom

###### Abstract

In this paper we discuss a hybridised method for FEM-BEM coupling. The
coupling from both sides use a Nitsche type approach to couple to the trace
variable. This leads to a formulation that is robust and flexible with respect
to approximation spaces and can easily be combined as a building block with
other hybridised methods. Energy error norm estimates and the convergence of
Jacobi iterations are proved and the performance of the method is illustrated
on some computational examples.

Key words. FEM-BEM coupling, Nitsche’s method and hybridised methods

##  1 Introduction

The coupling of finite element (FEM) and boundary element (BEM) methods is the
most widely used approach for solving multi-physical problems on an unbounded
domain. It allows to take advantage of both methods. On the one hand, the BEM
reduces the dimension of the problem by using the boundary integral equation,
hence it is commonly used in exterior unbounded domains. On the other hand,
the FEM is known from their robustness and universal applicability even for
problem of inhomogeneous or non-linear nature.

The first coupled procedure was introduced by Zienkiewicz, Kelly and Bettess
[ZKB77]. It has been analysed by Brezzi, Johnson and Nédélec [BJN78], [BJ79]
and [JN80] for problem in unbounded domains. It is often referred to as the
Johnson-Nédélec coupling. Extension for higher order equations was considered
by Wendland [Wen86]. The convergence analysis requires compactness of the
double layer potential that can be obtain on the smooth boundary. Furthermore,
even for a symmetric discretisation scheme, the coupling method produces a
system of equations with a non-symmetric coefficient matrix.

In order to avoid these disadvantages, a symmetric coupling of FEM and BEM was
devised by Costabel [Cos87] and Han [Han90]. The independence of the
compactness condition was obtain by using both equations of the Calderón
system, contrary to the previously introduced methods that employ only one of
the two equations of the Calderón system. Some years later, Sayas [Say09]
showed that the weaker assumption of a Lipschitz coupling interface is
sufficient for the Johnson-Nédélec coupling. His analysis has since been
simplified by Steinbach [Ste11].

More recent developments have focused on the coupling of BEM with mixed FEM.
In [BCS96] and [MVMP96] the authors analysed symmetric coupling of BEM and
mixed FEM that uses Raviart-Thomas elements. Further work on coupling BEM with
mixed FEM with such elements as Brezzi-Douglas-Marini or Brezzi-Douglas-
Marini-Fortin spaces was analysis by Carstensen and Funken [CF00].

Gatica, Heuer and Sayas [GS06] and [GHS10] introduced the first the coupling
of BEM and discontinuous Galerkin (DG) methods, in order to exploit the
possibility to easily use high order approximation in the latter. Another
coupling of interior penalty DG methods with BEM was presented by Of, Rodin,
Steinbach and Taus [ORST12]. A general approach using the unified
hybridisation technique was presented by Cockburn and Sayas [CS12]. The class
of FEM considered includes the mixed, the DG and the hybridisable
discontinuous Galerkin (HDG) methods. Further collaboration of these authors
with Gúzman led to a new convergence results published in [CGS12].

In this paper, we present the coupling of FEM and BEM using weak imposition of
coupling conditions. Nitsche’s method [Nit71] is widely used in the context of
FEM for imposing boundary conditions. In addition, methods based on Nitsche’s
method have been successfully utilized for BEM domain decomposition problems
in [GHH09] and [CH12], and more recently for weakly imposing boundary
conditions for BEM in [BBS19]. Merging these two approaches for FEM and BEM we
weakly impose both coupling conditions in a hybridised formulation on the
boundary. The hybridisation is made by introducing a trace variable and
imposing the coupling in the form of a Nitsche type Dirichlet condition on the
two systems. The use of Nitsche’s method allows us to use the Dirichlet trace
as the hybrid variable, ensuring continuity by a consistent penalty term. The
test function partner of the trace variable, then acts to ensure continuity of
fluxes. The global system can be constructed using arbitrary approximation
orders for the two sub systems and the trace variable and the sub problem can
be solved independently. The stability of the method poses no constraint on
the approximation spaces and mesh refinement does not require special
treatment as in the case of Johnson-Nédélec coupling. This means that the two
systems can have independent meshes, that both have to be integrated only
against the trace variable. We here consider the standard continuous FEM, but
the formulation is by and large agnostic to the choice of FEM used in the bulk
and the method can be applied with discontinuous FEM as well, such as DG, HDG
[CGL09], or HHO [DPEL14] using a hybridised coupling on the interior domain
boundary. In the case of using discontinuous FEM, our formulation can be
interpreted as a hybridised interior penalty formulation of the class of
methods discussed in [CS12]. Finally we note that, thanks to the use of
Nitsche type mortaring, the method proposed herein can be used in the
framework for unfitted hybridised methods introduced in [BEH+19]. In that case
a surface mesh is required for the definition of the BEM method, but the FEM
approximation on the interior domain can be computed on an unfitted bulk
discretisation.

As many existing approaches of coupling FEM and BEM, we use Finite Element
Tearing And Interconnecting (FETI) and Boundary Element Tearing And
Interconnecting (BETI) type of methods [LS05] to solve the reduced system for
the hybrid variable. FETI is formulated using a Schur complement formulation,
while BETI is usually formulated in terms of Steklov-Poincaré operators.
Although Nitsche’s method is an established framework for domain decomposition
for finite elements methods such as FETI, it was not recognised by BETI
community. In this paper, we demonstrate how the hybrid Nitsche approach can
be integrated into the BETI framework.

The rest of the paper is organized as follows. We introduce the model problem
in this section. In Section 2, we present on continuous level the symmetric
coupling of BEM and FEM formulation known from [Cos87] and [Han90]. The
discrete formulation including weakly imposed coupling condition is introduced
and analysed in Section 3. Although the formulation obtained is not symmetric,
we comment of how symmetry can be obtained. We discuss iterative domain
decomposition in the model case of a simple relaxed Jacobi algorithm in
Section 4 and prove its convergence. In Section 5, we present some numerical
results and we conclude with some remarks in Section 6.

###  1.1 Model problem

Let us consider the unbounded domain $\Omega=\mathbb{R}^{3}$. We divide
$\Omega$ into a bounded internal part $\Omega^{-}$ and an unbounded external
part $\Omega^{+}$ with common Lipschitz boundary $\Gamma$, with $n$ the outer
unit normal vector of the domain $\Omega^{-}$ on $\Gamma$. We let
$\partial_{n}u:=\tfrac{\partial u}{\partial_{n}}$ denote the outward normal
derivative, $f\in L^{2}(\Omega)$ be a function with support in $\Omega^{-}$
and introduce a function $\varepsilon\in L^{\infty}(\Omega)$, $\epsilon\geq
0$. Then we can formulate our model problem as follows

(1.1) $\left\\{\begin{array}[]{rcll}-\Delta u^{-}+\varepsilon u&=&f&\mbox{in
}\Omega^{-},\\\ -\Delta u^{+}&=&0&\mbox{in }\Omega^{+},\\\
u^{-}&=&u^{+}&\mbox{on }\Gamma,\\\
\partial_{n}u^{-}&=&\partial_{n}u^{+}&\mbox{on }\Gamma,\\\
|u^{+}|&\rightarrow&0&\mbox{while }|x|\rightarrow\infty.\end{array}\right.$

where $u^{i}=u|_{\Omega^{i}}$ , $i\in\\{-,+\\}$. The function $\varepsilon$ is
introduced to make the interior problem heterogeneous and hence unsuitable for
treatment using the boundary element method alone.

######  Remark 1.

It is straightforward to extend the discussion to the case with a smoothly
varying diffusion coefficient in $\Omega^{-}$ which has a jump over $\Gamma$.

##  2 Variational formulation

Let $\left<\cdot,\cdot\right>_{\Gamma}$ denote the $L^{2}(\Gamma)$-inner
product that can be extended to a duality pairing on
$H^{-\frac{1}{2}}(\Gamma)\times H^{\frac{1}{2}}(\Gamma)$. We recall the
following result.

######  Lemma 1 (Duality pairing relation).

For any $\lambda\in H^{-\frac{1}{2}}(\Gamma)$ and $v\in
H^{\frac{1}{2}}(\Gamma)$, the following relation holds

(2.2)
$\left<\lambda,v\right>_{\Gamma}\leq\left\|\lambda\right\|_{H^{-\frac{1}{2}}(\Gamma)}\left\|v\right\|_{H^{\frac{1}{2}}(\Gamma)}.$

###### Proof.

It is obvious for $v=0$ and immediately by the definition of the dual norm

$\left\|\lambda\right\|_{H^{-\frac{1}{2}}(\Gamma)}=\sup_{0\neq v\in
H^{\frac{1}{2}}(\Gamma)}\tfrac{\left<\lambda,v\right>_{\Gamma}}{\left\|v\right\|_{H^{\frac{1}{2}}(\Gamma)}}.$

∎

We start with the variational formulation of the internal problem. Applying
integration by parts for first equation of (1.1) for every $v\in
H_{0}^{1}(\Omega^{-})$ we have

(2.3) $\int_{\Omega^{-}}\nabla u\cdot\nabla v\leavevmode\nobreak\
dx+\int_{\Omega^{-}}\varepsilon uv\leavevmode\nobreak\
dx-\left<\partial_{n}u,\
v\right>_{\Gamma}=\int_{\Omega^{-}}fv\leavevmode\nobreak\ dx.$

We define the Green’s function for the Laplace operator in $\mathbb{R}^{3}$ as
follows

$G(x,y):=\tfrac{1}{4\pi|x-y|}.$

In this paper, we focus on the problem in $\mathbb{R}^{3}$. A similar analysis
can be used for problems in $\mathbb{R}^{2}$, in which case this definition
should be replaced by $G(x,y):=\frac{\log{|x-y|}}{2\pi}$. Following the
standard approach (see, e.g. [Ste08, Chapter 6]), we introduce single layer
and double layer operators $\mathcal{V}:H^{-\frac{1}{2}}(\Gamma)\rightarrow
H^{1}(\Omega^{+})$ and $\mathcal{K}:H^{\frac{1}{2}}(\Gamma)\rightarrow
H^{1}(\Omega^{+})$ respectively as

$\displaystyle(\mathcal{V}\varphi)(x)$
$\displaystyle:=\int_{\Gamma}G(x,y)\varphi(y)\leavevmode\nobreak\ dy$
$\displaystyle\mbox{for }\varphi\in H^{-\frac{1}{2}}(\Gamma),$
$\displaystyle(\mathcal{K}v)(x)$ $\displaystyle:=\int_{\Gamma}\tfrac{\partial
G(x,y)}{\partial n_{y}}v(y)\leavevmode\nobreak\ dy$ $\displaystyle\mbox{for
}v\in H^{\frac{1}{2}}(\Gamma),$

where $x\in\Omega^{+}\setminus\Gamma$ and $n_{y}$ is an outer unit normal
vector (for $\Omega^{i}$, $i\in\\{-,+\\}$) in the point $y$.

Following [Ste08, Chapter 1], we define the Dirichlet and Neumann traces

$\displaystyle\gamma_{D}^{i}:H^{1}(\Omega^{i})$ $\displaystyle\rightarrow
H^{\frac{1}{2}}(\Gamma)$ $\displaystyle\gamma_{D}^{i}f(x)$
$\displaystyle:=\lim_{\Omega^{i}\ni y\rightarrow x\in\Gamma}f(y),$
$\displaystyle\gamma_{N}^{i}:H^{1}(\Delta,\Omega^{i})$
$\displaystyle\rightarrow H^{-\frac{1}{2}}(\Gamma)$
$\displaystyle\gamma_{N}^{i}f(x)$ $\displaystyle:=\lim_{\Omega^{i}\ni
y\rightarrow x\in\Gamma}n_{x}\cdot\nabla f(y),$

where $H^{1}(\Delta,\Omega^{i}):=\left\\{v\in H^{1}(\Omega^{i}):\Delta v\in
L^{2}(\Omega^{i})\right\\}$, for $i\in\\{-,+\\}$, and $n_{x}$ is an outer (for
$\Omega^{-}$) normal vector to $\Gamma$ in the point $x$. The following
results will be useful in what follows.

######  Lemma 2 (Trace theorem).

Let $i\in\\{-,+\\}$, then for $\Omega^{i}\subset\mathbb{R}^{3}$ the trace
operator $\gamma_{D}^{i}:H^{1}(\Omega^{i})\rightarrow H^{\frac{1}{2}}(\Gamma)$
is bounded for all $v\in H^{1}(\Omega^{i})$

(2.4) $\|\gamma_{D}^{i}v\|_{H^{\frac{1}{2}}(\Gamma)}\leq
C_{T}\|v\|_{H^{1}(\Omega^{i})}.$

###### Proof.

See [Ste08, Theorem 2.21]. ∎

We use $\\{\cdot\\}_{\Gamma}$ to denote an average of the interior and
exterior traces of a function. Then, applying the trace mappings yields to the
single layer, double layer, adjoint double layer potentials and hypersingular
boundary integral operator

$\displaystyle V:H^{-\frac{1}{2}}(\Gamma)$ $\displaystyle\rightarrow
H^{\frac{1}{2}}(\Gamma)$ $\displaystyle V$
$\displaystyle:=\\{\gamma_{D}\mathcal{V}\\}_{\Gamma},$ $\displaystyle
K:H^{\frac{1}{2}}(\Gamma)$ $\displaystyle\rightarrow H^{\frac{1}{2}}(\Gamma)$
$\displaystyle K$ $\displaystyle:=\\{\gamma_{D}\mathcal{K}\\}_{\Gamma},$
$\displaystyle K^{\prime}:H^{-\frac{1}{2}}(\Gamma)$ $\displaystyle\rightarrow
H^{-\frac{1}{2}}(\Gamma)$ $\displaystyle K^{\prime}$
$\displaystyle:=\\{\gamma_{N}\mathcal{V}\\}_{\Gamma},$ $\displaystyle
W:H^{\frac{1}{2}}(\Gamma)$ $\displaystyle\rightarrow H^{-\frac{1}{2}}(\Gamma)$
$\displaystyle W$ $\displaystyle:=\\{\gamma_{N}\mathcal{K}\\}_{\Gamma}.$

For the solution $u$ of the problem (1.1), we have the following boundary
integral equations on $\Gamma$

(2.5) $\displaystyle\begin{pmatrix}\gamma_{D}^{-}u\\\
\gamma_{N}^{-}u\end{pmatrix}=C^{-}\begin{pmatrix}\gamma_{D}^{-}u\\\
\gamma_{N}^{-}u\end{pmatrix}$ ,
$\displaystyle\begin{pmatrix}\gamma_{D}^{+}u\\\
\gamma_{N}^{+}u\end{pmatrix}=C^{+}\begin{pmatrix}\gamma_{D}^{+}u\\\
\gamma_{N}^{+}u\end{pmatrix},$

where $C^{\pm}:H^{\frac{1}{2}}(\Gamma)\times
H^{-\frac{1}{2}}(\Gamma)\rightarrow H^{\frac{1}{2}}(\Gamma)\times
H^{-\frac{1}{2}}(\Gamma)$ denotes two Calderón projectors defined as follows

$C^{\pm}:=\begin{pmatrix}\tfrac{1}{2}Id\pm K&\mp V\\\ \mp W&\tfrac{1}{2}Id\mp
K^{\prime}\end{pmatrix}.$

From the relation (2.5) for external traces we can construct the following
exterior Dirichlet-to-Neumann operator

(2.6) $DtN^{+}:=-W+(\tfrac{1}{2}Id-K^{\prime})\circ
V^{-1}\circ(K-\tfrac{1}{2}Id).$

Obviously, it makes sense only if the inverse of the operator $V$ exists.

Using Dirichlet-to-Neumann operator (2.6) we introduce a new variable
$\lambda=\gamma_{N}^{+}u=\partial_{n}u^{+}$ as

(2.7) $\lambda:=\left(V^{-1}\circ(K-\tfrac{1}{2}Id)\right)\gamma_{D}^{-}u.$

The classical symmetric coupling that satisfies the transmission conditions of
(1.1) is as follow

Find $u\in H^{1}(\Omega^{-})$ and $\lambda\in H^{-\frac{1}{2}}(\Gamma)$ such
that for all $v\in H^{1}(\Omega^{-})$ and $\zeta\in H^{-\frac{1}{2}}(\Gamma)$

(2.8) $\left\\{\begin{array}[]{rcl}\int_{\Omega^{-}}\nabla u\cdot\nabla
v\leavevmode\nobreak\ dx+\int_{\Omega^{-}}\varepsilon uv\leavevmode\nobreak\
dx+\left<Wu,\
v\right>_{\Gamma}-\left<\left(\tfrac{1}{2}Id-K^{\prime}\right)\lambda,\
v\right>_{\Gamma}&=&\int_{\Omega^{-}}fv\leavevmode\nobreak\ dx,\\\\[8.53581pt]
\left<\left(\tfrac{1}{2}Id-K\right)u,\ \zeta\right>_{\Gamma}+\left<V\lambda,\
\zeta\right>_{\Gamma}&=&0.\end{array}\right.$

###  2.1 Well-posedness of the continuous problem

The following results are well known (see [Cos87] or [Han90]), but for
reader’s convenience we present them in the case of problem (1.1). Let us
propose a more compact formulation of (2.8).

Find $u\in H^{1}(\Omega^{-})$ and $\lambda\in H^{-\frac{1}{2}}(\Gamma)$ such
that for all $v\in H^{1}(\Omega^{-})$ and $\zeta\in H^{-\frac{1}{2}}(\Gamma)$

(2.9)
$A\left(\left(u,\lambda\right),\left(v,\zeta\right)\right)=\int_{\Omega^{-}}fv\leavevmode\nobreak\
dx,$

where

$\displaystyle A\left(\left(w,\lambda\right),\left(v,\zeta\right)\right):=$
$\displaystyle\int_{\Omega^{-}}\nabla w\cdot\nabla v\leavevmode\nobreak\
dx+\int_{\Omega^{-}}\varepsilon uv\leavevmode\nobreak\
dx-\tfrac{1}{2}\left<\lambda,\ v\right>_{\Gamma}+\tfrac{1}{2}\left<w,\
\zeta\right>_{\Gamma}$ $\displaystyle+\left<Ww,\
v\right>_{\Gamma}+\left<K^{\prime}\lambda,\ v\right>_{\Gamma}-\left<Kw,\
\zeta\right>_{\Gamma}+\left<V\lambda,\ \zeta\right>_{\Gamma}.$

For simplicity we introduce the space $\mathbb{V}:=H^{1}(\Omega^{-})\times
H^{-\frac{1}{2}}(\Gamma)$ and the associated norm

(2.10)
$\left\|(v,\varphi)\right\|_{\mathbb{V}}^{2}:=\left\|v\right\|_{H^{1}(\Omega^{-})}^{2}+\left\|\varphi\right\|_{H^{-\frac{1}{2}}(\Gamma)}^{2}.$

######  Lemma 3 (Continuity).

There exists constant $\beta>0$ such that for all $w,v\in H^{1}(\Omega^{-})$
and $\lambda,\varphi\in H^{-\frac{1}{2}}(\Gamma)$

(2.11)
$\left|A\left(\left(w,\lambda\right),\left(v,\varphi\right)\right)\right|\leq\beta\left\|(w,\lambda)\right\|_{\mathbb{V}}\left\|(v,\varphi)\right\|_{\mathbb{V}}.$

###### Proof.

We use the Cauchy-Schwarz inequality, the relation (2.2) and continuity of
boundary operators (see [Ste08, Section 6.2-6.5]) to obtain

$\displaystyle\left|A\left(\left(w,\lambda\right),\left(v,\varphi\right)\right)\right|\leq$
$\displaystyle\max\\{1,\|\varepsilon\|_{L^{\infty}(\Omega)}\\}\left\|w\right\|_{H^{1}(\Omega^{-})}\left\|v\right\|_{H^{1}(\Omega^{-})}$
$\displaystyle+\tfrac{1}{2}\left\|w\right\|_{H^{\frac{1}{2}}(\Gamma)}\left\|\varphi\right\|_{H^{-\frac{1}{2}}(\Gamma)}+\tfrac{1}{2}\left\|\lambda\right\|_{H^{-\frac{1}{2}}(\Gamma)}\left\|v\right\|_{H^{\frac{1}{2}}(\Gamma)}$
$\displaystyle+C_{K}\left\|w\right\|_{H^{\frac{1}{2}}(\Gamma)}\left\|\varphi\right\|_{H^{-\frac{1}{2}}(\Gamma)}+C_{V}\left\|\lambda\right\|_{H^{-\frac{1}{2}}(\Gamma)}\left\|\varphi\right\|_{H^{-\frac{1}{2}}(\Gamma)}$
$\displaystyle+C_{K^{\prime}}\left\|\lambda\right\|_{H^{-\frac{1}{2}}(\Gamma)}\left\|v\right\|_{H^{\frac{1}{2}}(\Gamma)}+C_{W}\left\|w\right\|_{H^{\frac{1}{2}}(\Gamma)}\left\|v\right\|_{H^{\frac{1}{2}}(\Gamma)}.$

We finished by applying the trace inequality (2.4) to terms including
$\left\|\cdot\right\|_{H^{\frac{1}{2}}(\Gamma)}$ norm,

$\left\|v\right\|_{H^{\frac{1}{2}}(\Gamma)}+\left\|w\right\|_{H^{\frac{1}{2}}(\Gamma)}\leq
C_{T}\left(\left\|v\right\|_{H^{1}(\Omega^{-})}+\left\|w\right\|_{H^{1}(\Omega^{-})}\right).$

∎

######  Lemma 4 (Coercivity).

There exists constant $\alpha>0$ such that for all $v\in H^{1}(\Omega^{-})$
and $\varphi\in H^{-\frac{1}{2}}(\Gamma)$

(2.12)
$A\left(\left(v,\varphi\right),\left(v,\varphi\right)\right)\geq\alpha\left\|(v,\varphi)\right\|_{\mathbb{V}}^{2}.$

###### Proof.

Using coercivity of $V$ (see [Ste08, Theorem 6.22]) and coercivity of $W$ (see
[Ste08, Theorem 6.24]) we obtain

$\displaystyle
A\left(\left(v,\varphi\right),\left(v,\varphi\right)\right)\geq$
$\displaystyle
c_{\epsilon}\left\|v\right\|_{H^{1}(\Omega^{-})}^{2}+\alpha_{V}\left\|\varphi\right\|_{H^{-\frac{1}{2}}(\Gamma)}^{2}+\alpha_{W}\left\|v-\bar{v}\right\|_{H^{\frac{1}{2}}(\Gamma)}^{2}$

where $\bar{v}$ denotes the average over $\Gamma$ of $v$ and
$c_{\varepsilon}=\min(1,\varepsilon)$. This shows (2.12) when $\varepsilon>0$.
For the case $\epsilon=0$ we need a Poincaré inequality of the form

(2.13) $c_{P}\|v\|_{H^{1}(\Omega)}\leq\|\nabla v\|_{L^{2}(\Omega)}+|g(v)|$

where $g(v)$ is some functional that is non-zero for constant non-zero $v$
[EG04, Lemma B.63]. We claim that this holds with
$g(v)=\left<V\varphi(v),\varphi(v)\right>_{\Gamma}$ where $\varphi(v)$ is
defined by the second equation of (2.8). We immediately see that if this is
true then there exists $\alpha>0$ such that

$\displaystyle
A\left(\left(v,\varphi\right),\left(v,\varphi\right)\right)\geq$
$\displaystyle\alpha(\left\|v\right\|_{H^{1}(\Omega^{-})}^{2}$

We need to show that for constant $v\neq 0$ there holds
$\left<V\varphi(v),\varphi(v)\right>_{\Gamma}\neq 0$. Let $v|_{\Gamma}=u_{0}$
be a non-zero constant. Then we need to study

$\left<\left(\tfrac{1}{2}Id-K\right)u_{0},\
\zeta\right>_{\Gamma}+\left<V\varphi(v),\
\zeta\right>_{\Gamma}=0,\quad\forall\zeta\in H^{-\frac{1}{2}}(\Gamma).$

We argue by contradiction. Assume that $\varphi(v)=0$. Then $u_{0}$ is the
trace of a solution to the homogeneous Neumann problem in $\Omega^{-}$.
However by the first line of the left relation of equation (2.5) there holds
for all such traces

(2.14) $\left<\left(\tfrac{1}{2}Id-K\right)u_{0},\
\zeta\right>_{\Gamma}=\left<Id\,u_{0},\
\zeta\right>_{\Gamma},\quad\forall\zeta\in H^{-\frac{1}{2}}(\Gamma).$

Hence the functional is defined by

$\left<V\varphi(v),\ \zeta\right>_{\Gamma}=-\left<Id\,u_{0},\
\zeta\right>_{\Gamma},\quad\forall\zeta\in H^{-\frac{1}{2}}(\Gamma).$

However since the operator on the left hand side, defined by, $V$ is injective
it follows that $\varphi(v)\neq 0$, which leads to a contradiction. Hence
$\left<V\varphi(v),\varphi(v)\right>_{\Gamma}\neq 0$ for $v$ constant. This
concludes the proof. ∎

The existence and uniqueness of the solution of problem (2.9) is achieved by
using the Lax-Milgram theorem.

##  3 The discrete problem

We assume that $\Omega^{-}$ is a polyhedral domain. The boundary $\Gamma$ may
be decomposed in a set of $n$ planar surfaces $\\{\Gamma_{j}\\}_{j=1}^{n}$. It
will be convenient to use following broken Sobolev spaces over the polyhedral
boundary $\Gamma$. For $s>1$ define

$\widetilde{H}^{s}(\Gamma):=\left\\{v\in H^{1}(\Gamma):\ v|_{\Gamma_{j}}\in
H^{s}(\Gamma_{j}),\,1\leq j\leq n\right\\}.$

When $s>1$ the norm on $\widetilde{H}^{s}(\Gamma)$, is defined as the broken
norm over the faces of the polyhedral boundary $\Gamma$

$\|v\|_{\widetilde{H}^{s}(\Gamma)}:=\left(\sum_{j=1}^{n}\|v\|^{2}_{H^{s}(\Gamma_{j})}\right)^{\frac{1}{2}}.$

When $0\leq s\leq 1$ the space $\widetilde{H}^{s}(\Gamma)$ coincides with the
usual space $H^{s}(\Gamma)$ and their norms are the same (for more details see
[SS11, Definition 4.1.48]).

Let $\mathcal{T}_{h}$ be a triangulation of $\overline{\Omega}^{-}$ made of
tetrahedrons. For each triangulation $\mathcal{T}_{h}$, $\mathcal{E}_{h}$
denotes the set of its facets. In addition, for each of element
$K\in\mathcal{T}_{h}$, $h_{K}:=\mbox{diam}(K)$, and
$h:=\max_{K\in\mathcal{T}_{h}}h_{K}$. Let $\mathcal{G}_{i}$, $i=1,2$ denote
two different surface triangulation of the boundary $\Gamma$. For notation
convenience we here assume that the trace mesh of $\mathcal{T}_{h}$ and the
$\mathcal{G}_{i}$ all have the similar local mesh size.

The following result will be useful in what follows.

######  Lemma 5 (Trace inequality).

There exists $C_{max}>0$, independent of $h_{K}$, such that for all
$K\in\mathcal{T}_{h}$ and polynomial function $v$ in $K$ the following
discrete trace inequality holds

(3.15) $\displaystyle h_{K}^{\frac{1}{2}}\|v\|_{L^{2}(\partial K)}$
$\displaystyle\leq C_{max}\|{v}\|_{L^{2}(K)}.$

###### Proof.

See [DPE12, Lemma 1.46] ∎

To discretise the problem (2.8) over the triangulation one can choose either
continuous or discontinuous finite elements. For simplicity of the analysis we
choose the following spaces

$\displaystyle V_{h}^{j}$ $\displaystyle:=\left\\{{v}_{h}\in
C^{0}\left(\Omega^{-}\right):\ {v}_{h}|_{K}\in\mathbb{P}_{j}\left(K\right)\
\forall\ {K\in\mathcal{T}_{h}}\right\\},$ $\displaystyle W_{h}^{k}$
$\displaystyle:=\left\\{{w}_{h}\in C^{0}\left(\Gamma\right):\
{w}_{h}|_{E}\in\mathbb{P}_{k}\left(E\right)\ \forall\
E\in\mathcal{G}_{1}\right\\},$ $\displaystyle\Lambda_{h}^{l}$
$\displaystyle:=\left\\{\lambda_{h}\in L^{2}\left(\Gamma\right):\
\lambda_{h}|_{E}\in\mathbb{P}_{l}\left(E\right)\ \forall\
E\in\mathcal{G}_{1}\right\\},$ $\displaystyle M_{h}^{m}$
$\displaystyle:=\left\\{\widetilde{v}_{h}\in L^{2}(\Gamma):\ \
\widetilde{v}_{h}|_{E}\in\mathbb{P}_{m}\left(E\right)\ \forall\
E\in\mathcal{G}_{2}\right\\}.$

Let us denote $\mathcal{V}_{h}:=V_{h}^{j}\times W_{h}^{k}$ and
$v_{h}=\left(v_{h}^{-},v_{h}^{+}\right)\in\mathcal{V}_{h}$. Using the above
spaces we propose the hybrid discrete formulation of the problem (2.8)

Find
$\left(u_{h},\lambda_{h},\widetilde{u}_{h}\right)\in\mathcal{V}_{h}\times\Lambda_{h}^{l}\times
M_{h}^{m}$ such that for all
$\left(v_{h},\varphi_{h},\widetilde{v}_{h}\right)\in\mathcal{V}_{h}\times\Lambda_{h}^{l}\times
M_{h}^{m}$

(3.16)
$a_{h}\left(\left(u_{h}^{-},\widetilde{u}_{h}\right),\left(v_{h}^{-},\widetilde{v}_{h}\right)\right)+b_{h}\left(\left(u_{h}^{+},\lambda_{h},\widetilde{u}_{h}\right),\left(v_{h}^{+},\varphi_{h},\widetilde{v}_{h}\right)\right)=\int_{\Omega^{-}}fv_{h}^{-}dx,$

where

$\displaystyle
a_{h}\left(\left(w_{h},\widetilde{w}_{h}\right),\left(v_{h},\widetilde{v}_{h}\right)\right):=$
$\displaystyle\int_{\Omega^{-}}\nabla w_{h}:\nabla
v_{h}dx+\int_{\Omega^{-}}\varepsilon w_{h}v_{h}dx$
$\displaystyle-\left<\partial_{n}w_{h},\
v_{h}-\widetilde{v}_{h}\right>_{\Gamma}-\left<w_{h}-\widetilde{w}_{h},\
\partial_{n}v_{h}\right>_{\Gamma}$
$\displaystyle+\tfrac{\tau}{h}\left<w_{h}-\widetilde{w}_{h},\
v_{h}-\widetilde{v}_{h}\right>_{\Gamma},$ $\displaystyle
b_{h}\left(\left(w_{h},\lambda_{h},\widetilde{w}_{h}\right),\left(v_{h},\varphi_{h},\widetilde{v}_{h}\right)\right):=$
$\displaystyle\left<\left(\tfrac{1}{2}Id-K\right)w_{h},\
\varphi_{h}\right>_{\Gamma}+\left<V\lambda_{h},\ \varphi_{h}\right>_{\Gamma}$
$\displaystyle+\left<Ww_{h},\
v_{h}\right>_{\Gamma}-\left<\left(\tfrac{1}{2}Id-K^{\prime}\right)\lambda_{h},\
v_{h}\right>_{\Gamma}$ $\displaystyle+\left<\lambda_{h},\
v_{h}-\widetilde{v}_{h}\right>_{\Gamma}-\left<w_{h}-\widetilde{w}_{h},\
\varphi_{h}\right>_{\Gamma}$
$\displaystyle+\tfrac{\tau}{h}\left<w_{h}-\widetilde{w}_{h},\
v_{h}-\widetilde{v}_{h}\right>_{\Gamma}.$

The stabilisation parameter $\tau>0$ has to be chosen appropriately. The
formulation of bilinear form $a_{h}$ is well known for example from [Egg09].
As we said before it is possible to use the discontinuous finite element
method for example symmetric interior penalty hybrid discontinuous Galerkin
method presented in [ES10].

######  Remark 2 (Impedance boundary condition).

The hybrid weakly imposed Dirichlet and Neumann boundary conditions is related
to an impedance boundary condition of the type

$\widetilde{u}=-\gamma\frac{\partial u}{\partial n}+u,$

with $\gamma=\frac{h}{\tau}$. This can be seen considering terms associated
with $\widetilde{v}_{h}$ in above definition of the bilinear forms.

######  Remark 3 (Relation to a standard Nitsche type method without
hybridisation).

The trace variable $\widetilde{u}_{h}$ can be eliminated by replacing it by a
linear combination of $u_{h}^{+}$ and $u_{h}^{-}$ with similarly the test
function $\widetilde{v}_{h}$ replaced by the same linear combination of the
test functions $v_{h}^{+}$ and $v_{h}^{-}$ (see [BHL21, Section 4.2]). The
below analysis carries over verbatim to this case. If in addition $v_{h}^{+}$
and $v_{h}^{-}$ are chosen in the same spaces (with the same trace meshes)
$u_{h}^{-}$ and $v_{h}^{-}$ can be substituted for $u_{h}^{+}$ and $v_{h}^{+}$
respectively, resulting in a (fully coupled) method with the unknowns
$u_{h}^{-}$ and $\lambda_{h}$ only.

###  3.1 Symmetric formulation

Despite using the symmetric Nitsche method, our whole system is not symmetric.
This is a consequence of the lack of symmetry of the boundary element method
with weak imposition. We can use the Steklov-Poincaré operator to eliminate
the flux variable, so that the non-symmetric method above is transformed into
a symmetric reduced system as we show below.

The following equations are associated with bilinear form $b_{h}$ from (3.16)
reads

$\displaystyle-\left<\left(\tfrac{1}{2}Id-K\right)w_{h},\varphi_{h}\right>_{\Gamma}-\left<V\lambda_{h},\varphi_{h}\right>_{\Gamma}$
$\displaystyle=-\left<w_{h}-\widetilde{w}_{h},\varphi_{h}\right>_{\Gamma},$
$\displaystyle-\left<Ww_{h},v_{h}\right>_{\Gamma}+\left<\left(\tfrac{1}{2}Id-K^{\prime}\right)\lambda_{h},v_{h}\right>_{\Gamma}-\left<\lambda_{h},v_{h}-\widetilde{v}_{h}\right>_{\Gamma}$
$\displaystyle=\tfrac{\tau}{h}\left<w_{h}-\widetilde{w}_{h},v_{h}-\widetilde{v}_{h}\right>_{\Gamma}.$

Similar to the continuous formulation we use the Dirichlet-to-Neumann operator
(2.6) to obtain

(3.17)
$\lambda_{h}:=\left(V^{-1}\circ(K-\tfrac{1}{2}Id)\right)w_{h}+V^{-1}\left(w_{h}-\widetilde{w}_{h}\right).$

Injecting this relation into the second equation leads to the formulation of
the new symmetric bilinear form

$\displaystyle\widehat{b}_{h}\left(\left(w_{h},\widetilde{w}_{h}\right),\left(v_{h},\widetilde{v}_{h}\right)\right):=$
$\displaystyle\left<Ww_{h},\
v_{h}\right>_{\Gamma}-\left<\left(\tfrac{1}{2}Id-K^{\prime}\right)V^{-1}(K-\tfrac{1}{2}Id)w_{h},\
v_{h}\right>_{\Gamma}$
$\displaystyle-\left<\left(\tfrac{1}{2}Id-K^{\prime}\right)V^{-1}\left(w_{h}-\widetilde{w}_{h}\right),\
v_{h}\right>_{\Gamma}$ $\displaystyle+\left<V^{-1}(K-\tfrac{1}{2}Id)w_{h},\
v_{h}-\widetilde{v}_{h}\right>_{\Gamma}$
$\displaystyle+\left<V^{-1}\left(w_{h}-\widetilde{w}_{h}\right),\
v_{h}-\widetilde{v}_{h}\right>_{\Gamma}+\tfrac{\tau}{h}\left<w_{h}-\widetilde{w}_{h},\
v_{h}-\widetilde{v}_{h}\right>_{\Gamma}.$

###  3.2 Well-posedness of the discrete problem

Let us consider the following norms

$\displaystyle\left\|\left({w_{h}},{\widetilde{w}_{h}}\right)\right\|_{\mathcal{F}_{*}}^{2}:=$
$\displaystyle\|{w_{h}}\|_{H^{1}(\Omega^{-})}^{2}+\tfrac{\tau}{h}\left\|{w_{h}}-{\widetilde{w}_{h}}\right\|_{L^{2}(\Gamma)}^{2},$
(3.18)
$\displaystyle\left\|\left({w_{h}},{\widetilde{w}_{h}}\right)\right\|_{\mathcal{F}}^{2}:=$
$\displaystyle\left\|\left({w_{h}},{\widetilde{w}_{h}}\right)\right\|_{\mathcal{F}_{*}}^{2}+h\left\|{\partial_{n}w_{h}}\right\|_{L^{2}(\Gamma)}^{2},$
$\displaystyle\left\|\left({w_{h}},\lambda_{h},{\widetilde{w}_{h}}\right)\right\|_{\mathcal{B}_{*}}^{2}:=$
$\displaystyle\|w_{h}\|_{H^{\frac{1}{2}}(\Gamma)}^{2}+\|\lambda_{h}\|_{H^{-\frac{1}{2}}(\Gamma)}^{2}+\tfrac{\tau}{h}\left\|{w_{h}}-{\widetilde{w}_{h}}\right\|_{L^{2}(\Gamma)}^{2},$
(3.19)
$\displaystyle\left\|\left({w_{h}},\lambda_{h},{\widetilde{w}_{h}}\right)\right\|_{\mathcal{B}}^{2}:=$
$\displaystyle\left\|\left({w_{h}},\lambda_{h},{\widetilde{w}_{h}}\right)\right\|_{\mathcal{B}_{*}}^{2}+h\left\|\lambda_{h}\right\|_{L^{2}(\Gamma)}^{2}.$

######  Lemma 6 (Equivalence of the norms).

For all
$\left({w_{h}},\lambda_{h},{\widetilde{w}_{h}}\right)\in\mathcal{V}_{h}\times\Lambda_{h}^{l}\times
M_{h}^{m}$ there exist positive constants $C_{\mathcal{F}},C_{\mathcal{B}}$
such that

(3.20)
$\displaystyle\left\|\left({w_{h}}^{-},{\widetilde{w}_{h}}\right)\right\|_{\mathcal{F}_{*}}\leq\left\|\left({w_{h}}^{-},{\widetilde{w}_{h}}\right)\right\|_{\mathcal{F}}$
$\displaystyle\leq
C_{\mathcal{F}}\left\|\left({w_{h}}^{-},{\widetilde{w}_{h}}\right)\right\|_{\mathcal{F}_{*}},$
(3.21)
$\displaystyle\left\|\left({w_{h}}^{+},\lambda_{h},{\widetilde{w}_{h}}\right)\right\|_{\mathcal{B}_{*}}\leq\left\|\left({w_{h}}^{+},\lambda_{h},{\widetilde{w}_{h}}\right)\right\|_{\mathcal{B}}$
$\displaystyle\leq
C_{\mathcal{B}}\left\|\left({w_{h}}^{+},\lambda_{h},{\widetilde{w}_{h}}\right)\right\|_{\mathcal{B}_{*}}.$

###### Proof.

For (3.20) we use the trace inequality (3.15) and for (3.21) we use the
inverse inequality
$h^{\frac{1}{2}}\left\|\lambda_{h}\right\|_{L^{2}(\Gamma)}\leq
C\left\|\lambda_{h}\right\|_{H^{-\frac{1}{2}}(\Gamma)}$ (see [SS11, Remark
4.4.4]). ∎

######  Lemma 7 (Continuity).

There exists positive constant $\beta_{\mathcal{F}}$ such that for all $w,v\in
H^{\frac{3}{2}+\delta}(\Omega^{-})$, for $\delta>0$, and
$\widetilde{w},\widetilde{v}\in L^{2}(\Gamma)$

(3.22)
$\left|a_{h}\left(\left(w,\widetilde{w}\right),\left(v,\widetilde{v}\right)\right)\right|\leq\beta_{\mathcal{F}}\left\|(w,\widetilde{w})\right\|_{\mathcal{F}}\left\|(v,\widetilde{v})\right\|_{\mathcal{F}}.$

There exists positive constant $\beta_{\mathcal{B}}$ such that for all $w,v\in
H^{\frac{1}{2}}(\Gamma)$, $\lambda,\varphi\in L^{2}(\Gamma)$ and
$\widetilde{w},\widetilde{v}\in L^{2}(\Gamma)$

(3.23)
$\left|b_{h}\left(\left(w,\lambda,\widetilde{w}\right),\left(v,\varphi,\widetilde{v}\right)\right)\right|\leq\beta_{\mathcal{B}}\left\|(w,\lambda,\widetilde{w})\right\|_{\mathcal{B}}\left\|(v,\varphi,\widetilde{v})\right\|_{\mathcal{B}}.$

###### Proof.

We use Cauchy-Schwarz inequality to obtain (3.22). In the case of equation
(3.23), we use Cauchy-Schwarz inequality and Lemma 3. ∎

To show the well-posedness of (3.16) we need the ellipticity of the bilinear
forms $a_{h}$ and $b_{h}$

######  Lemma 8 (Coercivity).

Assume that positive constant $\tau$ is large enough. Then, there exists
positive constant $\alpha$ such that for all
$\left({w_{h}},\lambda_{h},{\widetilde{w}_{h}}\right)\in\mathcal{V}_{h}\times\Lambda_{h}^{l}\times
M_{h}^{m}$

(3.24) $\displaystyle
a_{h}\left(\left(w_{h}^{-},\widetilde{w}_{h}\right),\left(w_{h}^{-},\widetilde{w}_{h}\right)\right)+$
$\displaystyle
b_{h}\left(\left(w_{h}^{+},\lambda_{h},{\widetilde{w}_{h}}\right),\left({w_{h}}^{+},\lambda_{h},{\widetilde{w}_{h}}\right)\right)$
$\displaystyle\geq\alpha\left(\left\|(w_{h},\widetilde{w}_{h})\right\|_{\mathcal{F}}^{2}+\left\|\left({w_{h}},\lambda_{h},{\widetilde{w}_{h}}\right)\right\|_{\mathcal{B}}^{2}\right).$

###### Proof.

Let us start with bilinear form $a_{h}$. First we assume that
$\varepsilon\geq\varepsilon_{min}>0$

$\displaystyle
a_{h}\left(\left(w_{h}^{-},\widetilde{w}_{h}\right),\left(w_{h}^{-},\widetilde{w}_{h}\right)\right)=$
$\displaystyle|w_{h}^{-}|^{2}_{H^{1}(\Omega^{-})}+\|\varepsilon^{1/2}w_{h}^{-}\|_{\Omega^{-}}^{2}-2\left<\partial_{n}w_{h}^{-},w_{h}^{-}-\widetilde{w}_{h}\right>_{\Gamma}$
$\displaystyle+\tfrac{\tau}{h}\left\|w_{h}^{-}-\widetilde{w}_{h}\right\|_{L^{2}(\Gamma)}^{2}.$

Using Cauchy-Schwarz and trace inequalities (3.15), followed by Young’s
inequality, we arrive at

$\displaystyle
a_{h}\left(\left(w_{h}^{-},\widetilde{w}_{h}\right),\left(w_{h}^{-},\widetilde{w}_{h}\right)\right)\geq$
$\displaystyle|w_{h}^{-}|^{2}_{H^{1}(\Omega^{-})}+\varepsilon_{\min}\|w_{h}^{-}\|_{\Omega^{-}}^{2}+\tfrac{\tau}{h}\left\|w_{h}^{-}-\widetilde{w}_{h}\right\|_{L^{2}(\Gamma)}^{2}$
$\displaystyle-2\left\|\partial_{n}w_{h}^{-}\right\|_{L^{2}(\Gamma)}\left\|w_{h}^{-}-\widetilde{w}_{h}\right\|_{L^{2}(\Gamma)}$
$\displaystyle\geq$
$\displaystyle|w_{h}^{-}|^{2}_{H^{1}(\Omega^{-})}+\varepsilon_{\min}\|w_{h}^{-}\|_{\Omega^{-}}^{2}+\tfrac{\tau}{h}\left\|w_{h}^{-}-\widetilde{w}_{h}\right\|_{L^{2}(\Gamma)}^{2}$
$\displaystyle-2\left|w_{h}^{-}\right|_{H^{1}(\Omega^{-})}\left(C_{max}h^{-\frac{1}{2}}\left\|w_{h}^{-}-\widetilde{w}_{h}\right\|_{L^{2}(\Gamma)}\right)$
$\displaystyle\geq$
$\displaystyle\tfrac{1}{2}|w_{h}^{-}|^{2}_{H^{1}(\Omega^{-})}+\varepsilon_{\min}\|w_{h}^{-}\|_{\Omega^{-}}^{2}+\tfrac{\tau-2C_{max}^{2}}{h}\left\|w_{h}^{-}-\widetilde{w}_{h}\right\|_{L^{2}(\Gamma)}^{2}.$

We finish by applying the equivalence of the norms (3.20) under the assumption
that $\tau>2C_{max}^{2}$.

In the case of bilinear form $b_{h}$, by using the results from Lemma 4 we
obtain, with
$\bar{w}_{h}=|\Gamma|^{-1}\int_{\Gamma}w^{+}_{h}\leavevmode\nobreak\
\mbox{d}s$,

$\displaystyle
b_{h}\left(\left(w_{h}^{+},\lambda_{h},{\widetilde{w}_{h}}\right),\left({w_{h}}^{+},\lambda_{h},{\widetilde{w}_{h}}\right)\right)\geq$
$\displaystyle\alpha_{V}\left\|\lambda_{h}\right\|_{H^{-\frac{1}{2}}(\Gamma)}^{2}+\alpha_{W}\left\|w_{h}^{+}-\bar{w}_{h}\right\|_{H^{\frac{1}{2}}(\Gamma)}^{2}$
$\displaystyle+\tfrac{\tau}{h}\left\|w_{h}^{+}-\widetilde{w}_{h}\right\|_{L^{2}(\Gamma)}^{2}.$

Observe that when $\varepsilon>0$ we may bound

$\displaystyle\|w_{h}^{+}\|_{H^{\frac{1}{2}}(\Gamma)}$
$\displaystyle\leq\|w_{h}^{+}-\bar{w}_{h}\|_{H^{\frac{1}{2}}(\Gamma)}+\|\bar{w}_{h}\|_{L^{2}(\Gamma)}$
$\displaystyle\leq\|w_{h}^{+}-\bar{w}_{h}\|_{H^{\frac{1}{2}}(\Gamma)}+\|w_{h}^{+}-\widetilde{w}_{h}\|_{L^{2}(\Gamma)}+\|w_{h}^{-}-\widetilde{w}_{h}\|_{L^{2}(\Gamma)}+\|w_{h}^{-}\|_{L^{2}(\Gamma)}$
$\displaystyle\leq\|w_{h}^{+}-\bar{w}_{h}\|_{H^{\frac{1}{2}}(\Gamma)}+\|w_{h}^{+}-\widetilde{w}_{h}\|_{L^{2}(\Gamma)}+\|w_{h}^{-}-\widetilde{w}_{h}\|_{L^{2}(\Gamma)}+C\|w_{h}^{-}\|_{H^{1}(\Omega^{-})}$

where we applied the trace inequality (2.4) in the last step. The right hand
side is controlled by the lower bounds on $a_{h}$ and $b_{h}$ above. Once
again, we finish by applying the equivalence of the norms (3.21).

In case $\varepsilon=0$ we need to show that a Poincaré inequality holds,
similar to (2.13), this time on the form

(3.25)
$c_{P}\|w_{h}^{-}\|_{H^{1}(\Omega^{-})}\leq|w_{h}^{-}|_{H^{1}(\Omega^{-})}+\tfrac{1}{h}\left\|w_{h}^{-}-\widetilde{w}_{h}\right\|_{L^{2}(\Gamma)}+\tfrac{1}{h}\left\|w_{h}^{+}-\widetilde{w}_{h}\right\|_{L^{2}(\Gamma)}\\\
+\left\|\lambda_{h}\right\|_{H^{-\frac{1}{2}}(\Gamma)}+\left\|w_{h}^{+}-\bar{w}_{h}\right\|_{H^{\frac{1}{2}}(\Gamma)}$

To this end, since coercivity holds up to a constant, we may assume that
$w^{-}_{h}=\tilde{w}_{h}=w^{+}_{h}=\bar{w}_{h}$ and proceed verbatim as in the
continuous case, since in that case the continuous and discrete expressions
corresponding to (2.14) are the same. ∎

The existence and uniqueness of the solution of problem (3.16) is achieved by
using the Lax-Milgram theorem. In addition, the proposed method is consistent
as the following result shows.

######  Lemma 9 (Consistency).

Let $\delta>0$, $u^{-}\in H^{\frac{3}{2}+\delta}(\Omega^{-})$, $u^{+}\in
H^{\frac{1}{2}}(\Gamma)$ and $\partial_{n}u^{+}=\lambda\in L^{2}(\Gamma)$ be
the solution of problem (1.1) and $\widetilde{u}=u^{-}=u^{+}$ on $\Gamma$. If
$\left(u_{h},\lambda_{h},\widetilde{u}_{h}\right)\in\mathcal{V}_{h}\times\Lambda_{h}^{l}\times
M_{h}^{m}$ solves (3.16) then, for all
$\left(v_{h},\varphi_{h},\widetilde{v}_{h}\right)\in\mathcal{V}_{h}\times\Lambda_{h}^{l}\times
M_{h}^{m}$ the following holds

$a_{h}\left(\left(u^{-}-u_{h}^{-},\widetilde{u}-\widetilde{u}_{h}\right),\left(v_{h}^{-},\widetilde{v}_{h}\right)\right)+b_{h}\left(\left(u^{+}-u_{h}^{+},\lambda-\lambda_{h},\widetilde{u}-\widetilde{u}_{h}\right),\left(v_{h}^{+},\varphi_{h},\widetilde{v}_{h}\right)\right)=0.$

###### Proof.

Because of the transmission conditions from (1.1) we have $u=u^{+}=u^{-}$ and
$\partial_{n}u=\partial_{n}u^{+}=\partial_{n}u^{-}$ on $\Gamma$

$\displaystyle
a_{h}\left(\left(u-u_{h}^{-},\widetilde{u}-\widetilde{u}_{h}\right),\left(v_{h}^{-},\widetilde{v}_{h}\right)\right)=$
$\displaystyle\left<\partial_{n}u,\widetilde{v}_{h}\right>_{\Gamma}-\left<u-\widetilde{u},\partial_{n}v_{h}^{-}\right>_{\Gamma}$
$\displaystyle+\tfrac{\tau}{h}\left<u-\widetilde{u},v_{h}^{-}-\widetilde{v}_{h}\right>_{\Gamma},$
$\displaystyle
b_{h}\left(\left(u-u_{h}^{+},\lambda-\lambda_{h},\widetilde{u}-\widetilde{u}_{h}\right),\left(v_{h}^{+},\varphi_{h},\widetilde{v}_{h}\right)\right)=$
$\displaystyle-\left<\lambda,\widetilde{v}_{h}\right>_{\Gamma}-\left<u-\widetilde{u},\varphi\right>_{\Gamma}$
$\displaystyle+\tfrac{\tau}{h}\left<u-\widetilde{u},v_{h}^{+}-\widetilde{v}_{h}\right>_{\Gamma}.$

By adding above expressions and using the facts that $\lambda=\partial_{n}u$
and $\widetilde{u}=u$ on $\Gamma$, we obtain consistency. ∎

###  3.3 Error analysis

In this section we present the error estimates for the method. These estimates
are proved using the following norm

(3.26)
$\big{\|}\big{(}u,\lambda,\widetilde{u}\big{)}\big{\|}_{h}:=\left\|\left(u^{-},\widetilde{u}\right)\right\|_{\mathcal{F}}+\left\|\left(u^{+},\lambda,\widetilde{u}\right)\right\|_{\mathcal{B}}.$

The first step is the following version of Cea’s lemma.

######  Lemma 10 (Cea’s Lemma).

Let $\delta>0$, $u^{-}\in H^{\frac{3}{2}+\delta}(\Omega^{-})$, $u^{+}\in
H^{\frac{1}{2}}(\Gamma)$ and $\partial_{n}u^{+}=\lambda\in L^{2}(\Gamma)$ be
the solution of problem (1.1), $\widetilde{u}=u^{-}=u^{+}$ on $\Gamma$, and
let
$\left(u_{h},\lambda_{h},\widetilde{u}_{h}\right)\in\mathcal{V}_{h}\times\Lambda_{h}^{l}\times
M_{h}^{m}$ solve (3.16). Then there exists $C>0$, independent of $h$, such
that

(3.27)
$\left\|\left(u-u_{h},\lambda-\lambda_{h},\widetilde{u}-{\widetilde{u}_{h}}\right)\right\|_{h}\leq
C\inf_{\left(v_{h},\varphi_{h},\widetilde{v}_{h}\right)\in\mathcal{V}_{h}\times\Lambda_{h}^{l}\times
M_{h}^{m}}\left\|\left(u-v_{h},\lambda-\varphi_{h},\widetilde{u}-{\widetilde{v}_{h}}\right)\right\|_{h}.$

###### Proof.

Let us denote

$\displaystyle
A_{h}\left(\left(w_{h},\phi_{h},\widetilde{w}_{h}\right),\left(v_{h},\varphi_{h},\widetilde{v}_{h}\right)\right):=$
$\displaystyle
a\left(\left(w_{h}^{-},\widetilde{w}_{h}\right),\left(v_{h}^{-},\widetilde{v}_{h}\right)\right)$
$\displaystyle+b\left(\left(w_{h}^{+},\phi_{h},\widetilde{w}_{h}\right),\left(v_{h}^{+},\varphi_{h},\widetilde{v}_{h}\right)\right).$

Using Lemma 8, we get that there exists $\alpha>0$, independent of $h$, such
that for all
$\left(v_{h},\varphi_{h},\widetilde{v}_{h}\right)\in\mathcal{V}_{h}\times\Lambda_{h}^{l}\times
M_{h}^{m}$ there exists
$\left(w_{h},\phi_{h},\widetilde{w}_{h}\right)\in\mathcal{V}_{h}\times\Lambda_{h}^{l}\times
M_{h}^{m}$ such that
$\left\|\left(w_{h},\phi_{h},\widetilde{w}_{h}\right)\right\|_{h}=1$, and

(3.28)
$A_{h}\left(\left(v_{h},\varphi_{h},\widetilde{v}_{h}\right),\left(w_{h},\phi_{h},\widetilde{w}_{h}\right)\right)\geq\alpha\left\|\left(v_{h},\varphi_{h},\widetilde{v}_{h}\right)\right\|_{h}.$

Now using Lemma 7, we get continuity of $A_{h}$, there exists $\beta>0$

(3.29)
$\displaystyle\left|A_{h}\left(\left(v,\varphi,\widetilde{v}\right)\right),\left(w,\phi,\widetilde{w}\right)\right|$
$\displaystyle\leq\beta\left\|\left(v,\varphi,\widetilde{v}\right)\right\|_{h}\left\|\left(w,\phi,\widetilde{w}\right)\right\|_{h}.$

Let
$\left(v_{h},\varphi_{h},\widetilde{v}_{h}\right)\in\mathcal{V}_{h}\times\Lambda_{h}^{l}\times
M_{h}^{m}$. Then, using the triangle inequality we see that

$\displaystyle\left\|\left(u-u_{h},\lambda-\lambda_{h},\widetilde{u}-{\widetilde{u}_{h}}\right)\right\|_{h}\leq$
$\displaystyle\left\|\left(u-v_{h},\lambda-\varphi_{h},\widetilde{u}-{\widetilde{v}_{h}}\right)\right\|_{h}$
$\displaystyle+\left\|\left(v_{h}-u_{h},\varphi_{h}-\lambda_{h},\widetilde{v}_{h}-\widetilde{u}_{h}\right)\right\|_{h}.$

and it follows from (3.28), Lemma 9 and (3.29) that

$\displaystyle\left\|\left(v_{h}-u_{h},\varphi_{h}-\lambda_{h},\widetilde{v}_{h}-\widetilde{u}_{h}\right)\right\|_{h}\leq$
$\displaystyle\tfrac{1}{\alpha}A_{h}\left(\left(v_{h}-u,\varphi_{h}-\lambda,{\widetilde{v}_{h}}-\widetilde{u}\right),\left(w_{h},\phi_{h},\widetilde{w}_{h}\right)\right)$
$\displaystyle+\tfrac{1}{\alpha}A_{h}\left(\left(u-u_{h},\lambda-\lambda_{h},\widetilde{u}-\widetilde{u}_{h}\right),\left(w_{h},\phi_{h},\widetilde{w}_{h}\right)\right)$
$\displaystyle\leq$
$\displaystyle\tfrac{\beta}{\alpha}\left\|\left(v_{h}-u,\varphi_{h}-\lambda,{\widetilde{v}_{h}}-\widetilde{u}\right)\right\|_{h}.$

Thus, we get (3.27) with $C:=1+\tfrac{\beta}{\alpha}$. ∎

######  Lemma 11 (Energy norm estimates).

For $s>\tfrac{3}{2}$, $r>1$ and $p\geq\tfrac{1}{2}$, let $u^{-}\in
H^{s}(\Omega^{-})$, $u^{+}\in\widetilde{H}^{r}(\Gamma)$ and
$\partial_{n}u^{+}=\lambda\in\widetilde{H}^{p}(\Gamma)$ be the solution of
problem (1.1). On $\Gamma$ there holds $u^{-}=u^{+}=\widetilde{u}$ on
$\Gamma$. Let
$\left(u_{h},\lambda_{h},\widetilde{u}_{h}\right)\in\mathcal{V}_{h}\times\Lambda_{h}^{l}\times
M_{h}^{m}$ solve (3.16). If the mesh is quasi uniform, then there exists
$C>0$, independent of $h$, such that

(3.30)
$\left\|\left(u-u_{h},\lambda-\lambda_{h},\widetilde{u}-{\widetilde{u}_{h}}\right)\right\|_{h}\leq\\\
C\left(h^{\sigma-1}\|u\|_{H^{\sigma}(\Omega^{-})}+h^{\phi-\frac{1}{2}}\|u\|_{\widetilde{H}^{\phi}(\Gamma)}+h^{\xi-\frac{1}{2}}\|u\|_{\widetilde{H}^{\xi}(\Gamma)}+h^{\psi+\frac{1}{2}}\|\lambda\|_{\widetilde{H}^{\psi}(\Gamma)}\right),$

where $\sigma=\min\\{j+1,s\\}$, $\phi=\min\\{k+1,r\\}$, $\xi=\min\\{m+1,r\\}$
and $\psi=\min\\{l+1,p\\}$.

###### Proof.

The result is a consequence of (3.27) and approximation. Applying triangle
inequality and trace inequality [BS08, Theorem 1.6.6] followed by Young’s
inequality, we obtain

$\displaystyle\left\|\left(u-v_{h},\widetilde{u}-\widetilde{v}_{h}\right)\right\|_{\mathcal{F}}^{2}:=$
$\displaystyle\left\|u-v_{h}\right\|_{H^{1}(\Omega^{-})}^{2}+\tfrac{\tau}{h}\left\|u-v_{h}-\left(\widetilde{u}-\widetilde{v}_{h}\right)\right\|_{L^{2}(\Gamma)}^{2}$
$\displaystyle+h\left\|\partial_{n}u-\partial_{n}v_{h}\right\|_{L^{2}(\Gamma)}^{2}$
$\displaystyle\leq$
$\displaystyle\left\|u-v_{h}\right\|_{H^{1}(\Omega^{-})}^{2}+C_{1}\left(\tfrac{\tau}{h^{2}}\left\|u-v_{h}\right\|_{L^{2}(\Omega^{-})}^{2}+\tau\sum_{K\in\mathcal{T}_{h}}\left\|u-v_{h}\right\|_{H^{1}(K)}^{2}\right)$
$\displaystyle+C_{2}\left(\sum_{K\in\mathcal{T}_{h}}\left\|u-v_{h}\right\|_{H^{1}(K)}^{2}+h^{2}\sum_{K\in\mathcal{T}_{h}}\left\|u-v_{h}\right\|_{H^{2}(K)}^{2}\right)+\tfrac{\tau}{h}\left\|\widetilde{u}-\widetilde{v}_{h}\right\|_{L^{2}(\Gamma)}^{2}.$

Using approximation results [EG04, Theorem 1.109] and [SS11, Theorem 4.3.19]
for the last term, we claim

$\inf_{\left(v_{h},\widetilde{v}_{h}\right)\in\mathcal{V}_{h}\times
M_{h}^{m}}\left\|\left(u-v_{h},\widetilde{u}-\widetilde{v}_{h}\right)\right\|_{\mathcal{F}}\leq
C_{\mathcal{F}}\left(h^{\sigma-1}\|u\|_{H^{\sigma}(\Omega^{-})}+h^{\xi-\frac{1}{2}}\|u\|_{\widetilde{H}^{\xi}(\Gamma)}\right).$

For the boundary part, by applying triangle inequality, we obtain

$\displaystyle\left\|\left(u-v_{h},\lambda-\lambda_{h},\widetilde{u}-\widetilde{v}_{h}\right)\right\|_{\mathcal{B}}^{2}:=$
$\displaystyle\|u-v_{h}\|_{H^{\frac{1}{2}}(\Gamma)}^{2}+\|\lambda-\lambda_{h}\|_{H^{-\frac{1}{2}}(\Gamma)}^{2}+h\left\|\lambda-\lambda_{h}\right\|_{L^{2}(\Gamma)}^{2}$
$\displaystyle+\tfrac{\tau}{h}\left\|u-v_{h}-\left(\widetilde{u}-\widetilde{v}_{h}\right)\right\|_{L^{2}(\Gamma)}^{2}$
$\displaystyle\leq$
$\displaystyle\|u-v_{h}\|_{H^{\frac{1}{2}}(\Gamma)}^{2}+\|\lambda-\lambda_{h}\|_{H^{-\frac{1}{2}}(\Gamma)}^{2}+h\left\|\lambda-\lambda_{h}\right\|_{L^{2}(\Gamma)}^{2}$
$\displaystyle+\tfrac{\tau}{h}\left\|u-v_{h}\right\|_{L^{2}(\Gamma)}^{2}+\tfrac{\tau}{h}\left\|\widetilde{u}-\widetilde{v}_{h}\right\|_{L^{2}(\Gamma)}^{2}.$

Using approximation results [SS11, Theorems 4.3.19, 4.3.20 and 4.3.22], we
claim

$\inf_{\left(v_{h},\lambda_{h},\widetilde{v}_{h}\right)\in\mathcal{V}_{h}\times\Lambda_{h}^{l}\times
M_{h}^{m}}\left\|\left(u-v_{h},\lambda-\lambda_{h},\widetilde{u}-\widetilde{v}_{h}\right)\right\|_{\mathcal{B}}\leq\\\
C_{\mathcal{B}}\left(h^{\phi-\frac{1}{2}}\|u\|_{\widetilde{H}^{\phi}(\Gamma)}+h^{\xi-\frac{1}{2}}\|u\|_{\widetilde{H}^{\xi}(\Gamma)}+h^{\psi+\frac{1}{2}}\|\lambda\|_{\widetilde{H}^{\psi}(\Gamma)}\right).$

We conclude the proof by applying Lemma 10. ∎

##  4 Iterative solution

For the solution of the linear system we will iterate on the Schur complement
for the trace variable, solving independently in the two sub domains. To
justify this split approach we here show that a simple relaxed Jacobi
iteration on the two systems will converge. The condition number of the Schur
complement can be analysed using the arguments of [BEH+19, Section 4].

  1. 1.

Given $\widetilde{u}^{n}$ solve for $u^{n+1}$ and $\lambda^{n+1}$ by solving
the linear system

$A_{h}[(u^{n+1},\lambda^{n+1},\widetilde{u}^{n}),(v,\varphi,0)]=\int_{\Omega^{-}}fv\leavevmode\nobreak\
dx.$

  2. 2.

Given $u^{n+1}$ and $\lambda^{n+1}$, solve for the new trace variable
$\widetilde{u}^{n+1}$, for $\sigma>0$,

$A_{h}[(u^{n+1},\lambda^{n+1},\widetilde{u}^{n+1}),(0,0,\widetilde{v})]+\sigma\tfrac{\tau}{h}\left<\widetilde{u}^{n+1}-\widetilde{u}^{n},\widetilde{v}\right>_{\Gamma}=0.$

To prove that the iterative algorithm converges we only need to show that if
$f=0$, $u^{n+1}$, $\widetilde{u}^{n+1}$ and $\lambda^{n+1}$ all go to zero as
$n\rightarrow\infty$.

We add and subtract $\widetilde{u}^{n+1}$ in the first equation and add the
second to obtain

$A_{h}[(u^{n+1},\lambda^{n+1},\widetilde{u}^{n+1}),(v,\varphi,\widetilde{v})]+\sigma\tfrac{\tau}{h}\left<\widetilde{u}^{n+1}-\widetilde{u}^{n},\widetilde{v}\right>_{\Gamma}=\\\
\left<(\widetilde{u}^{n+1}-\widetilde{u}^{n}),\partial_{n}v+\varphi-\tfrac{\tau}{h}((v^{-}+v^{+})\right>_{\Gamma}.$

Test this equation with $u^{n+1},\ \lambda^{n+1},\ \widetilde{u}^{n+1}$ and
use coercivity to obtain

$\tfrac{1}{2}\sigma\tfrac{\tau}{h}\|\widetilde{u}^{N}\|_{L^{2}(\Gamma)}^{2}+\sum_{n=0}^{N-1}(\alpha\|(u^{n+1},\lambda^{n+1},\widetilde{u}^{n+1})\|_{h}^{2}+\tfrac{1}{2}\sigma\tfrac{\tau}{h}\|\widetilde{u}^{n+1}-\widetilde{u}^{n}\|_{L^{2}(\Gamma)}^{2})\\\
\leq\tfrac{1}{2}\sigma\tfrac{\tau}{h}\|\widetilde{u}^{0}\|_{L^{2}(\Gamma)}^{2}+\sum_{n=0}^{N-1}\left<(\widetilde{u}^{n+1}-\widetilde{u}^{n}),\partial_{n}u^{n+1}+\lambda^{n+1}-\tfrac{\tau}{h}((u^{n+1})^{-}+(u^{n+1})^{+})\right>_{\Gamma}.$

Here we used the well known formula

(4.31)
$\sum_{n=0}^{N-1}\left<(\widetilde{u}^{n+1}-\widetilde{u}^{n}),\widetilde{u}^{n+1}\right>_{\Gamma}=\tfrac{1}{2}\|\widetilde{u}^{N}\|_{L^{2}(\Gamma)}^{2}-\tfrac{1}{2}\|\widetilde{u}^{0}\|_{L^{2}(\Gamma)}^{2}+\tfrac{1}{2}\sum_{n=0}^{N-1}\|\widetilde{u}^{n+1}-\widetilde{u}^{n}\|_{L^{2}(\Gamma)}^{2}.$

Considering the terms on the right hand side and using trace inequality (3.15)
we see that

$\left<(\widetilde{u}^{n+1}-\widetilde{u}^{n}),\partial_{n}(u^{n+1})^{-}\right>_{\Gamma}\leq\tfrac{1}{4}\sigma\tfrac{\tau}{h}\|\widetilde{u}^{n+1}-\widetilde{u}^{n}\|^{2}_{L^{2}(\Gamma)}+C_{max}\tau^{-1}\sigma^{-1}\|u^{n+1}\|_{H^{1}(\Omega^{-})}^{2}.$

Using the duality pairing between $H^{\frac{1}{2}}$ and $H^{-\frac{1}{2}}$
followed by the global inverse inequality
$\|\widetilde{u}^{n+1}-\widetilde{u}^{n}\|_{H^{\frac{1}{2}}(\Gamma)}\leq
C_{t}h^{-\frac{1}{2}}\|\widetilde{u}^{n+1}-\widetilde{u}^{n}\|_{L^{2}(\Gamma)}$
(see [SS11, Theorem 4.4.3]) and Young’s inequality we have

$\left<(\widetilde{u}^{n+1}-\widetilde{u}^{n}),\lambda^{n+1}\right>_{\Gamma}\leq\tfrac{1}{4}\sigma\tfrac{\tau}{h}\|\widetilde{u}^{n+1}-\widetilde{u}^{n}\|^{2}_{L^{2}(\Gamma)}+C_{t}^{-2}\tau^{-1}\sigma^{-1}\|\lambda^{n+1}\|_{H^{-\frac{1}{2}}(\Gamma)}^{2}.$

Finally

$\displaystyle\left<(\widetilde{u}^{n+1}-\widetilde{u}^{n}),\tfrac{\tau}{h}((u^{n+1})^{-}+(u^{n+1})^{+})\right>_{\Gamma}\leq$
$\displaystyle\
\tfrac{1}{4}\sigma\tfrac{\tau}{h}\|\widetilde{u}^{n+1}-\widetilde{u}^{n}\|^{2}_{L^{2}(\Gamma)}$
$\displaystyle+\sigma^{-1}\big{(}\tfrac{\tau}{h}\|(u^{n+1})^{+}-\widetilde{u}^{n+1}\|^{2}_{L^{2}(\Gamma)}$
$\displaystyle+\tfrac{\tau}{h}\|(u^{n+1})^{-}-\widetilde{u}^{n+1}\|^{2}_{L^{2}(\Gamma)}\big{)}$
$\displaystyle+2\tfrac{\tau}{h}\left<(\widetilde{u}^{n+1}-\widetilde{u}^{n}),\widetilde{u}^{n+1}\right>_{\Gamma}.$

Using the once again the telescoping property (4.31) we see that for

$\sigma>\tau^{-1}\alpha^{-1}\max\left\\{C_{max},C_{t}^{-2},\tau\right\\}+2,$

the right hand sides can all be absorbed in the left hand side to yield

$\sum_{n=0}^{N-1}\left(\widetilde{\alpha}\|(u^{n+1},\lambda^{n+1},\widetilde{u}^{n+1})\|_{h}^{2}+\tfrac{1}{4}\sigma\tfrac{\tau}{h}\|\widetilde{u}^{n+1}-\widetilde{u}^{n}\|_{L^{2}(\Gamma)}^{2}\right)\leq\tfrac{1}{2}(\sigma-2)\tfrac{\tau}{h}\|\widetilde{u}^{0}\|_{L^{2}(\Gamma)}^{2}.$

It follows that as $N\rightarrow\infty$ $u^{n+1}$, $\widetilde{u}^{n+1}$ and
$\lambda^{n+1}$ all go to zero, since the sum of the left hand side has to be
bounded by the constant of the right hand side.

##  5 Numerical experiments

In our experiment tests we consider $\varepsilon=1$ and
$\mathcal{V}_{h}\times\Lambda_{h}^{l}\times M_{h}^{m}$ with $j=k=m=1$. The
value $l$ varies depending on the geometry of domains considered. We let the
trace meshes $\mathcal{G}_{1}$ and $\mathcal{G}_{2}$ coincide with the trace
mesh of $\mathcal{T}_{h}$ on $\Gamma$. For our experiments we use two
numerical softwares: FEniCS [ABH+15] and Bempp [SBA+15]. We use the solution
of interior and exterior Dirichlet boundary value problems to construct a
Schur complement system solving the following equations

$\tfrac{\tau}{h}\left(u_{s}-\widetilde{u}\right)+\tfrac{\tau}{h}\left(u_{m}-\widetilde{u}\right)-\left(\partial_{n}u-\lambda\right)=0\mbox{
on }\Gamma.$

The solution $\widetilde{u}$ on $\Gamma$ of the Schur complement is obtained
using the nested conjugate gradient method (CG) [HS52]. Although one can use
direct solvers to solve the interior and exterior Dirichlet boundary value
problems, we here use iterative solvers to apply preconditioners. The interior
Dirichlet boundary value problem that is a symmetric system associated with
bilinear form $a_{h}$ (3.16) is solved by using FEniCS and CG without and with
algebraic multigrid preconditioner. The discrete exterior problem associated
with bilinear form $b_{h}$ (3.16) is not symmetric, however as we shown in
Section 3.1, we can apply the Steklov-Poincaré operator to the flux variable
and transform the equations into a symmetric system. For clarity of the code,
we here simply used the generalized minimal residual method (GMRES) [SS86]
without or with mass matrix preconditioner to solve in Bempp the external
Dirichlet boundary value problem. The tolerance of the iterative solvers is
chosen to be not greater the $10^{-8}$. A Jupyter notebook demonstrating the
functionality used in this paper will be made available at www.bempp.com.

###  5.1 Choice of parameter $\tau$

Thanks to Lemma 8 we know that the stabilisation parameter $\tau$ in the
discrete problem (3.16) must be large enough to assure the coercivity. We
start with an experiment showing how the value of the parameter $\tau$
influences the convergence and number of iterations. We consider $\Omega^{-}$
as a unit sphere with boundary $\Gamma$. We define

$\displaystyle u^{-}(x,y,z)$
$\displaystyle=\frac{1}{2\pi}\sin\left(\pi(x^{2}+y^{2}+z^{2})\right)+\frac{1}{2\pi}\cos\left(\pi(x^{2}+y^{2}+z^{2})\right)+\frac{2\pi+1}{2\pi},$
$\displaystyle u^{+}(x,y,z)$
$\displaystyle=\frac{1}{\sqrt{x^{2}+y^{2}+z^{2}}}.$

It is easy to check that for the unit sphere domain $\Omega^{-}$ the above
elementary functions are the solution of our problem (1.1).

(a) The error of the interior solution.

(b) The error of the exterior solution.

(c) Iteration taken by CG to solve the preconditioned system.

Figure 1: The dependence of the errors and iteration count on the value of
$\tau$ for $h<2^{-1}$ (solid line with circles), $h<2^{-2}$ (dash-dotted line
with diamonds), and $h<2^{-3.5}$ (dashed line with squares) for the problem on
the unit sphere subdomain, with $j=k=m=l=1$.

Figure 1 shows the error values for different values of $\tau$ and
$j=k=m=l=1$. In this case, $\Gamma$ is smooth, and so
$W_{h}^{1}=\Lambda_{h}^{1}$. In Figure 1a, we plot in log-log scale the error
of the interior solution $\|u^{-}-u_{h}^{-}\|_{L^{2}(\Omega^{-})}$ and on
Figure 1b the error of the exterior solutions
$\|u^{+}-u_{h}^{+}\|_{L^{2}(\Gamma)}+\|\lambda^{+}-\lambda_{h}^{+}\|_{L^{2}(\Gamma)}$
for $h<2^{-1}$ (solid line with circles), $h<2^{-2}$ (dash-dotted line with
diamonds), and $h<2^{-3.5}$ (dashed line with squares).

It can be seen from the Figures 1 that both errors stop decreasing when $\tau$
is around 10. Furthermore, for $\tau>10$ the iterations increase with growing
$\tau$, hence we fix $\tau=10$ for the next experiments.

###  5.2 Spherical subdomain

Let $\Omega^{-}$ once again be the unit sphere, $\Gamma$ its boundary and
consider the same exact solution as above.

(a) Error of the interior (dashed line with circles) and exterior solutions
(solid line with squares). The dotted line shows order 2 convergence.

(b) Iteration taken by CG to solve the non-preconditioned system (dashed
line), compared with the preconditioned system (solid line)

(c) Time to solve the non-preconditioned system (dashed line), compared with
the preconditioned system (solid line) of the whole discrete problem (3.16)
including solving interior and exterior systems.

Figure 2: The convergence (left), CG iteration counts (middle) and solving
time (right) for the problem on the unit sphere with $\tau=10$ and
$j=k=m=l=1$.

Figure 2 shows the convergence, CG iteration counts and solving time when
$\tau=10$ and $k=l=1$. In Figure 2a, we plot in log-log scale the error of the
interior solution $\|u^{-}-u_{h}^{-}\|_{L^{2}(\Omega^{-})}$ (dashed line with
circles) and the error of the exterior solutions
$\|u^{+}-u_{h}^{+}\|_{L^{2}(\Gamma)}+\|\lambda^{+}-\lambda_{h}^{+}\|_{L^{2}(\Gamma)}$
(solid line with squares).

In Figure 2b, we plot in log-log scale the number of iterations taken by CG to
solve the non-preconditioned system associated with exterior problem (dashed
line), compared with the preconditioned system (solid line). In addition,
Figure 2c shows the time required by solvers of interior and exterior systems.
The interior system is solved by CG with or without algebraic multigrid
preconditioner and the exterior system is solved by GMRES with or without mass
preconditioner. Preconditioning reduces both the iteration count and the CPU
time needed by the solver.

###  5.3 Cubical subdomain

Let $\Omega^{-}=(0,1)^{3}$ be a cube and we solve the problem (1.1) with
$f=1$. We choose $\tau=10$ and $j=k=m=l+1=1$, where $\Lambda_{h}^{0}$ is the
space of piece-wise constants per element in the trace space.

(a) Error between interior solution and exterior. The dashed line shows order
2 convergence.

(b) Iteration taken by CG to solve the non-preconditioned system (dashed
line), compared with the preconditioned system (solid line)

(c) Time to solve the non-preconditioned system (dashed line), compared with
the preconditioned system (solid line) of the whole discrete problem (3.16)
including solving interior and exterior systems.

Figure 3: The convergence (left), CG iteration counts (middle) and solving
time (right) for the problem on the cube with $\tau=10$ and $j=k=m=l+1=1$.

Figure 3a shows the convergence when $\tau=10$ and $j=k=m=l+1=1$. In this
case, the exact solution is not known, thus in Figure 3a, we plot in log-log
scale the error between interior solution and exterior
$\frac{\|u_{h}^{-}|_{\Gamma}-u_{h}^{+}\|_{L^{2}(\Gamma)}}{\|u_{h}^{+}\|_{L^{2}(\Gamma)}}$
(solid line).

In Figure 3, we plot as well in log-log scale the number of iterations and
solving time taken by CG to solve the non-preconditioned system (dashed line),
compared with the preconditioned system (solid line). Once again
preconditioning brings improvement in terms of iteration counts and time taken
to solve the problem.

##  6 Conclusions

We have analyzed and demonstrated the effectiveness of Nitsche type methods
for coupling finite element and boundary element formulations. Our approach
gives flexibility to choose a continuous or discontinuous finite element space
in the FEM solver, hence the interior problem can be solved essentially using
any method that allows for the hybridised Nitsche method for interdomain
coupling. We are also free to choose the trace variable minimising the
coupling degrees of freedom. In this paper, we focus on the technical aspects
and analysis to allow for flexibility within our framework. A work
demonstrating the applicability for large problems using parallel approach is
in preparation.

The method can be extended to other models such as the Helmholtz equations. In
this case it is known [SW11] that the use of impedance interface conditions is
advantageous and such an approach can be mimicked in the present framework by
letting the stabilisation constant have non-zero imaginary part and depend on
the wave number. Formulations of the presented FEM/BEM coupling method to the
Helmholtz and Maxwell problems are currently in preparation. For these cases
however more effective operator preconditioning techniques for exterior
problem are essential, especially for high frequency problems. Despite that,
we expect that their implementation will be similar to the presented Laplace
case.

## References

  * [ABH+15] M. S. Alnæs, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg, C. Richardson, J. Ring, M. E. Rognes, and G. N. Wells.  The fenics project version 1.5.  Archive of Numerical Software, 3(100), 2015.
  * [BBS19] T. Betcke, E. Burman, and M. W. Scroggs.  Boundary element methods with weakly imposed boundary conditions.  SIAM J. Sci. Comput., 41(3):A1357–A1384, 2019.
  * [BCS96] U. Brink, C. Carstensen, and E. Stein.  Symmetric coupling of boundary elements and Raviart-Thomas-type mixed finite elements in elastostatics.  Numer. Math., 75(2):153–174, 1996.
  * [BEH+19] E. Burman, D. Elfverson, P. Hansbo, M. G. Larson, and K. Larsson.  Hybridized CutFEM for elliptic interface problems.  SIAM J. Sci. Comput., 41(5):A3354–A3380, 2019.
  * [BHL21] E. Burman, P. Hansbo, and M. G. Larson.  CutFEM Based on Extended Finite Element Spaces.  arXiv e-prints, page arXiv:2101.10052, January 2021.
  * [BJ79] F. Brezzi and C. Johnson.  On the coupling of boundary integral and finite element methods.  Calcolo, 16(2):189–201, 1979.
  * [BJN78] F. Brezzi, C. Johnson, and J.-C. Nédélec.  On the coupling of boundary integral and finite element methods.  In Proceedings of the Fourth Symposium on Basic Problems of Numerical Mathematics (Plzeň, 1978), pages 103–114. Charles Univ., Prague, 1978.
  * [BS08] 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.
  * [CF00] C. Carstensen and S. A. Funken.  Coupling of mixed finite elements and boundary elements.  IMA J. Numer. Anal., 20(3):461–480, 2000.
  * [CGL09] B. Cockburn, J. Gopalakrishnan, and R. Lazarov.  Unified hybridization of discontinuous Galerkin, mixed, and continuous Galerkin methods for second order elliptic problems.  SIAM J. Numer. Anal., 47(2):1319–1365, 2009.
  * [CGS12] B. Cockburn, J. Guzmán, and F.-J. Sayas.  Coupling of Raviart-Thomas and hybridizable discontinuous Galerkin methods with BEM.  SIAM J. Numer. Anal., 50(5):2778–2801, 2012.
  * [CH12] F. Chouly and N. Heuer.  A Nitsche-based domain decomposition method for hypersingular integral equations.  Numer. Math., 121(4):705–729, 2012.
  * [Cos87] 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.
  * [CS12] B. Cockburn and F.-J. Sayas.  The devising of symmetric couplings of boundary element and discontinuous Galerkin methods.  IMA J. Numer. Anal., 32(3):765–794, 2012.
  * [DPE12] D. A. Di Pietro and A. Ern.  Mathematical aspects of discontinuous Galerkin methods, volume 69 of Mathématiques & Applications (Berlin) [Mathematics & Applications].  Springer, Heidelberg, 2012.
  * [DPEL14] D. A. Di Pietro, A. Ern, and S. Lemaire.  An arbitrary-order and compact-stencil discretization of diffusion on general meshes based on local reconstruction operators.  Comput. Methods Appl. Math., 14(4):461–472, 2014.
  * [EG04] A. Ern and J. L. Guermond.  Theory and practice of finite elements, volume 159 of Applied Mathematical Sciences.  Springer-Verlag, New York, 2004.
  * [Egg09] H. Egger.  A class of hybrid mortar finite element methods for interface problems with non-matching meshes, 2009.
  * [ES10] H. Egger and J. Schöberl.  A hybrid mixed discontinuous Galerkin finite-element method for convection-diffusion problems.  IMA J. Numer. Anal., 30(4):1206–1234, 2010.
  * [GHH09] G. N. Gatica, M. Healey, and N. Heuer.  The boundary element method with Lagrangian multipliers.  Numer. Methods Partial Differential Equations, 25(6):1303–1319, 2009.
  * [GHS10] G. N. Gatica, N. Heuer, and F.-J. Sayas.  A direct coupling of local discontinuous Galerkin and boundary element methods.  Math. Comp., 79(271):1369–1394, 2010.
  * [GS06] G. N. Gatica and F.-J. Sayas.  An a priori error analysis for the coupling of local discontinuous Galerkin and boundary element methods.  Math. Comp., 75(256):1675–1696, 2006.
  * [Han90] 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.
  * [HS52] M. R. Hestenes and E. Stiefel.  Methods of conjugate gradients for solving linear systems.  J. Research Nat. Bur. Standards, 49:409–436 (1953), 1952.
  * [JN80] 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.
  * [LS05] U. Langer and O. Steinbach.  Coupled boundary and finite element tearing and interconnecting methods.  In Domain decomposition methods in science and engineering, volume 40 of Lect. Notes Comput. Sci. Eng., pages 83–97. Springer, Berlin, 2005.
  * [MVMP96] S. Meddahi, J. Valdés, O. Menéndez, and P. Pérez.  On the coupling of boundary integral and mixed finite element methods.  J. Comput. Appl. Math., 69(1):113–124, 1996.
  * [Nit71] J. Nitsche.  Über ein Variationsprinzip zur Lösung von Dirichlet-Problemen bei Verwendung von Teilräumen, die keinen Randbedingungen unterworfen sind.  Abh. Math. Sem. Univ. Hamburg, 36:9–15, 1971.
  * [ORST12] G. Of, G. J. Rodin, O. Steinbach, and M. Taus.  Coupling of discontinuous Galerkin finite element and boundary element methods.  SIAM J. Sci. Comput., 34(3):A1659–A1677, 2012.
  * [Say09] 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.
  * [SBA+15] W. Śmigaj, T. Betcke, S. Arridge, J. Phillips, and M. Schweiger.  Solving boundary integral problems with BEM++.  ACM Transactions on Mathematical Software, 41(2):6:1–6:40, 2015.
  * [SS86] Y. Saad and M. H. Schultz.  GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems.  SIAM J. Sci. Statist. Comput., 7(3):856–869, 1986.
  * [SS11] S. A. Sauter and C. Schwab.  Boundary element methods, volume 39 of Springer Series in Computational Mathematics.  Springer-Verlag, Berlin, 2011.  Translated and expanded from the 2004 German original.
  * [Ste08] O. Steinbach.  Numerical approximation methods for elliptic boundary value problems.  Springer, New York, 2008.  Finite and boundary elements, Translated from the 2003 German original.
  * [Ste11] O. Steinbach.  A note on the stable one-equation coupling of finite and boundary elements.  SIAM J. Numer. Anal., 49(4):1521–1531, 2011.
  * [SW11] O. Steinbach and M. Windisch.  Stable boundary element domain decomposition methods for the Helmholtz equation.  Numer. Math., 118(1):171–195, 2011.
  * [Wen86] W. L. Wendland.  On asymptotic error estimates for the combined boundary and finite element method.  In Innovative numerical methods in engineering (Atlanta, Ga., 1986), pages 55–69. Comput. Mech., Southampton, 1986.
  * [ZKB77] O. C. Zienkiewicz, D. W. Kelly, and P. Bettess.  The coupling of the finite element method and boundary solution procedures.  Internat. J. Numer. Methods Engrg., 11(2):355–375, 1977.