We consider boundary element methods where the Calderón projector is
used for the system matrix and boundary conditions are weakly imposed using a
particular variational boundary operator designed using techniques
from augmented Lagrangian methods. Regardless of the boundary conditions,
both the primal trace variable and the flux are approximated. We focus on the imposition of Dirichlet,
mixed Dirichlet–Neumann, and Robin conditions. A salient feature of the
Robin condition is that the conditioning of the system is robust also
for stiff boundary conditions. The theory is illustrated by a series of numerical examples.

boundary element methods, Nitsche's method, Robin boundary conditions, mixed boundary conditions

65N38, 65R20

§ INTRODUCTION

Weak imposition of boundary conditions has been very successful in the
context of finite element methods. In particular, Nitsche's method [19] has
recently received increased interest in the scientific computation
Our aim in this paper is to discuss how the idea behind this type of
method can be applied in the context of boundary element methods to
impose different types of boundary condition in a unified

Weak imposition of boundary conditions here means that neither the Dirichlet trace nor the Neumann trace is imposed exactly,
instead an $h$-dependent boundary condition is imposed that is weighted in such a way that optimal error estimates may be derived
and the exact boundary condition is recovered in the asymptotic limit.
Methods based on Nitsche's method have been succesfully utilised for boundary element method domain decomposition problems,
where they have been used to impose interface conditions at 1D interfaces between segments of 2D screens embedded in 3D space [13, 10].
Our approach instead focusses on imposing
boundary conditions on the 2D boundary of a single domain problem through the addition of penalty terms to a general formulation written in terms of the multitrace
operator, in a similar vein to the method discussed in [1] for the finite element method.

The use of systems of boundary integral equations for problems with mixed boundary conditions is quite classical
[11, 25, 26, 27]. While these papers require the assembly of boundary operators on subsets of the boundary
mesh, the penalty method proposed in this paper requires only the addition of sparse mass matrices to the multitrace operator assembled on the
entire mesh. In addition to the greater simplicity of the resulting formulation, this method has the advantage that the sparse penalty terms
only affect the entries in the matrix for near interactions: this gives the resulting system a structure that can be utilised when
designing effective preconditioners.

This approach may not be competitive in the simple
case of pure Dirichlet or Neumann conditions due to the increase in the number of unknowns.
Therefore the main focus of this work is on more complex situations.
We will discuss the following four model cases:

* non-homogeneous Dirichlet conditions,

* non-homogeneous Neumann conditions,

* mixed Dirichlet–Neumann boundary conditions,

* generalised Robin conditions.

We consider the Laplace equation: Find $u$ such that

\begin{align}
-\Delta u &=0       &&\text{in } \Omega,\\
u &= g\D            &&\text{on } \Gamma\D, \\
\frac{\partial u}{\partial\B\nu} &= g\N
&&\text{on } \Gamma\N,\\
\frac{\partial u}{\partial\B\nu} &= \frac{1}{\param} (g\D-u) + g\N
&&\text{on }\Gamma\R.\label{eq:LaplaceRobin} %\\
\end{align}

Here $\Omega \subset \RR^3$ denotes a polyhedral domain with outward pointing normal
$\B\nu$ and boundary $\Gamma:= \Gamma\D \cup \Gamma\N \cup \Gamma\R$.
We assume for simplicity that the boundaries between $\Gamma\D$, $\Gamma\N$ and $\Gamma\R$ coincide with edges between the faces of $\Gamma$.
Whenever it is ambiguous, we will write $\B\nu_\x$ for the outward pointing normal at the point $\x$.
We assume that $g\D \in H^{1/2}(\Gamma\D \cup \Gamma\R)$ and $g\N \in
L^2(\Gamma\N \cup \Gamma\R)$. Observe that, by the Lax–Milgram lemma, there exists a unique solution
to <ref>. We assume that $u \in H^{3/2+\epsilon}(\Omega)$ for some $\epsilon>0$.

For the Robin boundary condition, we will use the ideas of Juntunen and
Stenberg [16]. A salient feature of this type of imposition
of the Robin condition is that it is robust under singular
perturbations. Indeed regardless of the Robin coefficient, the
conditioning of the resulting system matrix is no worse than
for the Neumann or the Dirichlet problem.

The proposed framework is flexible and allows for the design of a range of
different methods depending on the choice of weights and residuals. We
will present a sample of possible methods with the ambition of showing
the versatility of the framework rather than claiming
that for each case the choices are optimal.

An outline of the paper is as follows. First, we review some of the
basic elements of the theory of boundary operators in <ref>. Then, in <ref> we discuss
the design of formulations for the linear model problems in a formal setting. We propose the
corresponding boundary element methods in <ref> and give an
abstract analysis. The boundary elements obtained using the
formulations from <ref> are then shown to satisfy the
assumptions of the abstract theory.
Finally, we show some computational examples in <ref>.

While the present paper focuses on weak imposition of boundary conditions through Nitsche type coupling for BEM, ultimately the goal is to
develop a framework for complex BEM/BEM and FEM/BEM multiphysics coupling situations. Existing approaches here are often built upon FETI
and BETI type methods [17, 18]. While BETI is usually formulated in terms of Steklov–Poincaré operators, the
framework proposed in this paper builds directly upon Calderón projectors of the subdomains.

For the method proposed in the present work the multi-domain coupling will take a form similar to that using Nitsche's method in the FEM/FEM
coupling setting of [5]; see also the FEM/BEM coupling of [9] where a Nitsche's method for the coupling was
proposed, using the Steklov-Poincaré operator for the BEM system.

An important application area for the presented weak imposition of boundary conditions are inverse problems with unknown boundary
conditions. Since the boundary condition only enters through a sparse operator this can be easily updated in each step of a solver
iteration, while the boundary integral operators only need to be computed once.
In particular, for reconstruction of the coefficient in a Robin condition (see eg [15] for a finite element approach and [3]
for a detailed analysis of the stability of this problem), the robustness with respect to the coefficient of the present method is an

§ BOUNDARY OPERATORS

We define the Green's function for the Laplace operator in $\RR^3$ by

\begin{equation}
G(\x,\y) =
\frac{1}{4\uppi|\x-\y|}.
\end{equation}

In this paper, we focus on the problem in
$\RR^3$. Similar analysis can be used for problems in $\RR^2$, in which case this definition should be replaced by
$G(\x,\y) =-\log|\x-\y|/2\uppi$.

In the standard fashion (see eg <cit.>), we define
the single layer potential operator, $\pop{V}:H^{-1/2}(\Gamma)\to H^1(\Omega)$,
the double layer potential operator, $\pop{K}:H^{1/2}(\Gamma)\to H^1(\Omega)$,
for $v\in H^{1/2}(\Gamma)$, $\mu\in H^{-1/2}(\Gamma)$, and $\x\in\Omega\setminus\Gamma$ by

\begin{align}
(\popV\mu)(\x) &:= \int_{\Gamma} G(\x,\y) \mu(\y)\dx[\y]\label{eq:single},\\
(\popK v)(\x) &:= \int_{\Gamma} \frac{\partial G(\x,\y)}{\partial\B\nu_\y} v(\y)\dx[\y]\label{eq:double}.
\end{align}

We define the space $H^1(\Delta,\Omega):=\{v\in H^1(\Omega):\Delta v\in\Ltwo(\Omega)\}$, and
then we define the Dirichlet and Neumann traces, $\trace\D:H^1(\Omega)\to H^{1/2}(\Gamma)$
and $\trace\N:H^1(\Delta,\Omega)\to H^{-1/2}(\Gamma)$, by

\begin{align}
\trace\D f(\x)&:=\lim_{\Omega\ni\y\to\x\in\Gamma}f(\y),\\
\trace\N f(\x)&:=\lim_{\Omega\ni\y\to\x\in\Gamma}\B\nu_\x\cdot\nabla f(\y).
\end{align}

We recall that if the Dirichlet and Neumann traces of a
harmonic function are known, then the potentials <ref> and
<ref> may be used to reconstruct the function in $\Omega$ using the
following relation.

\begin{equation}\label{eq:represent}
u = -\popK(\trace\D u) + \popV(\trace\N u).
\end{equation}

It is also known <cit.> that for all $\mu \in H^{-1/2}(\Gamma)$, the function

\begin{equation}\label{eq:sing_rec}
u^{\popV}_\mu := \popV \mu
\end{equation}

satisfies $-\Delta u^{\popV}_\mu = 0$ and

\begin{equation}\label{eq:sing_rec_stab}
\Hnorm[\Omega]{1}{u^{\popV}_\mu} \leqslant c \Hnorm{-1/2}{\mu}.
\end{equation}

Similarly, for the double layer potential there holds <cit.> that for all $v \in H^{1/2}(\Gamma)$, the function

\begin{equation}\label{eq:doub_rec}
u^{\popK}_v := \popK v
\end{equation}

satisfies $-\Delta u^{\popK}_v = 0$ and

\begin{equation}\label{eq:doub_rec_stab}
\Hnorm[\Omega]{1}{u^{\popK}_v} \leqslant c \Hnorm{1/2}{v}.
\end{equation}

We define
$\average{\trace\D f}$ and $\average{\trace\N f}$
to be the averages of the interior and exterior Dirichlet and Neumann traces of $f$.
We define the single layer, double layer, adjoint double layer, and hypersingular boundary integral operators,
$\bopV:H^{-1/2}(\Gamma)\to H^{1/2}(\Gamma)$,
$\bopK:H^{1/2}(\Gamma)\to H^{1/2}(\Gamma)$,
$\bopKadj:H^{-1/2}(\Gamma)\to H^{-1/2}(\Gamma)$, and
$\bopW:H^{1/2}(\Gamma)\to H^{-1/2}(\Gamma)$,

\begin{align}
(\bopK v)(\x)&:=\average{\trace\D\popK v}(\x),&
(\bopV \mu)(\x)&:=\average{\trace\D\popV\mu}(\x),\\
(\bopW v)(\x)&:=-\average{\trace\N\popK v}(\x),&
\end{align}

$v\in H^{1/2}(\Gamma)$ and $\mu\in H^{-1/2}(\Gamma)$

The following coercivity results are known for the single layer and hypersingular operators in
$\RR^3$, where $\left\langle  \cdot, \cdot \right\rangle_\Gamma$ denotes the
$H^{1/2}(\Gamma)$–$H^{-1/2}(\Gamma)$ duality pairing.

There exists $\alpha_\bopV>0$ such that

\begin{align*}
\alpha_\bopV \Hnorm{-1/2}{\mu}^2 &\leqslant \left\langle \bopV \mu,\mu\right\rangle_{\Gamma},
&&\forall \mu \in H^{-1/2}(\Gamma).
\end{align*}

There exists $\alpha_\bopW>0$ such that

\begin{align*}
\alpha_\bopW \Hnorm{1/2}{v}^2 &\leqslant \left\langle\bopW v,v\right\rangle_{\Gamma},
&&\forall v \in H_*^{1/2}(\Gamma),
\end{align*}

where $H_*^{1/2}(\Gamma)$ denotes the set of functions
$v\in H^{1/2}(\Gamma)$ such that $\overline{v} = 0$,
where $\displaystyle\overline{v}:=\frac{\langle v,1\rangle_\Gamma}{\langle1,1\rangle_\Gamma}$ is the average value of $v$.
From this it follows that

\begin{align*}
\alpha_\bopW \Hstarnorm{1/2}{v}^2 &\leqslant \left\langle\bopW v,v\right\rangle_{\Gamma},
&&\forall v \in H^{1/2}(\Gamma),
\end{align*}

where $\Hstarnorm{1/2}{\cdot}$ is defined, for $v\in H^{1/2}(\Gamma)$, by

The following boundedness results are also known.

There exist $C_\bopV,C_\bopK,C_\bopKadj,C_\bopW>0$ such that

\begin{align*}
\text{i)}&&\Hnorm{1/2}{\bopV \mu} &\leqslant C_\bopV \Hnorm{-1/2}{\mu} &&\forall \mu \in H^{-1/2}(\Gamma),\\
\text{ii)}&&\Hnorm{1/2}{\bopK v} &\leqslant C_\bopK \Hnorm{1/2}{v} &&\forall v \in H^{1/2}(\Gamma),\\
\text{iii)}&&\Hnorm{-1/2}{\bopKadj \mu} &\leqslant C_\bopKadj \Hnorm{-1/2}{\mu} &&\forall \mu \in H^{-1/2}(\Gamma),\\
\text{iv)}&&\Hnorm{-1/2}{\bopW v} &\leqslant C_\bopW \Hnorm{1/2}{v} &&\forall v \in H^{1/2}(\Gamma).
\end{align*}

We define the Calderón projector by

\begin{equation}\label{eq:calder}
\bop{C}:=
\begin{pmatrix}
(1-\sigma)\bopI-\bopK & \bopV\\
\bopW & \sigma\bopI+\bopKadj
\end{pmatrix},
\end{equation}

where $\sigma$ is defined as in <cit.>,
and recall that if $u$ is a solution of <ref> then it satisfies

\begin{equation}\label{eq:calder_id}
\bop{C}\begin{pmatrix}\trace\D u\\\trace\N u\end{pmatrix}=\begin{pmatrix}\trace\D u\\\trace\N u\end{pmatrix}.
\end{equation}

Taking the product of <ref> with two test functions, and using the fact that $\sigma=\frac12$ almost everywhere,
we arrive at the following equations.

\begin{align}
\left\langle \trace\D u,\mu \right\rangle_{\Gamma} &= \left\langle (\tfrac12\bopI - \bopK)
\trace\D u,\mu\right\rangle_{\Gamma} + \left\langle \bopV\trace\N
u,\mu\right\rangle_{\Gamma} &&\forall \mu \in H^{-1/2}(\Gamma),\label{eq:Calderon_1}\\
\left\langle \trace\N u,v \right\rangle_{\Gamma} &= \left\langle (\tfrac12\bopI + \bopKadj)
\trace\N u,v \right\rangle_{\Gamma} + \left\langle \bopW \trace\D
u,v\right\rangle_{\Gamma} &&\forall v \in H^{1/2}(\Gamma).\label{eq:Calderon_2}
\end{align}

For a more compact notation, we introduce $\lambda=\trace\N u$ and $u=\trace\D u$ and the Calderón form

\begin{multline}\label{eq:compact_form}
\form{C}[(u,\lambda),(v,\mu)]:=
\left\langle(\tfrac12\bopI-\bopK)
u,\mu\right\rangle_{\Gamma} + \left\langle \bopV\lambda,\mu\right\rangle_{\Gamma}\\
\left\langle (\tfrac12\bopI + \bopKadj)
\lambda,v \right\rangle_{\Gamma} + \left\langle \bopW u,v\right\rangle_{\Gamma}.
\end{multline}

We may then rewrite <ref> and <ref> as

\begin{equation}\label{eq:realtion}
\form{C}[(u,\lambda),(v,\mu)]=\left\langle u,\mu\right\rangle_\Gamma+\left\langle\lambda,v\right\rangle_\Gamma.
\end{equation}

We will also frequently use the multitrace form, defined by

\begin{equation}\label{eq:mult_trace}
\form{A}[(u,\lambda),(v,\mu)]:=
u,\mu\right\rangle_{\Gamma} + \left\langle \bopV\lambda,\mu\right\rangle_{\Gamma}
\left\langle \bopKadj
\lambda,v \right\rangle_{\Gamma} + \left\langle \bopW u,v\right\rangle_{\Gamma}.
\end{equation}

Using this, we may rewrite <ref> as

\begin{equation}\label{eq:skewsym_relation}
\form{A}[(u,\lambda),(v,\mu)]=\tfrac12\left\langle u,\mu\right\rangle_\Gamma+\tfrac12\left\langle\lambda,v\right\rangle_\Gamma.
\end{equation}

To quantify the two traces we introduce the product space

\[
%\productspace{V} := H^{1/2}(\Gamma) \times (H^{-1/2}(\Gamma)\cap \Ltwo(\Gamma\N\cup\Gamma\R)).
\productspace{V} := \begin{cases}
H^{1/2}(\Gamma) \times H^{-1/2}(\Gamma)&\text{if }\Gamma\N\cup\Gamma\R=\emptyset,\\
H^{1/2}(\Gamma) \times \Ltwo(\Gamma)&\text{otherwise.}
\end{cases}
\]

The additional regularity on the flux variable is required later when imposing Neumann and Robin conditions.
We also introduce the associated norm

\[
\Vnorm{(v,\mu)}:= \Hnorm{1/2}{v}+\Hnorm{-1/2}{\mu}.
\]

Using the results in <ref>,
we obtain the continuity and coercivity of $\form{A}$.

There exists $C>0$ such that

\begin{align*}
\left|\form{A}[(w,\eta),(v,\mu)]\right| &\leqslant C \Vnorm{(w,\eta)}\Vnorm{(v,\mu)}&&\forall(w,\eta),(v,\mu)\in\productspace{V}.
\end{align*}

Use the stability results from <ref>.

There exists $\alpha>0$ such that

\begin{align*}
\alpha\left(
\Hstarnorm{1/2}{v}^2+\Hnorm{-1/2}{\mu}^2
\right)
\form{A}[(v,\mu),(v,\mu)]
\end{align*}

Use the coercivity of $\bopV$ and $\bopW$ from <ref> and let $\alpha=\min(\alpha_\bopW,\alpha_\bopV)$.

§ WEAK IMPOSITION OF BOUNDARY CONDITIONS

In this section, we will derive boundary integral formulations of the
problem <ref>, that we will then use for our boundary element formulations.
We assume that the boundary condition may be written as

\begin{equation}
R_\Gamma(u,\lambda) = 0.
\end{equation}

The idea that we will exploit in the following is simply to add a suitable weighted weak form of
this constraint to the Calderón form
<ref>. Formally, this leads to an expression of the

\begin{equation}\label{eq:cald_bc}
\form{C}[(u,\lambda),(v,\mu)]=\left\langle u,\mu
\right\rangle_\Gamma+\left\langle \lambda, v\right\rangle_{\Gamma} +
\left\langle R_\Gamma(u,\lambda), \betaparam_1 v + \betaparam_2 \mu \right\rangle_{\Gamma},
\end{equation}

or equivalently

\begin{equation}\label{eq:multi_bc}
\form{A}[(u,\lambda),(v,\mu)]=\tfrac12\left\langle u,\mu
\right\rangle_\Gamma+\tfrac12\left\langle \lambda, v\right\rangle_{\Gamma} +
\left\langle R_\Gamma(u,\lambda), \betaparam_1 v + \betaparam_2 \mu \right\rangle_{\Gamma},
\end{equation}

where $\betaparam_1$ and $\betaparam_2$ are problem dependent scaling operators that will be chosen as a
function of the physical parameters in order to obtain robustness of
the method.

§.§ Dirichlet boundary condition

In this section, we assume that $\Gamma\D\equiv\Gamma$ and consider the resulting Dirichlet problem.
We choose $\betaparam_1 = \betaparam\D^{1/2}$, $\betaparam_2 = \betaparam\D^{-1/2}$, where $\betaparam\D$ will
be identified with a mesh-dependent penalty parameter, and

\begin{equation}\label{eq:Dir_R}
R_{\Gamma\D}(u,\lambda) := \betaparam\D^{1/2} (g\D-u)
\end{equation}

where $g\D \in H^{1/2}(\Gamma)$ is the Dirichlet data.

Inserting this into <ref>, we obtain the formulation

\begin{multline}\label{eq:multi_D}
\form{A}[(u,\lambda),(v,\mu)] -\tfrac12 \left\langle  \lambda, v
\right\rangle_{\Gamma\D} + \tfrac12
\left\langle u, \mu \right\rangle_{\Gamma\D} +\left\langle \betaparam\D  u, v
\right\rangle_{\Gamma\D} = \left\langle g\D, \betaparam\D v + \mu
\right\rangle_{\Gamma\D}.
\end{multline}

One can compare the method with the classical (non-symmetric) Nitsche's method
by formally identifying $\lambda$ with $\partial_\B\nu u$ and $\mu$ with
$\partial_\B\nu v$ (up to the multiplicative factor $\tfrac12$).

For a more compact notation, we introduce the boundary operator
associated with the non-homogeneous Dirichlet condition

\begin{equation}\label{eq:operator_Dirichlet}
\form{B}\D[(u,\lambda),(v,\mu)]:=-\tfrac12\left\langle  \lambda, v \right\rangle_{\Gamma\D} +
\tfrac12\left\langle u, \mu \right\rangle_{\Gamma\D} +\left\langle \betaparam\D u, v
\right\rangle_{\Gamma\D},
\end{equation}

and the operator associated with the right hand side

\begin{equation}\label{eq:LDir}
\form{L}\D(v,\mu) := \left\langle  g\D,\betaparam\D v + \mu\right\rangle_{\Gamma\D}.
\end{equation}

Using these and <ref>, we arrive at the following problem: Find
$(u,\lambda) \in \productspace{V}$ such that

\begin{align}\label{eq:abstract_form_Dir}
\form{A}[(u,\lambda),(v,\mu)]+\form{B}\D[(u,\lambda),(v,\mu)] &=
\form{L}\D(v,\mu)&&\forall(v,\mu) \in \productspace{V}.
\end{align}

If we set $\betaparam\D=0$ in <ref>, we obtain a penalty-free formulation for the Dirichlet problem.

§.§ Neumann boundary condition

In this section, we assume that $\Gamma\N\equiv\Gamma$ and consider the resulting Neumann problem.
We choose $\betaparam_1 = \betaparam\N^{-1/2}$, $\betaparam_2 = \betaparam\N^{1/2}$, and define

\begin{equation}\label{eq:Neu_R}
R_{\Gamma\N}(u,\lambda) := \betaparam\N^{1/2} (g\N - \lambda),
\end{equation}

where $g\N \in \Ltwo(\Gamma\N)$, with $\int_\Gamma g\N=0$, is the Neumann data.

Proceeding as in the Dirichlet case, we obtain the formulation

\begin{multline}\label{eq:multi_N}
\form{A}[(u,\lambda),(v,\mu)] -\tfrac12 \left\langle  u, \mu
\right\rangle_{\Gamma\N} + \tfrac12
\left\langle \lambda, v \right\rangle_{\Gamma\N} +\left\langle \betaparam\N  \lambda, \mu
\right\rangle_{\Gamma\N} = \left\langle g\N, \betaparam\N \mu + v
\right\rangle_{\Gamma\N}.
\end{multline}

\begin{align}\label{eq:operator_Neumann}
\form{B}\N[(u,\lambda),(v,\mu)]&:=
-\tfrac12\left\langle u, \mu \right\rangle_{\Gamma\N}
+\tfrac12\left\langle  \lambda, v \right\rangle_{\Gamma\N}
+\left\langle \betaparam\N\lambda,\mu\right\rangle_{\Gamma\N},\\
\label{eq:LNeu}
\form{L}\N(v,\mu) &:= \left\langle  g\N,\betaparam\N \mu + v\right\rangle_{\Gamma\N},
\end{align}

we may write this as the variational problem:
Find $(u,\lambda) \in \zerostar{\productspace{V}}$ such that

\begin{align}\label{eq:abstract_form_Neu}
\form{A}[(u,\lambda),(v,\mu)]+\form{B}\N[(u,\lambda),(v,\mu)] &=
\form{L}\N(v,\mu)&&\forall(v,\mu) \in \zerostar{\productspace{V}}.
\end{align}

Here, we use the space $\zerostar{\productspace{V}}:=H_*^{1/2}(\Gamma\N)\times\Ltwo(\Gamma\N)$, as the solution to the Neumann problem can only be
determined up to a constant, so we include the extra condition that $\overline{u}=0$.

If we set $\betaparam\N=0$ in <ref>, we obtain a penalty-free formulation for the Neumann problem.
In this case, we may take
$\zerostar{\productspace{V}}=H_*^{1/2}(\Gamma\N)\times H^{-1/2}(\Gamma\N)$ and
$g\N\in H^{-1/2}(\Gamma\N)$.

When $\betaparam\N>0$, observe that for the terms imposing the Neumann condition to be well
defined, we need $\lambda \in \Ltwo(\Gamma\N)$. This can be avoided by
replacing $\betaparam\N$
with a regularising operator $\bop{R}:H^{-1/2}(\Gamma\N) \to H^{1/2}(\Gamma\N)$. For example,
we could take $\bop{R}=\betaparam_\bopV \bopV$, where
$\betaparam_\bopV \in \RR$ and $\bopV$ is the single layer boundary operator on $\Gamma\N$.
This formulation with the operator $\bop{R}$ is given in <cit.>, where it was derived
using a domain decomposition approach where a Robin condition was used to weakly impose a Neumann condition.

The resulting formulations using $\betaparam\N$ are in general easier to analyse, since
they give control of $\lambda$ on the Neumann boundary in the natural
norm $\Hnorm[\Gamma\N]{-1/2}{\lambda}$.

§.§ Mixed Dirichlet–Neumann boundary condition

We now consider the case of mixed Dirichlet–Neumann boundary conditions, when $\Gamma = \Gamma\D \cup \Gamma\N$. We note that in this
case, and in the Robin case, we take $\productspace{V} = H^{1/2}(\Gamma) \times \Ltwo(\Gamma)$.

Let $R_{\Gamma\D}$ and $R_{\Gamma\N}$ be defined by
<ref> and <ref>.
Using the abstract form <ref>, we obtain

\begin{multline}\label{eq:mixed_bc}
\form{A}[(u,\lambda),(v,\mu)]=\tfrac12\left\langle u,\mu
\right\rangle_\Gamma+\tfrac12\left\langle \lambda, v\right\rangle_{\Gamma} \\
\left\langle R_{\Gamma\D}(u,\lambda), \betaparam\D^{1/2} v + \betaparam\D^{-1/2} \mu
\right\rangle_{\Gamma\D} + \left\langle R_{\Gamma\N}(u,\lambda), \betaparam\N^{-1/2} v + \betaparam\N^{1/2} \mu
\right\rangle_{\Gamma\N}.
\end{multline}

Developing <ref>, and defining

\begin{multline}\label{eq:operator_Mixed}
\form{B}\ND[(u,\lambda),(v,\mu)] :=
\tfrac12\left\langle u, \mu \right\rangle_{\Gamma\D}
- \tfrac12\left\langle\lambda, v \right\rangle_{\Gamma\D}
+ \left\langle\betaparam\D u, v \right\rangle_{\Gamma\D} \\
\tfrac12\left\langle\lambda, v \right\rangle_{\Gamma\N}
- \tfrac12\left\langle u, \mu \right\rangle_{\Gamma\N}
+ \left\langle\betaparam\N \lambda, \mu \right\rangle_{\Gamma\N},
\end{multline}

\begin{equation}\label{eq:LMix}
\form{L}\ND(v,\mu) := \left\langle g\D, \betaparam\D v + \mu
\right\rangle_{\Gamma\D}+ \left\langle  g\N,\betaparam\N\mu + v\right\rangle_{\Gamma\N},
\end{equation}

we arrive the variational formulation: Find $(u,\lambda) \in
\productspace{V}$
such that

\begin{align}\label{eq:MND_compact}
\form{A}[(u,\lambda),(v,\mu)] +
\form{B}\ND[(u,\lambda),(v,\mu)] &= \form{L}\ND(v,\mu)
&&\forall(v,\mu) \in
\productspace{V}.
\end{align}

If we set $\betaparam\D=0$ and $\betaparam\N=0$ in <ref>, we obtain a penalty-free formulation for the mixed Dirichlet–Neumann problem.
By taking $\Gamma\N=\emptyset$ or $\Gamma\D=\emptyset$, formulations for both Dirichlet and Neumann problems can be obtained from <ref>.

§.§ Robin conditions

For simplicity, we consider the case where $\Gamma=\Gamma\R$.
Considering the Robin condition <ref>, we may write, for some $\param>0$,

\begin{equation}
R_{\Gamma\R}(u,\lambda) := \betaparam\R^{1/2}\left(\param^{1/2}(g\N-\lambda) + \param^{-1/2}(g\D-u)\right).
\end{equation}

This function is a linear combination of the Dirichlet and the Neumann conditions.

\begin{equation}\label{eq:Robin_bc}
R_{\Gamma\R}(u,\lambda) = \alpha\D R_{\Gamma\D}(u,\lambda) + \alpha\N R_{\Gamma\N}(u,\lambda),
\end{equation}

where $\alpha\N=\betaparam\R^{1/2}\betaparam\N^{-1/2}\param^{1/2}$ and $\alpha\D=\betaparam\R^{1/2}\betaparam\D^{-1/2}\param^{-1/2}$.

We take $\betaparam_1=\betaparam\R^{1/2}$ and $\betaparam_2=\betaparam\R^{-1/2}$, and look for a term of the form

\begin{equation}
\left\langle \phi R_{\Gamma\R}(u,\lambda) ,
\betaparam\R^{1/2} v + \betaparam\R^{-1/2} \mu \right\rangle_{\Gamma\R},
\end{equation}

where the $\phi$ and $\betaparam\R$ must have the following properties to ensure that
the formulation degenerates into the formulation for the Dirichlet and Neumann problems as $\param\to0$ and $\param\to\infty$.

\begin{align*}
&&\alpha\D \phi \to 1,
&&\alpha\N \phi \to 0
&&\text{as }\param\to 0,\\
&&\alpha\N \phi \to 1,
&&\alpha\D \phi \to 0
&&\text{as }\param\to \infty.
\end{align*}

It is straightforward to verify that these conditions are satisfied
for the choices

\begin{align}
\phi&:= \frac{\param^{1/2}}{\param\betaparam\R + 1},\\
\betaparam\R&:=\frac{\param\betaparam\N^{-1}+\betaparam\D}{\param+1}.
\end{align}

Later, we will use
$\betaparam\D=\betaparam h^{-1}$
$\betaparam\N=\betaparam h$,
where $\betaparam$ is a constant, as in the mixed Dirichlet–Neumann case.

Collecting the above considerations, we arrive at the formulation

\begin{multline}\label{eq:cald_Rob}
\form{A}[(u,\lambda),(v,\mu)] =
\tfrac12\left\langle u,\mu \right\rangle_{\Gamma}+\tfrac12\left\langle \lambda, v\right\rangle_{\Gamma} \\
+ \left\langle
\param(g\N-\lambda) + (g\D-u),
\frac{\betaparam\R}{\param \betaparam\R+1} v +  \frac{1}{\param \betaparam\R+1} \mu
\right\rangle_{\Gamma\R}.
\end{multline}

Taking $\param\to0$, we recover the Dirichlet formulation <ref>;
and taking $\param \to \infty$ results in the Neumann formulation  <ref>.

By introducing

\begin{multline*}
\form{B}\R[(u,\lambda),(v,\mu)]:=
\frac12\left\langle \frac{\param
\betaparam\R -1}{\param
\betaparam\R +1} \lambda, v\right\rangle_{\Gamma\R}
- \frac12\left\langle \frac{\param
\betaparam\R -1}{\param
\betaparam\R +1} u, \mu\right\rangle_{\Gamma\R} \\
+\left\langle \frac{\param}{\param
\betaparam\R +1} \lambda, \mu \right\rangle_{\Gamma\R}
+\left\langle \frac{\betaparam\R}{\param
\betaparam\R +1} u,v\right\rangle_{\Gamma\R}
\end{multline*}

\[
\form{L}\R(v,\mu)
\left\langle g\D+ \param g\N, \frac{\betaparam\R}{\param \betaparam\R +1} v  +
\frac{1}{\param \betaparam\R +1} \mu\right\rangle_{\Gamma\R},
\]

we may write this as the variational problem: Find $(u,\lambda) \in\productspace{V}$ such that

\begin{align}\label{eq:cald_Rob_compact}
\form{A}[(u,\lambda),(v,\mu)] +
\form{B}\R[(u,\lambda),(v,\mu)] &= \form{L}\R(v,\mu) &&\forall
(v,\mu) \in\productspace V.
\end{align}

§ BOUNDARY ELEMENT METHOD FOR THE SINGLE DOMAIN PROBLEM

All the methods introduced above are written as the sum of the multitrace operator $\form{A}$ and a
boundary condition operator $\form{B}$. We write this generally as: Find
$(u,\lambda) \in \productspace{V}$ such that

\begin{align}\label{eq:abstract_form}
\form{A}[(u,\lambda),(v,\mu)]+\form{B}[(u,\lambda),(v,\mu)] &=
\form{L}(v,\mu) &&\forall (v,\mu) \in \productspace{V}.
\end{align}

In this section, we analyse this general problem, then show that the analysis is applicable to the
boundary conditions discussed in <ref>.

For the sake of example and to fix the ideas, we introduce a family of conforming, shape regular triangulations of
$\Gamma$, $\{\Th\}_{h>0}$, indexed by the largest element
diameter of the mesh, $h$.
We assume that the triangulations are fitted to the different
boundary sets $\Gamma\D$, $\Gamma\R$ and $\Gamma\N$.
We then consider the following finite element spaces.

\begin{align*}
\Vee^k_h &:= \{v_h \in C^0(\Gamma): v_h\vert_{T} \in \mathbb{P}_k(T) \text{, for every }T\in\Th\},\\
\Lambda^l_h &:= \{v_h \in \Ltwo(\Gamma): v_h\vert_{T} \in \mathbb{P}_l(T) \text{, for every }T\in\Th\},\\
\tilde \Lambda^l_h &:= \{v_h \in \Lambda^l_h: v_h\vert_{\Gamma_i} \in
C^0(\Gamma_i)\text{, for }i=1,\hdots,M\},
\end{align*}

where $\mathbb{P}_k(T)$ denotes the space of polynomials of order less
than or equal to $k$,
and $\{\Gamma_i\}_{i=1}^M$ are the polygonal faces of $\Gamma$.

We observe that $\Vee^k_h \subset H^{1/2}(\Gamma)$, $\Lambda^l_h
\subset \Ltwo(\Gamma)$ and $\tilde \Lambda^l_h
\subset \Ltwo(\Gamma)$. We now introduce the discrete product
space $\productspace{V}_h :=  \Vee^k_h \times \Lambda^l_h$. The space $\tilde
\Lambda^l_h$ may be used in the place of $\Lambda^l_h$ without any modifications of the arguments below.

The boundary element formulation of the generic problem <ref> then takes the form: Find
$(u_h,\lambda_h) \in  \productspace{V}_h$ such that

\begin{align}\label{BEMform}
\form{A}[(u_h,\lambda_h),(v_h,\mu_h)] +
\form{B}[(u_h,\lambda_h),(v_h,\mu_h)] &= \form{L}(v_h,\mu_h)
&&\forall (v_h, \mu_h) \in \productspace{V}_h.
\end{align}

If we assume that
$(u,\lambda) \in \productspace{V}$ and $(u_h,\lambda_h) \in \productspace{V}_h$
<ref> and <ref>,
it immediately follows that the following Galerkin orthogonality relation holds.

\begin{multline}\label{galortho}
\form{A}[(u-u_h,\lambda-\lambda_h),(v_h,\mu_h)] +
\form{B}[(u-u_h,\lambda-\lambda_h),(v_h,\mu_h)] = 0
\\
\forall (v_h, \mu_h) \in \productspace{V}_h.
\end{multline}

We also get the following representation formula for the approximation in the
bulk using <ref>.

\begin{equation}\label{eq:rec_disc}
\tilde u_h = -\popK u_h + \popV\lambda_h.
\end{equation}

We will now proceed to derive some estimates for the solution of
<ref> and the reconstruction <ref>.

Let $\productspace{W}$ be a product Hilbert space for the primal and flux
variables, such that $\productspace{W} \subset \productspace{V}$.
Let $\Bnorm{\cdot}$ be a norm defined on $\productspace{W}$, such that for all $(v,\mu)\in\productspace{W}$,
$\Bnorm{(v,\mu)} \geqslant \Vnorm{(v,\mu)}$.

To reduce the number of constants that appear, especially when proving that <ref> holds, we introduce the following notation.

* If $\exists C>0$, independent of $h$, such that $a\leqslant Cb$, then we write $a\lesssim b$.

* If $a\lesssim b$ and $b\lesssim a$, then we write $a\eqsim b$.

For the abstract analysis, we will make use of the following standard assumptions.

[Weak coercivity]
There exists $\alpha>0$ such that $\forall (v,\mu) \in \productspace{W}$

\[\alpha \Bnorm{(v,\mu)} \leqslant \sup_{(w,\eta) \in \productspace{W}\setminus\{0\}}\frac{\form{A}[(v,\mu),(w,\eta)] +
\form{B}[(v,\mu),(w,\eta)]}{\Bnorm{(w,\eta)}},\]

and $\forall (w,\eta) \in \productspace{W}\setminus\{0\}$

\[\sup_{(v,\mu) \in \productspace{W}}\left|\form{A}[(v,\mu),(w,\eta)] +
\form{B}[(v,\mu),(w,\eta)]\right|>0.\]

[Discrete coercivity]
There exists $\alpha>0$ such that $\forall (v_h,\mu_h) \in \productspace{V}_h$

\[\alpha \Bnorm{(v_h,\mu_h)} \leqslant \sup_{(w_h,\eta_h) \in \productspace{V}_h\setminus\{0\}}\frac{\form{A}[(v_h,\mu_h),(w_h,\eta_h)] +
\form{B}[(v_h,\mu_h),(w_h,\eta_h)]}{\Bnorm{(w_h,\eta_h)}},\]

and $\forall (w_h,\eta_h) \in \productspace{V}_h\setminus\{0\}$

\[\sup_{(v_h,\mu_h) \in \productspace{V}_h}\left|\form{A}[(v_h,\mu_h),(w_h,\eta_h)] +
\form{B}[(v_h,\mu_h),(w_h,\eta_h)]\right|>0.\]

There exists an auxiliary norm $\starnorm{(v,\mu)}$ defined on $\productspace{W}$, and
there exists
$M>0$ such that $\forall(w,\eta),(v,\mu) \in \productspace{W}$

\[\left|\form{A}[(w,\eta),(v,\mu)] +
\form{B}[(w,\eta),(v,\mu)]\right| \leqslant M\starnorm{(w,\eta)}\Bnorm{(v,\mu)}
\]

$\forall (v,\mu) \in H^s(\Gamma)\times H^r(\Gamma)$,

\begin{equation*}
\inf_{(w_h,\eta_h) \in \productspace{V}_h}
\starnorm{(v-w_h,\mu - \eta_h)}
\lesssim h^{\zeta-1/2} \Hseminorm{\zeta}{v} +
h^{\xi+1/2} \Hseminorm{\xi}{\mu},
\end{equation*}

where $\zeta = \min(k+1,s)$, $\xi = \min(l+1,r)$,
$s\geqslant\frac12$ and $r\geqslant-\frac12$.

Typically, we use approximation spaces with $k=l+1$, where the polynomial spaces used for $\lambda$ are one order lower than those for $u$,
or spaces with $k=l$, where equal order spaces are used for both variables.

We note that if the form $\form{A}+\form{B}$ is coercive,
that is there exists $\alpha>0$ such that $\forall (v,\mu) \in \productspace{W}$

\[\alpha \Bnorm{(v,\mu)}^2 \leqslant \form{A}[(v,\mu),(v,\mu)] + \form{B}[(v,\mu),(v,\mu)],\]

then <ref> hold.

We now proceed to prove some results about the abstract problem.

Assume that <ref> holds, then the linear system defined by <ref> is invertible.
If, in addition, we assume that

* <ref> holds,

* there exists $L>0$ such that $\form{L}(w,\eta) \leqslant L \Bnorm{(w,\eta)}\quad\forall(w,\eta) \in \productspace{W}$,

* and $\starnorm{\cdot}$ is equivalent to $\Bnorm{\cdot}$,

then the formulation <ref> admits a unique solution
in $\productspace{W}$.

Note that <ref> implies the inf-sup condition,

\begin{align}
\inf_{(v,\mu)\in\productspace{W}\setminus\{0\}}\sup_{(w,\eta)\in\productspace{W}\setminus\{0\}}
\frac{\form{A}[(v,\mu),(w,\eta)] + \form{B}[(v,\mu),(w,\eta)]}{
\Bnorm{(v,\mu)}\Bnorm{(w,\eta)}
\end{align}

Therefore we may apply the Babuška–Lax–Milgram theorem <cit.>.

Assume that $(u,\lambda) \in \productspace{V}$ is the solution to a boundary
value problem of the form <ref> satisfying the
abstract form <ref>. Let $(u_h,\lambda_h) \in
\productspace{V}_h$ be the solution of <ref>.
If <ref> are

\begin{equation}\label{eq:best}
\Bnorm{(u-u_h,\lambda-\lambda_h)} \leqslant\frac{M}\alpha
\inf_{(v_h,\mu_h)\in \productspace{V}_h}
\starnorm{(u-v_h,\lambda-\mu_h)}.
\end{equation}

See <cit.>.

Let $(u,\lambda) \in H^s(\Gamma) \times H^r(\Gamma)$, for some $s\geqslant\frac12$ and $r\geqslant-\frac12$, satisfy the
abstract form <ref>.
Under the assumptions of <ref> and <ref>,

\[
\Bnorm{(u-u_h,\lambda-\lambda_h)} \lesssim
h^{\zeta-1/2} \Hseminorm{\zeta}{u}
+ h^{\xi+1/2}\Hseminorm{\xi}{\lambda},
\]

where $\zeta = \min(k+1,s)$ and $\xi = \min(l+1,r)$.

Apply <ref> to the right hand side of <ref>.

Assume that $(u,\lambda) \in \productspace{V}$ is the solution to a boundary
value problem of the form <ref> satisfying the
abstract form <ref> and that the assumptions of
<ref> are satisfied. Let $(u_h,\lambda_h) \in
\productspace{V}_h$. Let $\tilde u:\Omega \to \RR$ be the
reconstruction obtained using <ref>, with $\trace\N
u=\lambda$ and $\trace\D u=u$; and $\tilde u_h:\Omega \to \RR$
be the reconstruction obtained using <ref>. Then there

\[
\Hnorm[\Omega]{1}{\tilde u-\tilde u_h} \lesssim \frac{M}{\alpha} \inf_{v_h,\mu_h
\in \productspace{V}_h} \starnorm{(u-v_h,\lambda-\mu_h)}.
\]

Using <ref> and <ref>, we may write

\[
\tilde u - \tilde u_h = (u^\popV_\lambda - u^\popV_{\lambda_h})+(u^\popK_u - u^\popK_{u_h}).
\]

Using the triangle inequality, we have

\begin{equation}\label{eq:triangle}
\Hnorm[\Omega]{1}{\tilde u-\tilde u_h} \leqslant \Hnorm[\Omega]{1}{u^\popV_\lambda - u^\popV_{\lambda_h}}+\Hnorm[\Omega]{1}{u^\popK_u - u^\popK_{u_h}}.
\end{equation}

By <ref> and <ref>, there exist $c_1,c_2>0$ such that

\begin{align}
\Hnorm[\Omega]{1}{u^\popV_\lambda - u^\popV_{\lambda_h}} &\leqslant c_1\Hnorm{-1/2}{\lambda - \lambda_h},\label{eq:staberr1}\\
\Hnorm[\Omega]{1}{u^\popK_u- u^\popK_{u_h}} &\leqslant c_2\Hnorm{1/2}{u - u_h}.\label{eq:staberr2}
\end{align}

Collecting <ref>, we see that there exists $C>0$ such that

\begin{equation}\label{eq:final_line_of_proof}
\Hnorm[\Omega]{1}{\tilde u-\tilde u_h} \leqslant C\Vnorm{\lambda - \lambda_h,u-u_h}\leqslant C\Bnorm{\lambda - \lambda_h,u-u_h}.
\end{equation}

The statement now follows from <ref>.

Under the same assumptions of <ref> and <ref>,

\[
\Hnorm[\Omega]{1}{\tilde u-\tilde u_h} \lesssim
h^{\zeta-1/2} \Hseminorm{\zeta}{u}
+ h^{\xi+1/2}\Hseminorm{\xi}{\lambda},
\]

where $\zeta = \min(k+1,s)$ and $\xi = \min(l+1,r)$.

Apply <ref> to <ref> in the proof of <ref>.

§.§ Application of the theory to the Dirichlet problem

For the finite element spaces defined above, the Dirichlet problem
takes the form: Find
$(u_h,\lambda_h) \in\productspace{V}_h$ such that

\begin{align}\label{BEMform_Dir}
\form{A}[(u_h,\lambda_h),(v_h,\mu_h)] +
\form{B}\D[(u_h,\lambda_h),(v_h,\mu_h)] &= \form{L}\D(v_h,\mu_h)
&&\forall(v_h,\mu_h) \in \productspace{V}_h.
\end{align}

We introduce the following $\form{B}\D$-norm.

\[
\Bnorm[\D]{(v,\mu)} := \Vnorm{(v,\mu)} + \betaparam\D^{1/2} \Ltwonorm[\Gamma\D]{v},
\]

and we let $\starnorm{\cdot}=\Bnorm[\D]{\cdot}$.
We now proceed to verify that <ref> hold.

<Ref> are satisfied for the Dirichlet problem if $\exists\betamin>0$, independent of $h$,
such that $\betaparam\D>\betamin$.

Using the fact that $\Hstarnorm[\Gamma\D]{1/2}{v}^2+\Ltwonorm[\Gamma\D]{\overline{v}}^2\gtrsim\Hnorm[\Gamma\D]{1/2}{v}^2$, we deduce from <ref>
that for every positive $\alpha'\leqslant\alpha$,

\begin{align*}
\alpha'
\Vnorm{(v,\mu)}^2
\form{A}[(v,\mu),(v,\mu)]
\end{align*}

Using the definition of $\form{B}\D$, we see that

\begin{align*}
\form{B}\D[(v,\mu),(v,\mu)]&=\betaparam\D\langle v,v\rangle_{\Gamma\D}=\betaparam\D\Ltwonorm[\Gamma\D]{v}^2
\end{align*}

Taking $\alpha'=\min(\alpha,\betamin/2)$, we see that

\begin{align*}
\form{A}[(v,\mu),(v,\mu)]+\form{B}\D[(v,\mu),(v,\mu)]
\end{align*}

for some $\alpha''>0$.
Therefore, in this case the form $\form{A}+\form{B}\D$ is coercive, and so <ref> hold.

<Ref> are satisfied for the Dirichlet problem with $\betaparam\D=0$.

Taking $w=v$ and $\eta=\mu+c\overline{v}$, for some $c\in\RR$ to be fixed, we obtain

\begin{align}
\form{B}\D[(v,\mu),(w,\eta)]\nonumber\\
\left\langle \bopV \mu,\mu\right\rangle_{\Gamma}
+ c\left\langle \bopV \mu,\overline{v}\right\rangle_{\Gamma}
- c\left\langle \bopK v,\overline{v} \right\rangle_{\Gamma}
+ \left\langle \bopW v,v\right\rangle_{\Gamma}
+ \frac{c}2\left\langle v,\overline{v}\right\rangle_\Gamma\label{eq:dirichlet_zero_first_bit}
\end{align}

By <ref>, we know that

\begin{align}\label{diricostep1}
\left\langle \bopV \mu,\mu\right\rangle_{\Gamma}
+ \left\langle \bopW v,v\right\rangle_{\Gamma}
\geqslant
\alpha_\bopV\Hnorm{-1/2}{\mu}^2
\end{align}

By <ref>, we see that

\begin{align*}
c\left|\left\langle \bopV\mu, \overline{v}\right\rangle_{\Gamma}\right|
c\Hnorm{1/2}{\bopV \mu}\Hnorm{-1/2}{\overline{v}}\\
\end{align*}

Using the fact that for $a,b\geqslant0$, $ab\leqslant(a^2+b^2)/2$, we obtain

\begin{align}
c\left|\left\langle \bopV\mu, \overline{v}\right\rangle_{\Gamma}\right|
\frac{c^2C_\bopV^2}{2\alpha_\bopV}\Ltwonorm{\overline{v}}^2+\frac{\alpha_\bopV}{2}\Hnorm{-1/2}{\mu}^2.
\label{diricostep2}
\end{align}

We note that $u=\overline{v}$ is a solution to <ref>, $\trace\D\overline{v}=\overline{v}$ and $\trace\N \overline{v}=0$. Using this and applying <ref>, we see that
$\forall\mu\in H^{-1/2}(\Gamma)$, $\left\langle\bopK\overline{v},\mu\right\rangle_{\Gamma}=-\tfrac12\left\langle\overline{v},\mu\right\rangle_{\Gamma}$.
Therefore, using $\mu=\overline{v}$,

\begin{align*}
c\left\langle \bopK v, \overline{v} \right\rangle_{\Gamma}
&=c\left\langle \bopK (v-\overline{v}), \overline{v} \right\rangle_{\Gamma}
+c\left\langle \bopK \overline{v}, \overline{v} \right\rangle_{\Gamma}
\\
&=c\left\langle \bopK (v-\overline{v}), \overline{v} \right\rangle_{\Gamma}
-\frac{c}2\left\langle \overline{v}, \overline{v} \right\rangle_{\Gamma}.
\end{align*}

Using the fact that $\Hnorm{1/2}{v-\overline{v}}=\Hstarnorm{1/2}{v}$, and
proceeding in the same way as we did for the single layer term above, we obtain

\begin{align}
c\left\langle \bopK v,\overline{v} \right\rangle_{\Gamma}
%&\leqslant c\Hnorm{1/2}{\bopK(v-\overline{v})}\Hnorm{-1/2}{\overline{v}}-\frac{c}2\Ltwonorm{\overline{v}}^2
%&\leqslant cC_\bopK\Hnorm{1/2}{v-\overline{v}}\Ltwonorm{\overline{v}}-\frac{c}2\Ltwonorm{\overline{v}}^2
%&= cC_\bopK\Hstarnorm{1/2}{v}\Ltwonorm{\overline{v}}-\frac{c}2\Ltwonorm{\overline{v}}^2
\frac{\alpha_\bopW}{2}\Hstarnorm{1/2}{v}^2
\label{diricostep3}
\end{align}

We also have that

\begin{align}\label{diricostep4}
\frac{c}2\left\langle v,\overline{v}\right\rangle
\end{align}

Taking $\alpha=\min(\alpha_\bopV,\alpha_\bopK)$ and $C=\max(C_\bopV,C_\bopK)$, and
putting <ref> together, we obtain

\begin{align*}
\frac{\alpha}2\Hnorm{-1/2}{\mu}^2
\right)\Ltwonorm{\overline{v}}^2.
\end{align*}

Letting $\displaystyle c=\frac{\alpha}{2C^2}$ gives

\begin{align*}
\frac{\alpha}2\Hnorm{-1/2}{\mu}^2
\Hnorm{-1/2}{\mu}^2
\end{align*}

Finally, we show that

\begin{align*}
\Vnorm{(v,\mu)}
\Vnorm{(w,\eta)}
\end{align*}

\begin{align*}
\Vnorm{(v,\mu)}\Vnorm{(w,\eta)}
\Hnorm{-1/2}{\mu}^2
&\lesssim L.
\end{align*}

We obtain the first part of <ref> by dividing through by $\Vnorm{(w,\eta)}$ and taking the supremum.

To show the second part of <ref>, we let $(w,\eta)\in\productspace{W}\setminus\{0\}$ and proceed as follows.

\begin{align*}
L&:=\sup_{(v,\mu) \in \productspace{W}}\left|\form{A}[(v,\mu),(w,\eta)] + \form{B}\D[(v,\mu),(w,\eta)]\right|\\
\form{A}[(w,\eta-\overline{w}),(w,\eta)] + \form{B}\D[(w,\eta-\overline{w}),(w,\eta)]\\
+\langle\bopW w,w\rangle_\Gamma
\end{align*}

This is of the same form as <ref>, so we proceed as above to obtain

\[
\]

This is greater than zero for all $(w,\eta)\not=0$, and so we have proven the second part of <ref>.

<Ref> can be proven in the same way as above using the discrete space $\productspace{V}_h$ in the place of $\productspace{W}$.

<Ref> is satisfied for the Dirichlet problem.

Applying <ref>, the relation

$$\left\langle \eta,v \right\rangle_\Gamma \leqslant \Hnorm{-1/2}{\eta} \Hnorm{1/2}{v},$$

and the Cauchy–Schwarz inequality,

$$\betaparam\D \left\langle  w,v \right\rangle_{\Gamma} \leqslant \betaparam\D^{1/2} \Ltwonorm{w}
\betaparam\D^{1/2} \Ltwonorm{v},$$

to the form
$\form{A}+\form{B}\D$ yields the desired continuity result.

<Ref> is satisfied for the Dirichlet problem if $0\leqslant\betaparam\D\lesssim h^{-1}$.

Using standard approximation results
(see eg <cit.>), we see that

\begin{align*}
\inf_{(w_h,\eta_h)\in\productspace{V}_h}\Vnorm{(v-w_h,\mu-\eta_h)}
&=\inf_{w_h\in \Vee^k_h}\Hnorm{1/2}{v-w_h}+\inf_{\eta_h\in\Lambda^l_h}\Hnorm{-1/2}{\mu-\eta_h}\\
&\lesssim h^{\zeta-1/2}\Hseminorm{\zeta}{v}
+ h^{\xi+1/2}\Hseminorm{\xi}{\mu},\\
\inf_{w_h\in\Vee^k_h}\Ltwonorm[\Gamma\D]{v-w_h}&\lesssim h^{\zeta}\Hseminorm{\zeta}{v}.
\end{align*}

Applying these to the definition of $\starnorm{\cdot}$ gives

\begin{equation*}
\inf_{(w_h,\eta_h)\in\productspace{V}_h}\starnorm{(v-w_h,\mu-\eta_h)}
\lesssim h^{\zeta-1/2}\Hseminorm{\zeta}{v}
+ h^{\xi+1/2}\Hseminorm{\xi}{\mu}
+ \betaparam\D^{1/2}h^{\zeta}\Hseminorm{\zeta}{v}.
\end{equation*}

If $\betaparam\D=0$, <ref> holds.
If $0<\betaparam\D\lesssim h^{-1}$, then $\betaparam\D^{1/2}h^{\zeta}\lesssim h^{\zeta-1/2}$, and so <ref> holds.

We have shown that <ref> are satisfied. Additionally the extra assumptions in <ref>
are satisfied, so we conclude that the results of
<ref> apply to the Dirichlet problem.
This is summarised in the following result.

The Dirichlet problem <ref> has a unique solution $(u,\lambda)\in H^{s}(\Gamma)\times H^{r}(\Gamma)$, for some $s\geqslant\tfrac12$ and $r\geqslant-\tfrac12$.
The discrete Dirichlet problem <ref> is invertible. If
$\exists\betamin>0$ such that $\betamin<\betaparam\D\lesssim h^{-1}$
or $\betaparam\D=0$, its solution $(u_h,\lambda_h)\in\Vee_h^k\times\Lambda_h^l$ satisfies

\[
\Bnorm[\D]{(u-u_h,\lambda-\lambda_h)} \lesssim
h^{\zeta-1/2} \Hseminorm{\zeta}{u} + h^{\xi+1/2}
\Hseminorm{\xi}{\lambda},
\]

where $\zeta = \min(k+1,s)$ and $\xi = \min(l+1,r)$. Additionally,

\[
\Hnorm[\Omega]{1}{\tilde u-\tilde u_h} \lesssim
h^{\zeta-1/2} \Hseminorm{\zeta}{u} + h^{\xi+1/2}
\Hseminorm{\xi}{\lambda},
\]

where $\tilde u$ and $\tilde u_h$ are the solutions in $\Omega$ computed using <ref>.

§.§ Application of the theory to the Neumann problem

The Neumann problem takes the form: Find
$(u_h,\lambda_h) \in\zerostar{\productspace{V}}_h$ such that

\begin{align}\label{BEMform_Neu}
\form{A}[(u_h,\lambda_h),(v_h,\mu_h)] +
\form{B}\N[(u_h,\lambda_h),(v_h,\mu_h)] &= \form{L}\N(v_h,\mu_h)
&&\forall(v_h,\mu_h) \in \zerostar{\productspace{V}}_h.
\end{align}

Here $\zerostar{\productspace{V}}_h:=\zerostar{\Vee}^k_h(\Gamma)\times\Lambda^l_h(\Gamma)$
and $\zerostar{\Vee}^k_h(\Gamma):=\{v\in\Vee^k_h:\overline{v}=0\}$.

We introduce the following $\form{B}\N$-norm.

\[
\Bnorm[\N]{(v,\mu)} := \Vnorm{(v,\mu)} + \betaparam\N^{1/2} \Ltwonorm[\Gamma\N]{\mu},
\]

and we let $\starnorm{\cdot}=\Bnorm[\N]{\cdot}$.

We now proceed to verify that <ref> hold.

<Ref> are satisfied for the Neumann problem with $\betaparam\N\geqslant0$.

As $v\in H_*^{1/2}(\Gamma\N)$, we may immediately apply <ref> to show that the form is coercive.

<Ref> is satisfied for the Neumann problem.

The proof is the same as in the Dirichlet case.

<Ref> is satisfied for the Neumann problem if $0\leqslant\betaparam\N\lesssim h$.

The proof is the same as in the Dirichlet case.

As in the Dirichlet case,
the extra assumptions in <ref> are satisfied. We therefore conclude with the following result.

The Neumann problem <ref> has a unique solution $(u,\lambda)\in H_*^{s}(\Gamma)\times H^{r}(\Gamma)$,
for some $s\geqslant\tfrac12$ and $r\geqslant0$ if $\betaparam\N>0$.
If $\betaparam\N=0$, this holds for some $r\geqslant-\tfrac12$.
The discrete Neumann problem <ref> is invertible. If $0\leqslant\betaparam\N\lesssim h$, its solution $(u_h,\lambda_h)\in\zerostar{\Vee}_h^k\times\Lambda_h^l$ satisfies

\[
\Bnorm[\N]{(u-u_h,\lambda-\lambda_h)} \lesssim
h^{\zeta-1/2} \Hseminorm{\zeta}{u} + h^{\xi+1/2}
\Hseminorm{\xi}{\lambda},
\]

where $\zeta = \min(k+1,s)$ and $\xi = \min(l+1,r)$. Additionally,

\[
\Hnorm[\Omega]{1}{\tilde u-\tilde u_h} \lesssim
h^{\zeta-1/2} \Hseminorm{\zeta}{u} + h^{\xi+1/2}
\Hseminorm{\xi}{\lambda},
\]

where $\tilde u$ and $\tilde u_h$ are the solutions in $\Omega$ computed using <ref>.

§.§ Application of the theory to the mixed Dirichlet–Neumann problem

For the mixed problem, the boundary element method takes the form: Find $(u_h,\lambda_h) \in
\productspace{V}_h$ such that

\begin{multline}\label{eq:Nitsche_mixed_BEM}
\form{A}[(u_h,\lambda_h),(v_h,\mu_h)]
+\form{B}\ND[(u_h,\lambda_h),(v_h,\mu_h)] = \form{L}\ND(v_h,\mu_h)
\\\forall (v_h,\mu_h) \in
\productspace{V}_h.
\end{multline}

We now show that the assumptions for the abstract error estimate are
satisfied for the formulation <ref>. First, we
introduce the following norms.

\begin{align*}
\Bnorm[\ND]{(v,\mu)} &:= \Vnorm{(v,\mu)} + \betaparam\D^{1/2} \Ltwonorm[\Gamma\D]{v} + \betaparam\N^{1/2} \Ltwonorm[\Gamma\N]{\mu}\\
\starnorm{(v,\mu)}&:=\Vnorm{(v,\mu)}+ \betaparam\D^{1/2} \Ltwonorm{v}+\betaparam\N^{1/2}\Ltwonorm{\mu}.
\end{align*}

Observe that in this case the two norms are not the same, nor are they
equivalent, so the below results cannot be used to prove existence of
a unique solution to <ref>. Nevertheless, it is easy to verify
that if the exact solution to the mixed Dirichlet–Neumann problem is in $\productspace{V}$
then it satisfies <ref>.

<Ref> are satisfied for the mixed Dirichlet–Neumann problem if $\exists\betamin>0$,
independent of $h$, such that $\betaparam\D>\betamin$.

We obtain using <ref> that for $(v,\mu) \in \productspace{W}$,

\begin{align*}
\form{A}[(v,\mu),(v,\mu)] +
\form{B}\ND[(v,\mu),(v,\mu)]
\\
\alpha\Hnorm{-1/2}{\mu}^2
+ \alpha\Hstarnorm{1/2}{v}^2
+ \betaparam\D \Ltwonorm[\Gamma\D]{v}^2
+ \betaparam\N \Ltwonorm[\Gamma\N]{\mu}^2.
\end{align*}

Taking $\alpha'=\min(\alpha,\betamin/2)$, we get

\begin{multline*}
\alpha'\Hnorm{-1/2}{\mu}^2
+ \alpha'\left(\Hstarnorm{1/2}{v}^2+\Ltwonorm[\Gamma\D]{v}^2\right)
\\+ (\betaparam\D-\alpha') \Ltwonorm[\Gamma\D]{v}^2
+ \betaparam\N \Ltwonorm[\Gamma\N]{\mu}^2.
\end{multline*}

By <cit.>,
$\left(\Hstarnorm{1/2}{\cdot}^2+\Ltwonorm[\Gamma\D]{\cdot}^2\right)^{1/2}$ is an equivalent norm to
$\Hnorm{1/2}{\cdot}$. Therefore

\begin{align*}
\alpha'\Hnorm{-1/2}{\mu}^2
+ \alpha'\Hnorm{1/2}{v}^2
+ \betaparam\D\left(1-\frac{\alpha'}{\betamin}\right) \Ltwonorm[\Gamma\D]{v}^2
+ \betaparam\N \Ltwonorm[\Gamma\N]{\mu}^2
\\
\Hnorm{-1/2}{\mu}^2
+ \Hnorm{1/2}{v}^2
+ \betaparam\D \Ltwonorm[\Gamma\D]{v}^2
+ \betaparam\N \Ltwonorm[\Gamma\N]{\mu}^2
\end{align*}

Coercivity follows using the definition of $\Bnorm[\ND]{\cdot}$.

<Ref> is satisfied for the mixed Dirichlet–Neumann problem if
$\exists\betamin>0$, independent of $h$, such that $\betaparam\D^{1/2}\betaparam\N^{1/2}>\betamin$.

Using the fact that $\left\langle v,\mu\right\rangle_\Gamma=\left\langle v,\mu\right\rangle_{\Gamma\D}+\left\langle v,\mu\right\rangle_{\Gamma\N}$,
we see that

\begin{align*}
\form{B}\ND[(w,\eta),(v,\mu)]\vspace{-200mm}&=
\tfrac12\left\langle w,\mu\right\rangle_{\Gamma\D}
+\betaparam\D\left\langle w,v\right\rangle_{\Gamma\D}
\\&\qquad\qquad\qquad
-\tfrac12\left\langle w,\mu\right\rangle_{\Gamma\N}
\tfrac12\left\langle w,\mu\right\rangle_\Gamma
+\betaparam\D\left\langle w,v\right\rangle_{\Gamma\D}
\\&\qquad\qquad\qquad\quad
-\left\langle w,\mu\right\rangle_{\Gamma\N}
\tfrac12\left\langle w,\mu\right\rangle_\Gamma
+\betaparam\D\left\langle w,v\right\rangle_{\Gamma\D}
\\&\qquad\qquad\qquad\quad
-\betaparam\D^{1/2}\betaparam\N^{1/2}\left\langle w,\mu\right\rangle_{\Gamma\N}
\end{align*}

Proceeding as in <ref> leads to the desired result.

<Ref> is satisfied for the mixed Dirichlet–Neumann problem if $0<\betaparam\D\lesssim h^{-1}$ and $0<\betaparam\N\lesssim h$.

Proceeding as in the Dirichlet case, we see that

\begin{multline*}
\inf_{(w_h,\eta_h)\in\productspace{V}_h}\starnorm{(v-w_h,\mu-\eta_h)}
\lesssim h^{\zeta-1/2}\Hseminorm{\zeta}{v} + h^{\xi+1/2}\Hseminorm{\xi}{\mu}\\
+ \betaparam\D^{1/2}h^{\zeta}\Hseminorm{\zeta}{v}
+ \betaparam\N^{1/2}h^{\xi}\Hseminorm{\xi}{\mu}
\end{multline*}

If $0<\betaparam\D\lesssim h^{-1}$ and $0<\betaparam\N\lesssim h$, then

\begin{equation*}
\betaparam\D^{1/2}h^{\zeta}\Hseminorm{\zeta}{v} + \betaparam\N^{1/2}h^{\xi}\Hseminorm{\xi}{\mu}
\lesssim
h^{\zeta-1/2}\Hseminorm{\zeta}{v} + h^{\xi+1/2}\Hseminorm{\xi}{\mu},
\end{equation*}

and so <ref> holds.

Motivated by the bounds on $\betaparam\D$ and $\betaparam\N$ in this proposition, we will later take
$\betaparam\D=\betaparam h^{-1}$
$\betaparam\N=\betaparam h$, where $\betaparam$ is a constant.

If $k=l$, $\betaparam\N\lesssim h^{-1}$, and the solution is smooth enough, then

\begin{equation*}
\betaparam\N^{1/2}h^{\xi}=
\betaparam\N^{1/2}h^{\zeta}\lesssim h^{\zeta-1/2}.
\end{equation*}

Therefore the same order of convergence will be observed when the bounds on $\betaparam\N$ here and in the
theorem below may be replaced by $\betaparam\N\lesssim h^{-1}$ without loss of convergence. In this case, both $\betaparam\N$ and $\betaparam\D$
may be taken to be constants independent of $h$.

We conclude that the best approximation result of <ref>
and the error estimate of <ref> hold for the discrete
solutions of <ref>, as given in the following theorem.

Let $(u,\lambda)\in H^{s}(\Gamma)\times H^{r}(\Gamma)$, for some $s\geqslant\tfrac12$ and $r\geqslant0$,
be the unique solution to the mixed Dirichlet–Neumann problem. This solution satisfies <ref>.
Let $(u_h,\lambda_h)\in\Vee_h^k\times\Lambda_h^l$ be the solution of <ref>.
If $0<\betaparam\D\lesssim h^{-1}$, $0<\betaparam\N\lesssim h$ and
$\exists\betamin>0$ such that $\betaparam\D^{1/2}\betaparam\N^{1/2}>\betamin$ and $\betaparam\D>\betamin$, then

\[
\Bnorm[\ND]{(u-u_h,\lambda-\lambda_h)} \lesssim
h^{\zeta-1/2} \Hseminorm{\zeta}{u} + h^{\xi+1/2}
\Hseminorm{\xi}{\lambda},
\]

where $\zeta = \min(k+1,s)$ and $\xi = \min(l+1,r)$.

If we set $\betaparam\D=0$ and $\betaparam\N=0$, we arrive at a penalty-free formulation for the mixed Dirichlet–Neumann problem.
We conjecture based on numerical experiments that this result also holds for the penalty-free formulation. The analysis for this case would take
a similar form as in the Dirichlet and Neumann penalty-free cases.

§.§ Application of the theory to the Robin problem

The formulation for Robin conditions was proposed in <ref>.
To simplify the notation we introduce a function $\robindenom:\Gamma \to \RR_+$ defined by

\[
\robindenom(\x) := \frac{1}{\param(\x) \betaparam\R(\x) + 1},
\]

and we assume that $\param$ and $\betaparam\R$ are sufficiently regular so that

\begin{equation}\label{eq:phi_smooth}
\robindenom \in W^{1,2}(\Gamma) \cap \LL^{\infty}(\Gamma).
\end{equation}

This will be true if the mesh has some local quasi-uniformity and
$\param$ is smooth enough. Noting that

\[
\robindenom-\tfrac12=\frac{2-(\param\betaparam\R+1)}{2(\param\betaparam\R+1)}=-\tfrac{1}{2}\frac{\param\betaparam\R-1}{\param\betaparam\R+1},
\]

we may then write the operators $\form{B}\R$ and $\form{L}\R$ as

\begin{align}\label{eq:BR}
\form{B}\R[(u,\lambda),(v,\mu)]&=
\left\langle (\robindenom-\tfrac12) u,\mu \right\rangle_{\Gamma\R}
-\left\langle (\robindenom-\tfrac12) \lambda, v \right\rangle_{\Gamma\R}
+\left\langle \robindenom \betaparam\R u, v \right\rangle_{\Gamma\R}
+\left\langle  \robindenom \param\lambda, \mu \right\rangle_{\Gamma\R},\\
\form{L}\R[(v,\mu)]&=
\left\langle(g\D +\param g\N)\robindenom,\betaparam\R v + \mu \right\rangle_{\Gamma\R}.
\end{align}

The boundary element method for the Robin problem reads: Find $(u_h,\lambda_h)
\in \productspace{V}_h$ such that

\begin{align}\label{eq:Robin_BEM}
\form{A}[(u_h,\lambda_h),(v_h,\mu_h)] +\form{B}\R[(u_h,\lambda_h),(v_h,\mu_h)]
\form{L}\R[(v_h,\mu_h)]
&&\forall (v_h,\mu_h) \in \productspace{V}_h.
\end{align}

For the analysis the following technical lemmas will be useful.

If $\varphi \in W^{1,2}(\Gamma) \cap \LL^{\infty}(\Gamma)$ and $f \in
H^{1/2}(\Gamma)$, then $\varphi f \in H^{1/2}(\Gamma)$ and

\[
\Hnorm{1/2}{\varphi f } \leqslant
C\left(\norm{\varphi}_{\LL^\infty(\Gamma)} + \norm{\varphi}_{W^{1,2}(\Gamma)}\right) \Hnorm{1/2}{f}.
\]

The proof is a consequence of <cit.> which shows that

\begin{equation}\label{Gagliardo_Nirenberg}
\Hnorm{1/2}{\varphi f } \leqslant
C\left(\norm{\varphi}_{\LL^\infty(\Gamma)} \Hnorm{1/2}{f} +
\norm{f}_{\LL^4(\Gamma)} \norm{\varphi}_{W^{1,2}(\Gamma)}^{1/2} \norm{\varphi}_{\LL^\infty(\Gamma)}^{1/2}\right).
\end{equation}

We then recall the Sobolev injection $\norm{f}_{\LL^4(\Gamma)} \leqslant C
\Hnorm{1/2}{f}$ from <cit.> and
conclude using this result and an arithmetic-geometric inequality of
the right hand side of <ref>.

If $\varphi,f\in\LL^{2}(\Gamma)$ and $\varphi(\x)>0$ for all $\x\in\Gamma$, then there exists $C>0$ such that

\[
\Ltwonorm{\varphi f}^2 \geqslant
\]

Let $a=\inf_{\x\in\Gamma}\varphi(\x)$. Since $\Gamma$ is closed, there exists $\y\in\Gamma$ such that $\varphi(\y)=a$.
Therefore $a>0$. We now see that

\begin{align*}
\Ltwonorm{\varphi f}^2&=\int_\Gamma \varphi^2f^2\\
&\geqslant a^2\int_\Gamma f^2\\
&= C\Ltwonorm{f}^2,
\end{align*}

where $C = a^2$.

We introduce the norm

\[
\Bnorm[\R]{(v,\mu)} :=\Vnorm{(v,\mu)} +
\Ltwonorm{(\param \robindenom)^{1/2}\mu}+\Ltwonorm{(\robindenom \betaparam\R)^{1/2} v}
\]

and set $\starnorm{\cdot}=\Bnorm[\R]{\cdot}$.
We note that if $\param\to0$ or $\param\to\infty$, then $\Bnorm[\R]\cdot$ converges to $\Bnorm[\D]\cdot$ or $\Bnorm[\N]\cdot$ respectively.
We now proceed to show that <ref> hold.

<Ref> are satisfied for the Robin problem.

Let $(v,\mu) \in \productspace{W}$, and let
$L:=\form{A}[(v,\mu),(v,\mu)] +\form{B}\R[(v,\mu),(v,\mu)]$.
Using <ref>, we see that

\begin{multline*}
\alpha\Hnorm{-1/2}{\mu}^2
\\+\Ltwonorm{(\param\robindenom)^{1/2}\mu}^2
\end{multline*}

for any $\alpha\leqslant\min(\alpha_\bopV,\alpha_\bopW)$.

By <ref>, we have

\begin{equation}
\end{equation}

Taking $\alpha=\min(\alpha_\bopV,\alpha_\bopW,C/2)$, we obtain

\begin{equation*}
\alpha\Hnorm{-1/2}{\mu}^2
\end{equation*}

Using the definition of $\Bnorm[\R]{\cdot}$, we see that the form is coercive.

<Ref> is satisfied for the Robin problem if $\exists\betamin>0$, independent of $h$, such that

Using <ref>, we see that for $g \in
H^{-1/2}(\Gamma)$, $\varphi \in W^{1,2}(\Gamma) \cap
\LL^\infty(\Gamma)$, and $f \in H^{1/2}(\Gamma)$,

\[
\left\langle  \robindenom g, f \right\rangle_\Gamma \leqslant  C\left(\norm{\varphi}_{\LL^\infty(\Gamma)} + \norm{\varphi}_{W^{1,2}(\Gamma)}\right) \Hnorm{-1/2}{g}\Hnorm{1/2}{f}.
\]

Let $\param_{\min}:=\inf_{\x\in\Gamma}\param(\x)$. As in the proof of <ref>, we see that $\param_{\min}>0$.

\begin{equation*}
\end{equation*}

and so

\begin{align*}
\norm{\robindenom-\tfrac12}_{\LL^\infty(\Gamma)} + \norm{\robindenom-\tfrac12}_{W^{1,2}(\Gamma)}
\max\left(\tfrac12,\frac1{\betamin\param_{\min}+1}\right)
\left(\norm{1}_{\LL^\infty(\Gamma)} + \norm{1}_{W^{1,2}(\Gamma)}\right).
\end{align*}

Applying these two results to the first two boundary terms in $\form{B}\R[(w,\eta),(v,\mu)]$, we obtain

\begin{align*}
\left\langle  (\robindenom-\tfrac12) w, \mu \right\rangle_\Gamma-\left\langle (\robindenom-\tfrac12) v,\eta \right\rangle_\Gamma
C\Vnorm{(w,\eta)} \Vnorm{( v,\mu)}.
\end{align*}

By the Cauchy–Schwarz inequality, we obtain for the remaining terms

\begin{multline*}
\left\langle \robindenom \param \eta, \mu \right\rangle_{\Gamma} +
\left\langle  \robindenom \betaparam\R w, v \right\rangle_{\Gamma}
\\\leqslant \Ltwonorm{(\robindenom \param)^{1/2} \eta}\Ltwonorm{(\robindenom\param)^{1/2} \mu}
+ \Ltwonorm{(\robindenom \betaparam\R)^{1/2} w} \Ltwonorm{(\robindenom \betaparam\R)^{1/2}  v}.
\end{multline*}

Collecting the terms, we then have

\[
\form{B}\R[(w,\eta),(v,\mu)] \lesssim \Bnorm[\R]{(w,\eta)} \Bnorm[\R]{(v,\mu)}.
\]

<Ref> is satisfied for the Robin problem if $\betaparam\R\eqsim h^{-1}$.

First note that $\robindenom<1$ and

\begin{align}\label{eq:robin_term_bounds}
\Ltwonorm{(\robindenom\betaparam\R)^{1/2}v}&\leqslant\betaparam\R^{1/2}\Ltwonorm{v}&&\text{and}&
\Ltwonorm{(\robindenom\param)^{1/2}\mu}&\leqslant\betaparam\R^{-1/2}\Ltwonorm{\mu}.
\end{align}

If $\betaparam\R\eqsim h^{-1}$, then <ref> can be shown to hold.

When using equal order approximation, the same order of convergence
will be observed when the bounds on $\betaparam\R$ here and in the theorem below may be replaced by $h\lesssim\betaparam\R\lesssim h^{-1}$
for sufficiently smooth solutions. Note that the condition $h^{-1}\lesssim\betaparam\R$ implies the existance of $\betamin$, as required
by <ref>. The condition $h\lesssim\betaparam\R$ does not imply this, so in this case the additional requirement that $\exists\betamin>0$
such that $\betamin<\betaparam\R$ is necessary to ensure continuity.

The extra assumptions in <ref> are satisfied for the Robin problem.

As a consequence of the coercivity and continuity above and observing
that by the Cauchy–Schwarz inequality and the definition of $\robindenom$,
there exists $C$ such that

\[
\left\langle \robindenom(g\D+ \param g\N), \betaparam\R v + \mu \right\rangle_{\Gamma} \leqslant
C(\Ltwonorm{g\D} + \Ltwonorm{g\N}) \Bnorm[\R]{(v,\mu)}
\]

We conclude that <ref> hold for the Robin problem.
This is summarised in the following result.

The Robin problem <ref> has a unique solution $(u,\lambda)\in H^{s}(\Gamma)\times H^{r}(\Gamma)$, for some $s\geqslant\tfrac12$ and $r\geqslant0$.
The discrete Robin problem <ref> is invertible. If $\betaparam\R\eqsim h^{-1}$, its solution $(u_h,\lambda_h)\in\Vee_h^k\times\Lambda_h^l$ satisfies

\[
\Bnorm[\R]{(u-u_h,\lambda-\lambda_h)} \leqslant C\left(
h^{\zeta-1/2} \Hseminorm{\zeta}{u} + h^{\xi+1/2}
\Hseminorm{\xi}{\lambda} \right),
\]

for some $C>0$,
where $\zeta = \min(k+1,s)$ and $\xi = \min(l+1,r)$. Additionally,

\[
\Hnorm[\Omega]{1}{\tilde u-\tilde u_h} \leqslant C\left(
h^{\zeta-1/2} \Hseminorm{\zeta}{u} + h^{\xi+1/2}
\Hseminorm{\xi}{\lambda} \right),
\]

where $\tilde u$ and $\tilde u_h$ are the solutions in $\Omega$ computed using <ref>.

Again, we could set $\betaparam\R=0$ to arrive at a penalty-free formulation for Robin problems.
In this case, our numerical experiments show large errors for some values of the parameter $\param$, which leads us to conclude that this result
does not hold for the penalty-free formulation.

As $\param\to0$ and $\param\to\infty$, we obtain the Dirichlet and Neumann formulations analysed in <ref>.
We expect the condition number of the discrete system for the Robin problem to be no worse than in either extreme case, and observe this in

§ NUMERICAL RESULTS

Drawing inspiration from [16], we define

\begin{align*}
u(x,y,z)&=\sin(\pi x)\sin(\pi y)\sinh(\sqrt2\pi z)\\
g\D(x,y,z)&=\sin(\pi x)\sin(\pi y)\sinh(\sqrt2\pi z),\\
\pi\cos(\pi x)\sin(\pi y)\sinh(\sqrt2\pi z)\\
\pi\sin(\pi x)\cos(\pi y)\sinh(\sqrt2\pi z)\\
\sqrt2\pi\sin(\pi x)\sin(\pi y)\cosh(\sqrt2\pi z)
\end{pmatrix}\cdot\B\nu.
\end{align*}

It is easy to check that for any bounded domain $\Omega\subset\RR^3$ with boundary $\Gamma=\Gamma\D\cup\Gamma\N\cup\Gamma\R$
and any fixed $\param\in\RR$,
$u$ is the solution of

\begin{align}
-\Delta u &=0       &&\text{in } \Omega,\\
u &= g\D              &&\text{on } \Gamma\D, \\
\frac{\partial u}{\partial\B\nu} &= g\N
&&\text{on } \Gamma\N,\\
\frac{\partial u}{\partial\B\nu} &= \frac{1}{\param} (u - g\D) + g\N
&&\text{on }\Gamma\R.
\end{align}

In the examples presented here, we let $\Omega$ be the unit sphere, and $\Gamma$ its boundary.
In the computations presented, a series of approximations of the sphere by plane triangles are used.
The results in this section were computed using the boundary element library Bempp [22], an open source boundary element library
developed by the authors of this paper.
All examples in this paper were computed with version 3.3.2 of the Bempp library.
Jupyter notebooks demonstrating the functionality used in this paper will be made available at <www.bempp.com>.

§.§ Dirichlet boundary conditions

First, we look at the case where $\Gamma=\Gamma\D$, in which the problem reduces to the Dirichlet problem:

\begin{align}
-\Delta u &=0       &&\text{in } \Omega,\\
u &= g\D              &&\text{on } \Gamma.
\end{align}

For this problem, we compare the penalty method proposed in this paper <ref> to the standard single layer formulation:
Find $\lambda\in \Lambda_h$ such that

\begin{align}\label{usual_Dir}
\left\langle\bopV \lambda,\mu\right\rangle&=
\left\langle(\tfrac12\bopI+\bopK) g\D,\mu\right\rangle
&&\forall\mu\in \Lambda_h.
\end{align}

[small,axis equal,axis on top,
xmin=0.011,xmax=1.9,xlabel=$h$,x dir=reverse,xmode=log,
ylabel near ticks,ylabel=,ymode=log
0.5 4.5
0.0221 0.0087890625
table img/data/linlin/diri1usual_convergence.dat;
table img/data/linlin/diri1_convergence.dat;

[small,axis on top,
xlabel=$h$,x dir=reverse,xmode=log,
ylabel near ticks,ylabel=Number of GMRES iterations,ymin=0,ylabel near ticks,ymax=100
table img/data/linlin/diri1usual_iterations.dat;
table img/data/linlin/diri1_iterations.dat;
table img/data/linlin/diri1nop_iterations.dat;

The convergence (left) and GMRES iteration counts (right) of the penalty method with $\betaparam\D=0.1$ ()
compared to the standard single layer method <ref> (),
for the Dirichlet problem on the unit sphere, with $k=l=1$.
The iteration count plot shows the number of iterations taken to solve the mass matrix preconditioned system ()
and the non-preconditioned system ().
The dashed line shows order 2 convergence.

<Ref> shows the convergence and iteration counts when $\betaparam\D=0.1$ and $k=l=1$, and so we look for
We note that as $h$ decreases, $h^{-1}$ increases, so $0.1\lesssim h^{-1}$.
In this case, $\Gamma$ is smooth, and so $\Vee_h^1=\tilde\Lambda_h^1$.
The iteration count plot (right) shows the number of iterations taken to solve the non-preconditioned system (), compared with the system
with mass matrix preconditioned applied blockwise from the left (), as described in [6].
Mass matrix preconditioning greatly reduces the number of iterations required,
so for the remainder of this paper, we precondition all linear systems using mass matrix preconditioning.

For larger and more complex geometries, however, more specialised preconditioners are required.
With systems of boundary element equations, it is common to use operator preconditioning or Calderón preconditioning [8], where
properties of the boundary operators at the continuous level are used to derive a preconditioned equation of a form known to be well conditioned.
In our case, it is not clear how to apply this approach, although further investigation of this warrants future work.

An alternative avenue of investigation leads to hierarchical LU based preconditioners, or even direct solvers of this type [4].
The penalty terms in this paper are all sparse matrices that have non-zero entries only for neighbouring triangles, and so adding these
terms only affects the entries in the matrix arising from near interactions; the far interactions—which are exactly those that are approximated
in a hierarchical matrix compression—are not affected by these terms. Therefore H-matrix methodst can be applied to this method
with few algorithmic changes required.

<Ref> shows the dependence of the error and iteration count on the chosen value of $\betaparam\D$, for a range of values of $h$.
It can be seen that the number of iterations increases when $\betaparam\D$ is above around 0.1, and the error increases when $\betaparam\D$ is
above 100. This motivates our earlier choice of $0.1$ as the value of $\betaparam\D$, although anything smaller than this
appears to be a good choice of $\betaparam\D$.

In <ref>, it can be seen that the penalty method proposed here gives comparable convergence to the standard method in a similar
number of iterations. However, the system in the penalty method contains around twice the number of unknowns, and so each iteration will be more

Additionally, the discrete systems for the penalty method are non-symmetric, so are solved using GMRES [21]. The discrete systems for the
standard method <ref> are symmetric, so CG [14] or MINRES [20] could be used: these methods are typically less expensive
than GMRES, so this is a further disadvantage of the penalty method for pure Dirichlet and Neumann problems and justifies our focus on more complex
boundary conditions.

[small,axis on top,
xlabel=$\betaparam\D$,xmin=1e-6*1.0001, xmax=1e6/1.0001,xmode=log,
ylabel near ticks,ylabel=,ymax=100,ymode=log
table img/data/linlin/diriA_convergence_4.dat;
table img/data/linlin/diriA_convergence_7.dat;
table img/data/linlin/diriA_convergence_10.dat;

[small,axis on top,
xlabel=$\betaparam\D$,xmin=1e-6*1.0001, xmax=1e6/1.0001,xmode=log,
ylabel near ticks,ylabel=Number of GMRES iterations,ymin=0,ymax=40
table img/data/linlin/diriA_iterations_4.dat;
table img/data/linlin/diriA_iterations_7.dat;
table img/data/linlin/diriA_iterations_10.dat;

The dependence of the error (left) and iteration count (right) on the value of $\betaparam\D$ for
$h=2^{-2}$ (),
$h=2^{-3.5}$ (), and
$h=2^{-5}$ (),
for the Dirichlet problem on the unit sphere, with $k=l=1$.

§.§ Mixed Dirichlet–Neumann boundary conditions

We now consider the case where $\Gamma=\Gamma\D\cup\Gamma\N$ and the problem reduces to a mixed Dirichlet–Neumann problem:

\begin{align}
-\Delta u &=0       &&\text{in } \Omega,\\
u &= g\D              &&\text{on } \Gamma\D, \\
\frac{\partial u}{\partial\B\nu} &= g\N
&&\text{on } \Gamma\N.
\end{align}

Let $\Gamma\N:=\{(x,y,z)\in\Gamma:x>0\}$ and $\Gamma\D:=\Gamma\setminus\Gamma\N$.
We use the same $g\D$ and $g\N$ as above.

We compare the method proposed in this paper with the standard method for mixed Dirichlet–Neumann problems <cit.>:
Find $(u,\lambda)\in \tilde{H}^{1/2}(\Gamma\N)\times \tilde{H}^{-1/2}(\Gamma\D)$ such that

\begin{multline}\label{classicalformulation2}
\left\langle\bopW_{\text{N}\text{N}}u,v\right\rangle
\\=
\\
\forall(v,\mu)\in \tilde{H}^{1/2}(\Gamma\N)\times \tilde{H}^{-1/2}(\Gamma\D),
\end{multline}

where for a given boundary operator $\bop{B}$, $\bop{B}_{ij}$ is the corresponding boundary operator with the integral taken over $\Gamma_i$
and the point $\x\in\Gamma_j$. For example, $\bop{V}_{\text{N}\text{D}}$ is defined by

\begin{align}
&&\text{for }\vec{x}\in\Gamma\D.
\end{align}

We first let $k=l+1=1$, and so look for
$(u_h,\lambda_h)\in\Vee_h^1\times\Lambda_h^0$. As motivated above by <ref>, we set
$\betaparam\D=\betaparam h^{-1}$
$\betaparam\N=\betaparam h$,
where $\betaparam$ is a constant.
The dependence of the error and iteration count on $\betaparam$ is shown in <ref>.
We observe that $\betaparam=0.01$ is a good choice, as this gives a small error and iteration count.

The convergence of the error as we reduce $h$ is shown in <ref>. Here we observe order 1.5 convergence, and
the same rate of convergence as the standard method <ref>, with a marginally lower error in the standard method.
The iteration count for the penalty method increases more gradually than the standard method, although this issue could be removed through better preconditioning of the standard

[small,axis on top,
xlabel=$\betaparam$,xmin=1e-6*1.0001, xmax=1e6/1.0001,xmode=log,
ylabel near ticks,ylabel=,ymax=100,ymode=log
table img/data/conlin/neudiriA_convergence_4.dat;
table img/data/conlin/neudiriA_convergence_7.dat;
table img/data/conlin/neudiriA_convergence_10.dat;

[small,axis on top,
xlabel=$\betaparam$,xmin=1e-6*1.0001, xmax=1e6/1.0001,xmode=log,
ylabel near ticks,ylabel=Number of GMRES iterations,ymin=0,ymax=350
table img/data/conlin/neudiriA_iterations_4.dat;
table img/data/conlin/neudiriA_iterations_7.dat;
table img/data/conlin/neudiriA_iterations_10.dat;

The dependence of the error (left) and iteration count (right) on the value of $\betaparam$ for
$h=2^{-2}$ (),
$h=2^{-3.5}$ (), and
$h=2^{-5}$ (),
for the mixed Dirichlet–Neumann problem on the unit sphere, with $k=l+1=1$.
Here we use $\betaparam\D=\betaparam h^{-1}$ and $\betaparam\N=\betaparam h$.

[small,axis equal,axis on top,
xlabel=$h$,xmode=log,x dir=reverse,
ylabel near ticks,ylabel=,ymode=log
0.4 32.925
0.02 0.3681126907959028
table img/mixed/conlin/hyp_convergence.dat;
table img/data/conlin/neudiri1_convergence.dat;

[small,axis on top,
xlabel=$h$,xmode=log,x dir=reverse,
ylabel near ticks,ylabel=Number of GMRES iterations,ymin=0,ymax=350
table img/mixed/conlin/hyp_iterations.dat;
table img/data/conlin/neudiri1_iterations.dat;

The convergence (left) and iterations counts (right) of the penalty method with $\betaparam=0.01$ ()
compared to the standard method <ref> (),
for the mixed Dirichlet–Neumann problem on the unit sphere, with $k=l+1=1$.
The dashed line shows order 1.5 convergence.
Here we use $\betaparam\D=\betaparam h^{-1}$ and $\betaparam\N=\betaparam h$.

We next consider the case where $k=l=1$. In this case, as remarked in <ref>, we may replace
the bound on $\betaparam\N$ by $\betaparam\N\lesssim h^{-1}$, and so we may take both $\betaparam\D$ and $\betaparam\N$ to be constant:
we set $\betaparam\D=\betaparam\N=\betaparam$. The dependence of the error and iteration count on $\betaparam$ for this choice of
parameters is shown in <ref>.

The convergence to the solution when $\betaparam=0.01$ is shown in <ref>.
It can be seen here that order 2 convergence is observed, higher than the expected order 1.5 convergence.
In this case, the standard method <ref> only achieves order 1 convergence, with a much higher iteration
count that the penalty method. For this choice of discrete spaces, we also compared the our method with the formulation
given in <cit.>: this formulation is better conditioned than <ref> but still achieves
only order 1 convergence.

[small,axis on top,
xlabel=$\betaparam$,xmin=1e-6*1.0001, xmax=1e6/1.0001,xmode=log,
ylabel near ticks,ylabel=,ymax=10,ymode=log
table img/data/linlin/neudiriA_convergence_4.dat;
table img/data/linlin/neudiriA_convergence_7.dat;
table img/data/linlin/neudiriA_convergence_10.dat;

[small,axis on top,
xlabel=$\betaparam$,xmin=1e-6*1.0001, xmax=1e6/1.0001,xmode=log,
ylabel near ticks,ylabel=Number of GMRES iterations,ymin=0,ymax=100
table img/data/linlin/neudiriA_iterations_4.dat;
table img/data/linlin/neudiriA_iterations_7.dat;
table img/data/linlin/neudiriA_iterations_10.dat;

The dependence of the error (left) and iteration count (right) on the value of $\betaparam$ for
$h=2^{-2}$ (),
$h=2^{-3.5}$ (), and
$h=2^{-5}$ (),
for the mixed Dirichlet–Neumann problem on the unit sphere, with $k=l=1$.
Here we use $\betaparam\D=\betaparam\N=\betaparam$.

[small,axis equal,axis on top,width=63mm,
xlabel=$h$,xmode=log,x dir=reverse,
ylabel near ticks,ylabel=,ymode=log
0.4 20
0.02 1
0.35355 4.585
0.02 0.014672281414357528
table img/mixed/linlin/hyp_convergence.dat;
table img/mixed/linlin/noh_convergence.dat;
table img/data/linlin/neudiri1_convergence.dat;

[small,axis on top,width=63mm,
xlabel=$h$,xmode=log,x dir=reverse,
ylabel near ticks,ylabel=Number of GMRES iterations,ymin=0,ymax=100
table img/mixed/linlin/hyp_iterations.dat;
table img/mixed/linlin/noh_iterations.dat;
table img/data/linlin/neudiri1_iterations.dat;

The convergence (left) and iterations counts (right) of the penalty method with $\betaparam=0.01$ ()
compared to the standard method <ref> () and the method given in <cit.> (),
for the mixed Dirichlet–Neumann problem on the unit sphere, with $k=l=1$.
The dashed lines show order 2 and order 1 convergence.
Here we use $\betaparam\D=\betaparam\N=\betaparam$.

In <ref>, the error and iteration count remain steady as $\betaparam\to0$.
In numerical experiments on a sphere and cube with $\betaparam=0$, we see similar convergence to that observed in this section.
This leads us to conjecture that <ref> will hold for the penalty-free formulation, when $\betaparam=0$.

§.§ Robin problem

We now consider the case where $\Gamma=\Gamma\R$ and the problem reduces to a Robin problem:

\begin{align}
-\Delta u &=0       &&\text{in } \Omega,\\
\frac{\partial u}{\partial\B\nu} &= \frac{1}{\param} (u - g\D) + g\N
&&\text{on }\Gamma,
\end{align}

for some $\param\in\RR$.

In this section, we compare the method proposed in this paper with the standard method: Find $\lambda\in H^{-1/2}(\Gamma)$ such that

\begin{align}\label{classicalRobin}
\left\langle\bopW u,v\right\rangle+\left\langle\frac1\param\left(\tfrac12\bopI-\bopKadj\right)u,v\right\rangle
&=\left\langle\left(\tfrac12\bopI-\bopKadj\right)\left(\frac1\param g\D+g\N\right),v\right\rangle&&\forall\mu\in H^{-1/2}(\Gamma).
\end{align}

table img/data/conlin/robin3d_error.dat;

zlabel=Number of GMRES iterations,zmin=0,zmax=200
table img/data/conlin/robin3d_iterations_capped.dat;

The dependence of the error on $\param$ and $\betaparam$
for the Robin problem on the unit sphere with $h=0.1$, with $k=l+1=1$.
Here we use $\betaparam\D=\betaparam h^{-1}$ and $\betaparam\N=\betaparam h$.

Again, we begin letting $k=l+1=1$.
Here we use

\begin{equation*}
\betaparam\R:=\frac{\param\betaparam\N+\betaparam\D}{\param+1},
\end{equation*}

where $\betaparam\D=\betaparam h^{-1}$
and $\betaparam\N=\betaparam h$, for some constant $\betaparam$, as in the mixed Dirichlet–Neumann case.

The dependence of the error and iteration count on both $\param$ and $\betaparam$, on a grid with $h=0.1$, is shown in <ref>.
The convergence as $h$ is reduced
for $\param=\frac1{300}$, $\param=1$, and $\param=300$, and using $\betaparam=0.01$,
is shown in <ref>. In this case, order 1.5 convergence is observed.

[small,axis equal,axis on top,
xlabel=$h$,xmode=log,x dir=reverse,
ylabel near ticks,ylabel=,ymode=log,
0.4 32.925
0.02 0.3681126907959028
table img/data/conlin/robin1_convergence_DD.dat;
table img/data/conlin/robin1_convergence_DN.dat;
table img/data/conlin/robin1_convergence_NN.dat;

[small,axis on top,
xlabel=$h$,x dir=reverse,xmode=log,
ylabel near ticks,ylabel=Number of GMRES iterations,ymin=0,ymax=100
table img/data/conlin/robin1_iterations_DD.dat;
table img/data/conlin/robin1_iterations_DN.dat;
table img/data/conlin/robin1_iterations_NN.dat;

The convergence (left) and iteration counts (right) of the penalty method for the Robin problem with
$\param=300$ (),
$\param=1$ () and
$\param=1/300$ ()
on the unit sphere, using $k=l+1=1$ and $\betaparam=0.01$.
The dashed line shows order 1.5 convergence.
Here we use $\betaparam\D=\betaparam h^{-1}$ and $\betaparam\N=\betaparam h$.

As in the mixed Dirichlet–Neumann case, when $k=l=1$, we may replace the bound on $\betaparam\N$ with $\betaparam\N\lesssim h^{-1}$.
Again, we take $\betaparam\D=\betaparam\N=\betaparam$ for some constant $\betaparam$. The dependence of the error
and iteration count on both $\betaparam$ and $\param$ is shown in <ref>. As in the previous case, $\betaparam=0.01$ looks to
be a suitable choice for the parameter.

The convergence as we reduce $h$
for $\param=\frac1{300}$, $\param=1$, and $\param=300$, and using $\betaparam=0.01$,
is shown in <ref>. In this case, order 2 convergence is observed. For the standard method <ref>, the
same order of convergence and errors of almost identical size are oberved. For the standard method, the number of iterations required to
solve the system is higher for smaller values of $\param$; for the penalty method, the number of iterations is less affected by the value of
$\param$, leading to lower iteration counts than the standard method for small values of $\param$.

table img/data/linlin/robin3d_error.dat;

zlabel=Number of GMRES iterations,zmin=0,zmax=200
table img/data/linlin/robin3d_iterations_capped.dat;

The dependence of the error on $\param$ and $\betaparam$
for the Robin problem on the unit sphere with $h=0.1$, with $k=l=1$.
Here we use $\betaparam\D=\betaparam\N=\betaparam$.

[small,axis equal,axis on top,
xlabel=$h$,xmode=log,x dir=reverse,
ylabel near ticks,ylabel=,ymode=log,
0.35355 4
0.02 0.01280024550870886
0.707106781187 2.87013094897
0.5 1.609804802
0.353553390593 0.815730862548
0.25 0.409036605736
0.176776695297 0.182266780281
0.125 0.0741445440987
0.0883883476483 0.0397988294805
0.0625 0.018207955525
0.0441941738242 0.00817243204977
0.03125 0.00563246334918
0.707106781187 3.08850637497
0.5 1.71174947496
0.353553390593 0.866335489945
0.25 0.428882478169
0.176776695297 0.190517185996
0.125 0.0766963763636
0.0883883476483 0.0376485201851
0.0625 0.0169008482716
0.0441941738242 0.00799559314468
0.03125 0.00474823694878
0.707106781187 3.13459257777
0.5 1.72374221257
0.353553390593 0.872095478648
0.25 0.429629000915
0.176776695297 0.190529491281
0.125 0.0762845891653
0.0883883476483 0.037174727955
0.0625 0.0165910353303
0.0441941738242 0.00783462775204
0.03125 0.00457475941363

table img/data/linlin/robin1_convergence_DD.dat;
table img/data/linlin/robin1_convergence_DN.dat;
table img/data/linlin/robin1_convergence_NN.dat;

[small,axis on top,
xlabel=$h$,x dir=reverse,xmode=log,
ylabel near ticks,ylabel=Number of GMRES iterations,ymin=0,ymax=40
0.707106781187 5
0.5 5
0.353553390593 4
0.25 4
0.176776695297 4
0.125 4
0.0883883476483 3
0.0625 4
0.0441941738242 4
0.03125 4
0.707106781187 10
0.5 11
0.353553390593 9
0.25 11
0.176776695297 12
0.125 13
0.0883883476483 14
0.0625 17
0.0441941738242 20
0.03125 24
0.707106781187 19
0.5 20
0.353553390593 23
0.25 13
0.176776695297 13
0.125 15
0.0883883476483 16
0.0625 18
0.0441941738242 21
0.03125 26

table img/data/linlin/robin1_iterations_DD.dat;
table img/data/linlin/robin1_iterations_DN.dat;
table img/data/linlin/robin1_iterations_NN.dat;

The convergence (left) and iteration counts (right) of the penalty method ()
compared to the standard method <ref> (),
for the Robin problem with
$\param=300$ (),
$\param=1$ () and
$\param=1/300$ ()
on the unit sphere, using $k=l=1$ and $\betaparam=0.01$.
The dashed line shows order 2 convergence.
Here we use $\betaparam\D=\betaparam\N=\betaparam$.

Again, we could consider the penalty-free formulation for the Robin problem. However, <ref> suggest that as
$\betaparam\to0$, the error increases for some values of $\param$. This increased error can also be observed in the numerical
experiments we have run with $\betaparam=0$. Hence in the Robin case, the penalty term is necessary and <ref> does not hold
for $\betaparam\R=0$.

§ CONCLUSIONS
We have analysed and demonstrated the effectiveness of Nitsche type coupling methods for
boundary element formulations. In particular, for Robin and mixed Neumann/Dirichlet boundary conditions these
are simpler than the strong imposition of boundary conditions since the boundary condition only enters the
equations through a sparse operator.

An open problem is preconditioning. While the iteration counts in the presented examples were already
practically useful, for large and complex structures preconditioning is still essential. The hope is to use
the properties of the Calderón projector to build effective operator preconditioning techniques for the
presented Nitsche type frameworks.

An extension of the presented method to FEM/BEM formulations is currently in preparation. Other directions
are the Helmholtz and Maxwell problems. Although the analysis for these cases is more involved, we expect that
their implementation will be structurally similar to the presented Laplace case.

Acknowledgements. Erik Burman and Timo Betcke were supported by EPSRC
grants EP/P01576X/1 and EP/P012434/1.

[1]
I. Babus̆ka, The finite element method with penalty, Mathematics
of Computation, 27 (1973), pp. 221–228.

[2]
I. Babuška and A. K. Aziz, Survey lectures on the mathematical
foundations of the finite element method, in The mathematical foundations of
the finite element method with applications to partial differential equations
(Proc. Sympos., Univ. Maryland, Baltimore, Md., 1972), Academic
Press, New York, 1972, pp. 1–359.
With the collaboration of G. Fix and R. B. Kellogg.

[3]
L. Baratchart, L. Bourgeois, and J. Leblond, Uniqueness results for
inverse Robin problems with bounded coefficient, Journal of Functional
Analysis, 270 (2016), pp. 2508–2542.

[4]
M. Bebendorf, Hierarchical LU decomposition-based preconditioners
for BEM, Computing, 74 (2005), pp. 225–247.

[5]
R. Becker, P. Hansbo, and R. Stenberg, A finite element method for
domain decomposition with non-matching grids, M2AN. Mathematical Modelling
and Numerical Analysis, 37 (2003), pp. 209–225.

[6]
T. Betcke, M. W. Scroggs, and W. Śmigaj, Product algebras for
Galerkin discretizations of boundary integral operators and their
in preparation, 2017.

[7]
H. Brezis and P. Mironescu, Gagliardo-Nirenberg, composition and
products in fractional Sobolev spaces, J. Evol. Equ., 1 (2001),
pp. 387–404.
Dedicated to the memory of Tosio Kato.

[8]
A. Buffa and S. H. Christiansen, A dual finite element complex on
the barycentric refinement, Comptes Rendus Mathematique, 340 (2005),
pp. 461–464.

[9]
A. Chernov and P. Hansbo, Nonconforming hp-FE/BE coupling for
interface problems based on nitsche's method.

[10]
F. Chouly and N. Heuer, A Nitsche-based domain decomposition
method for hypersingular integral equations, Numerische Mathematik, 121
(2012), pp. 705–729.

[11]
M. Costabel and E. Stephan, Boundary integral equations for mixed
boundary value problems in polygonal domains and Galerkin approximation,
Banach Center Publications, 15 (1985), pp. 175–251.

[12]
E. Di Nezza, G. Palatucci, and E. Valdinoci, Hitchhiker's guide to
the fractional Sobolev spaces, Bulletin des Sciences Mathématiques, 136
(2012), pp. 521–573.

[13]
G. Gatica, M. Healey, and N. Heuer, The boundary element method with
Lagrangian multipliers, Numerical Methods for Partial Differential
Equations, 25 (2009), pp. 1303–1319.

[14]
M. R. Hestenes and E. Stiefel, Methods of conjugate gradients for
solving linear systems, Journal of research of the National Bureau of
Standards, 49 (1952), pp. 409–436.

[15]
B. Jin and J. Zou, Numerical estimation of the Robin coefficient
in a stationary diffusion equation, IMA Journal of Numerical Analysis, 30
(2010), pp. 677–701.

[16]
M. Juntunen and R. Stenberg, Nitsche's method for general boundary
conditions, Math. Comp., 78 (2009), pp. 1353–1374.

[17]
U. Langer and O. Steinbach, Boundary element tearing and
interconnecting methods, Computing. Archives for Scientific Computing, 71
(2003), pp. 205–228.

[18]
U. Langer and O. Steinbach, Coupled boundary and finite element
tearing and interconnecting methods, in Domain decomposition methods in
science and engineering, vol. 40 of Lecture Notes in Computational Science
and Engineering, Springer, Berlin, 2005, pp. 83–97.

[19]
J. Nitsche, Über ein Variationsprinzip zur Lösung von
Dirichlet-Problemen bei Verwendung von Teilräumen, die keinen
Randbedingungen unterworfen sind, Abhandlungen aus dem Mathematischen
Seminar der Universität Hamburg, 36 (1971), pp. 9–15.
Collection of articles dedicated to Lothar Collatz on his sixtieth

[20]
C. Paige and M. Saunders, Solution of sparse indefinite systems of
linear equations, SIAM Journal on Numerical Analysis, 12 (1975),
pp. 617–629.

[21]
Y. Saad and M. H. Schultz, GMRES: A generalized minimal residual
algorithm for solving nonsymmetric linear systems, SIAM Journal on
Scientific and Statistical Computing, 7 (1986), pp. 856–869.

[22]
W. Śmigaj, T. Betcke, S. Arridge, J. Phillips, and M. Schweiger, Solving boundary integral problems with BEM++, ACM Transactions on
Mathematical Software, 41 (2015), pp. Art. 6, 40.

[23]
O. Steinbach, Numerical approximation methods for elliptic boundary
value problems, Springer, New York, 2008.
Finite and boundary elements, Translated from the 2003 German

[24]
O. Steinbach and M. Windisch, Stable boundary element domain
decomposition methods for the Helmholtz equation, Numerische Mathematik,
118 (2011), pp. 171–195.

[25]
E. Stephan, Boundary integral equations for mixed boundary value
problems in $\mathbb{R}^3$, Mathematische Nachrichten, 134 (1987),
pp. 21–53.

[26]
T. Von Petersdorff and R. Leis, Boundary integral equations for
mixed Dirichlet, Neumann and transmission problems, Mathematical Methods
in the Applied Sciences, 11 (1989), pp. 185–213.

[27]
T. Von Petersdorff and E. P. Stephan, Regularity of mixed boundary
value problems in $\mathbb{R}^3$ and boundary element methods on graded
meshes, Mathematical Methods in the Applied Sciences, 12 (1990),
pp. 229–249.

[28]
J. Xu and L. Zikatanov, Some observations on Babuška and
Brezzi theories, Numerische Mathematik, 94 (2003), pp. 195–202.