#  IgA-BEM for 3D Helmholtz problems on multi-patch domains using B-spline
tailored numerical integration

Antonella Falini  Tadej Kanduč<EMAIL_ADDRESS>Maria Lucia Sampoli
Alessandra Sestini  Dipartimento di Informatica, Università degli Studi di
Bari Aldo Moro,
Via Orabona 4, Bari, Italy  Faculty of Mathematics and Physics, University of
Ljubljana,
Jadranska ulica 19,1000 Ljubljana, Slovenia  Department of Information
Engineering and Mathematics, University of Siena,
Via Roma 56, Siena, Italy  Department of Mathematics and Computer Science,
University of Firenze,
Viale Morgagni 67, Firenze, Italy

###### Abstract

An Isogeometric Boundary Element Method (IgA-BEM) is considered for the
numerical solution of Helmholtz problems on 3D bounded or unbounded domains,
admitting a smooth conformal multi-patch representation of their finite
boundary surface. The discretization space is formed by $C^{0}$ inter-patch
continuous basis functions whose restriction to a patch simplifies to the span
of tensor product B-splines composed with the given patch parameterization.
For both regular and singular integration, the proposed model utilizes a
numerical procedure defined on the support of each trial B-spline function,
which makes possible a function–by–function implementation of the matrix
assembly phase. Spline quasi-interpolation is the common ingredient of all the
considered quadrature rules; in the singular case it is combined with a
B-spline recursion over the spline degree and with a singularity extraction
technique, extended to the multi-patch setting for the first time. A threshold
selection strategy is proposed to automatically distinguish between nearly
singular and regular integrals. Numerical examples on relevant benchmarks show
that the expected convergence orders are achieved with uniform discretization
and a small number of uniformly spaced quadrature nodes.

###### keywords:

Helmholtz equation, Isogeometric Analysis (IgA), Boundary Element Method
(BEM), singular and nearly singular integral, numerical integration, B-spline
quasi-interpolation.

##  1 Introduction

Isogeometric Analysis (IgA) [1] is a powerful tool to obtain a numerical
solution of problems governed by partial differential equations, introduced in
the literature at the beginning of the new millennium. This new paradigm was
motivated by the observation that in engineering applications a domain is
described by its boundary parametric representation generated by Computer
Aided Design (CAD) software. For this aim CAD relies on flexible forms, often
based on multi-patch tensor product B-spline or rational B-spline (NURBS)
spaces, suited also for applications [2]. The IgA idea consists in adopting
the CAD functional spaces also for approximating the solution of the
differential problem, taking into account the well-known optimal approximation
power of spline spaces [3, 4]. The increasing success of the IgA paradigm in
the context of both domain and boundary element methods is easily explained
considering that, besides being capable to keep an exact representation of
complex domains described in a multi-patch CAD form, IgA formulations
guarantee a certain level of accuracy with considerably less degrees of
freedom than traditional Finite Element Analysis (FEA), which relies on larger
piecewise polynomial spaces with lower inter-element smoothness within each
patch [5]. The attractiveness of IgA is even higher for Boundary Element
Methods (BEMs), see for example [6, 7, 8, 9, 10, 11] and references therein,
because, in contrast to domain methods like FEM, they do not require any
preliminary volumetric parameterization, which can be a remarkably time
consuming task. BEMs rely on a boundary integral formulation of the problem,
which can be derived whenever the fundamental solution of the associated
differential operator is known. Hence, they require only the definition of a
mesh on the boundary of the domain, which is a much easier task especially in
the IgA context, where the CAD representation of the boundary is available.

The precursor of this work, dealing with 3D IgA-BEM, is presented in [12],
where exterior Laplace problems with Dirichlet boundary conditions on single
patch domains are considered and, in order to deal also with screen problems,
an indirect boundary integral formulation is adopted. In this paper the IgA
multi-patch formulation of BEMs (multi-patch IgA-BEMs) based on a direct
boundary integral formulation is considered for the numerical solution of
problems governed by the Helmholtz equation on 3D bounded or unbounded
domains, equipped with either Dirichlet or Neumann boundary conditions. Such
problems are of interest in acoustics, since they model radiation and rigid
scattering in the frequency domain [13], where radiation refers to the
pressure field produced by an object vibrating within a fluid (usually air or
water), while rigid scattering relates to the disturbance caused by an
obstacle immersed within an existing acoustic field. The adopted Boundary
Integral Equation (BIE) discretized with a collocation approach is the so-
called conventional BIE (CBIE), which involves only weakly singular integrals
and is augmented with the Sommerfield radiation condition at infinity for
exterior problems. The Helmholtz operator can thus be simply written as a sum
of the Laplacian $\Delta$ with an identity operator scaled with a positive
parameter $\kappa^{2}$. However, even when dealing with exterior Helmholtz
problems admitting a unique solution, such BIE generates spurious solutions
when $\kappa^{2}$ is an eigenvalue of $-\Delta$ for the considered boundary
conditions. This is the reason why other papers also relying on IgA-BEMs for
the numerical solution of the Helmholtz equation introduce its alternative
Burton-Miller (BM) boundary integral formulation, which does not suffer from
such drawback. Obtained as a linear combination of the mentioned CBIE with the
hypersingular BIE (HBIE) also paired with the Helmholtz operator, the BM
integral formulation is clearly more difficult, see for example [14], where it
is combined with a regularization technique to avoid to deal with
hypersingular integrals in the assembly phase. However, successive researches
have given numerical evidence that, when $\kappa^{2}$ is not an eigenvalue of
$-\Delta,$ better numerical results are achieved by using the simpler CBIE,
see for example [15, 16]. Since clearly only a finite number of different
values of the parameter $\kappa^{2}$ can be considered and this can be anyway
sufficient for a reasonable reconstruction of an associated time-dependent
acoustic field, we have opted for applying the IgA-BEM approach just to the
CBIE associated with Helmholtz, with the implicit assumption to avoid values
of $\kappa^{2}$ equal to eigenvalues of $-\Delta.$

The IgA-BEM implementation in this study relies on the B-spline tailored
cubature formulas for regular and weakly singular integrals, introduced in
[17] in the setting of weakly singular surface integrals. The theoretical
convergence order of such rules is analyzed in this work for the first time.
Since the proposed quadrature schemes are constructed on the support of every
bivariate B-spline, a function-by-function implementation of the matrix
assembly phase can be adopted, similarly to work for 2D Laplace problems in
[18]. This attractive feature markedly reduces the amount of data accessing in
the assembly phase, compared to a more common element–by–element strategy.
Another remarkable feature of the proposed quadratures is the freedom of
choice for the quadrature nodes – no reparameterization of nodes near the
singularities is needed and, by choosing uniformly spaced nodes in the
parametric space for all the required regular and singular numerical
integrations, it ensures the possibility of using a small overall number of
unique nodes. By applying a singularity subtraction technique on singular
integrals, the singular and regular part of the integrand are effectively
separated. Thanks to the thorough study developed in [19] for the weakly
singular kernel expansions, the continuity of the regularized integrands can
be controlled and the developed integration rules have consequently better
approximation properties for general smooth geometry parameterizations. A
novelty in this work is a careful decoupling of integrals appearing in the
Helmholtz BIE into simpler, real and imaginary, regular and singular
integrals. An original strategy for the automatic detection of near
singularity is also introduced. Another new key feature lies in coupling the
integration routines with multi-patch geometry representations; when the
source point and the integration domain are separated by a patch interface,
the unmatching parameterization across it poses an additional challenge. To
overcome difficulties of this type, for nearly singular inter-patch integrals
a new approximated singularity subtraction technique is introduced. It
replaces the source point with its projection on an extrapolation of the patch
surface, where the integration has to be performed.

The paper is organized as follows. The next section introduces the Helmholtz
problem and its conventional boundary integral formulation for both Dirichlet
and Neumann boundary conditions. Then, in Section 3 the multi-patch IgA
setting for the standard CAD representation of a bounded or unbounded domain
is described, as well as the finite dimensional functional spaces adopted for
the analysis. Section 4 describes the developed quadrature formulas, based on
the tensor-product formulation of a discrete spline quasi-interpolation
approach. The formulas are combined with a subtraction regularization
technique and applied to the single and double layer Helmholtz kernels. In
Section 5 the results obtained for some benchmark problems are reported and
commented, while some conclusive remarks are given in Section 6. Finally, some
technical details related to considered singular kernels are available in
Appendix.

##  2 The Helmholtz problem and its boundary integral formulation

We study 3D potential problems described by the Helmholtz equation with
Dirichlet or Neumann boundary conditions on domains
$\Omega\subset{\mathop{{\rm I}\kern-1.99997pt{\rm R}}\nolimits}^{3}$ admitting
a connected closed boundary surface $\Gamma=\partial\Omega,$

$\displaystyle\left\\{\begin{array}[]{ll}\Delta u+\kappa^{2}u=0&{\rm
in}\;\Omega,\\\ u=u_{\rm D}\quad{\rm or}\;\quad\frac{\partial u}{\partial{\bf
n}}=u_{\rm N}&{\rm on}\;\Gamma,\end{array}\right.$ (3)

where $u:\Omega\rightarrow\mathbb{C}$ denotes the unknown potential, ${\bf n}$
the unit normal on $\Gamma$ pointing outward from $\Omega$ and $\kappa>0.$
When the domain is a finite volume the problem is interior and we use the
notation $\Omega^{[i]}$ to underline the type of this domain when it is
useful. In the opposite case, we set $\Omega^{[e]}$ and we study the exterior
problem. Concerning the regularity of the boundary surface $\Gamma,$ for
simplicity we assume to be smooth, at least without self-intersections and
with a tangent plane well defined at each point and varying continuously. Note
that for exterior problems, equation (3) has to be augmented with an
additional condition, the so called Sommerfeld radiation condition,

$\displaystyle\frac{\partial u}{\partial r}-i\kappa
u=o\left(\frac{1}{r}\right),$ (4)

where $r$ denotes the point distance from the origin of the reference system
and $i$ is the imaginary unit; see for example [13]. This condition at
infinity is necessary to ensure for any positive $\kappa$ the existence and
uniqueness of the solution $u$ for exterior Helmholtz problems for both
Dirichlet and Neumann boundary conditions. Uniqueness of $u$ for interior
problems is ensured when $\kappa^{2}$ is not an eigenvalue of the reversed
Laplacian operator $-\Delta$ in $\Omega$ 111The scalar $\lambda$ is an
eigenvalue for $-\Delta$ if there exists a non vanishing function
$u_{\lambda}$ such that $-\Delta u_{\lambda}=\lambda u_{\lambda}$, also
fulfilling homogeneous boundary condition of the assigned type. It is possible
to prove that the operator $-\Delta$ admits only positive eigenvalues, which
define an unbounded infinite sequence of positive numbers depending on the
considered kind of boundary conditions and on the shape of the finite domain
$\Omega,$ see for example [20].. Regarding the regularity of the weak solution
$u$ of (3), we observe that $u$ belongs to the Sobolev space $H^{1}(\Omega),$
provided that, for Dirichlet (Neumann) boundary conditions, $u_{\rm D}\in
H^{1/2}(\Gamma)\,,\,\left(u_{\rm N}\in H^{-1/2}(\Gamma)\right)$.
The Helmholtz equation is of particular interest in acoustic because the
solution $u=u({\bf x},\kappa)$, ${\bf x}\in\Omega$ of (3) can be interpreted
as the inverse Fourier transform [21] of the time-dependent scattered pressure
field $p=p({\bf x}\,,\,t)$ generated by the scatterer in a given fluid domain,
which fulfills the wave equation,

$\frac{1}{c_{f}^{2}}\ \frac{\partial^{2}p}{\partial t^{2}}=\Delta p.$ (5)

Namely, setting

$\breve{p}=\breve{p}({\bf x}\,,\,\nu):=\int_{-\infty}^{+\infty}p({\bf
x}\,,\,t)e^{i\nu t}\,dt$

it is easy to check that $\breve{p}$ verifies the Helmholtz equation with
$\kappa=\nu/c_{f}$ and that, conversely, it is

$p=\frac{1}{2\pi}\int_{-\infty}^{+\infty}\breve{p}({\bf x}\,,\,\nu)e^{i\nu
t}\,d\nu\,.$

This is of interest because it allows us to deal with a differential problem
formulated just in terms of spatial variables.

In order to approximate the solution of (3) with IgA-BEMs, we consider the
Conventional Boundary Integral Equation (CBIE) associated with the Helmholtz
equation and derived by direct approach,

$\displaystyle\int_{\Gamma}{{\mathcal{G}}}_{\kappa}({\bf x},{\bf
y})\frac{\partial u}{\partial{\bf n}}({\bf y})d\Gamma_{{\bf y}}=c({\bf
x})u({\bf
x})+\int_{\Gamma}\frac{\partial{{\mathcal{G}}}_{\kappa}}{\partial{\bf n}_{{\bf
y}}}({\bf x},{\bf y})\ u({\bf y})\ d\Gamma_{{\bf y}},\qquad{\bf x}\in\Gamma,$
(6)

where $r:=\|{\bf r}\|_{2}$ with ${\bf r}:={\bf x}-{\bf y}$. The kernel
${{\mathcal{G}}}_{\kappa}$ and its derivative with respect to the exterior
unit normal ${\bf n}$ to $\Gamma$ computed in the point ${\bf y}\in\Gamma$ are
defined as

$\displaystyle{{\mathcal{G}}}_{\kappa}({\bf x},{\bf y}):=\frac{1}{4\pi
r}e^{{\rm i}\kappa
r},\qquad\qquad\frac{\partial{{\mathcal{G}}}_{\kappa}}{\partial{\bf n}_{{\bf
y}}}({\bf x},{\bf y})\,=\,\frac{1}{4\pi r}e^{{\rm i}\kappa
r}(-\frac{1}{r}+{\rm i}\kappa)\ \frac{\partial r}{\partial{\bf n}_{{\bf
y}}}\,,$

with

$\displaystyle\frac{\partial r}{\partial{\bf n}_{{\bf y}}}=-\frac{{\bf
r}\cdot{\bf n}_{{\bf y}}}{r},\qquad\qquad c({\bf
x}):=-\int_{\Gamma}\frac{\partial{\mathcal{G}}_{0}}{\partial{\bf n}_{{\bf
y}}}({\bf x},{\bf y})d\Gamma_{{\bf y}}\,,$

being $c({\bf x})\equiv 1/2$ for the assumed smooth surface $\Gamma$. The
restrictions of $u$ and $({\partial u}/{\partial{\bf n}})({\bf y})$ to
$\Gamma$ are usually called Cauchy data, being one of them coincident with the
available boundary datum and the other the unknown. By separating the real and
imaginary parts of ${\mathcal{G}}_{\kappa},$ we obtain

$\displaystyle{\mathcal{G}}_{\kappa}=\frac{1}{4\pi}\left(\frac{\cos(\kappa
r)}{r}+{\rm i}\frac{\sin(\kappa r)}{r}\right).$ (7)

Note that the real part of ${\mathcal{G}}_{\kappa}$ goes to infinity as $1/r$
when $r\rightarrow 0^{+}\,,$ while its imaginary part tends to
${\kappa}/{(4\pi)}$. The other kernel involved in (6) is
${\partial{{\mathcal{G}}}_{\kappa}}/{\partial{\bf n}_{\bf y}}$, which can be
written as follows, again separating its real and imaginary parts,

$\displaystyle\frac{\partial{{\mathcal{G}}}_{\kappa}}{\partial{\bf n}_{\bf
y}}({\bf x},{\bf y})=\frac{1}{4\pi}\ \left[\left(\frac{{\bf r}\cdot{\bf
n}_{\bf y}}{r^{3}}\left(\cos(\kappa r)+\kappa r\sin(\kappa
r)\right)\right)\,+\,{\rm i}\frac{{\bf r}\cdot{\bf n}_{\bf
y}}{r^{2}}\left(\frac{\sin(\kappa r)}{r}-\kappa\cos(\kappa r)\right)\right].$
(8)

For a regular surface $\Gamma$ described by a regular smooth parameterization,
the quantity ${{\bf r}\cdot{\bf n}_{\bf y}}/{r^{2}}$ is bounded for
$r\rightarrow 0^{+}$ but in general it is not continuous at $r=0$; see the
Appendix for the proof. The bounded behavior at $r=0$ of such quantity implies
that, analogously to ${{\mathcal{G}}}_{\kappa}$, the kernel
${\partial{{\mathcal{G}}}_{\kappa}}/{\partial{\bf n}_{\bf y}}$ is just weakly
singular, see Section 9.1 in [22] and also the analytical related insights
explicitly reported in the Appendix. Indeed, its real part goes to infinity as
$1/r$ when $r\rightarrow 0^{+}$ while its imaginary part tends to zero because
${{\bf r}\cdot{\bf n}_{\bf y}}/{r^{2}}$ remains bounded and it is multiplied
by a factor going to zero.

Clearly, when the boundary conditions in (3) are of Dirichlet type, the
unknown in (6) is the flux $\phi:=({\partial u}/{\partial{\bf n}})|\Gamma,$
belonging to $H^{-1/2}(\Gamma)$ and the CBIE becomes

$\displaystyle\int_{\Gamma}{{\mathcal{G}}}_{\kappa}({\bf x},{\bf y})\phi({\bf
y})d\Gamma_{{\bf y}}\,=\,c({\bf x})u_{\rm D}({\bf
x})+\int_{\Gamma}\frac{\partial{{\mathcal{G}}}_{\kappa}}{\partial{\bf n}_{\bf
y}}({\bf x},{\bf y})u_{\rm D}({\bf y})d\Gamma_{{\bf y}},\qquad{\bf
x}\in\Gamma.$ (9)

Note that this is an integral equation of the general type

$V_{\kappa}\phi({\bf x})=f({\bf x}),\qquad{\bf x}\in\Gamma,$ (10)

where $f$ denotes a known function and the single layer operator
$V_{\kappa}:H^{-1/2}(\Gamma)\rightarrow H^{1/2}(\Gamma)$ is an elliptic
isomorphism defined as,

$\displaystyle V_{\kappa}\phi({\bf
x}):=\int_{\Gamma}{\mathcal{G}}_{\kappa}({\bf x},{\bf y})\phi({\bf
y})\,d\Gamma_{{\bf y}}.$

Conversely, when Neumann conditions are dealt with, the unknown in (6) is the
acoustic potential $\phi:=u|\Gamma,$ with $\phi\in H^{1/2}(\Gamma)$ and the
CBIE becomes

$\displaystyle\int_{\Gamma}\frac{\partial{{\mathcal{G}}}_{\kappa}}{\partial{\bf
n}_{\bf y}}({\bf x},{\bf y})\phi({\bf y})\,d\Gamma_{{\bf y}}+c({\bf
x})\phi({\bf x})\,=\,+\int_{\Gamma}{{\mathcal{G}}}_{\kappa}({\bf x},{\bf
y})u_{\rm N}({\bf y})d\Gamma_{{\bf y}}\qquad{\bf x}\in\Gamma.$ (11)

This is an integral equation of the general type

$(c(x)I+K_{\kappa})\phi({\bf x})=f({\bf x}),\qquad{\bf x}\in\Gamma,$ (12)

associated with the operator $cI+K_{\kappa},$ where $I$ denotes the identity
operator and $K_{\kappa}:H^{1/2}(\Gamma)\rightarrow H^{1/2}(\Gamma)$ is the
following double layer operator,

$\displaystyle K_{\kappa}\phi({\bf
x}):=\int_{\Gamma}\frac{\partial{{\mathcal{G}}}_{\kappa}}{\partial{\bf n}_{\bf
y}}({\bf x},{\bf y})\ \phi({\bf y})\,d\Gamma_{{\bf y}}.$

Note that in our case the right hand side $f$ in both (10) and (12) is just
the right hand side respectively of (9) and (11) and also that clearly it is
in any case $\phi=\phi({\bf x},\kappa)$.

If $\phi$ is available, the solution $u=u({\bf x},\kappa)$, ${\bf x}\in\Omega$
of the BVP in (3) is obtained by using the so-called representation formula,

$\displaystyle u({\bf
x},\kappa)=\pm\left(\int_{\Gamma}{{\mathcal{G}}}_{\kappa}({\bf x},{\bf y})\
\frac{\partial u}{\partial{\bf n}}({\bf y})\ d\Gamma_{{\bf
y}}-\int_{\Gamma}\frac{\partial{{\mathcal{G}}}_{\kappa}}{\partial{\bf n}_{\bf
y}}({\bf x},{\bf y})\ u({\bf y})\ d\Gamma_{{\bf y}}\right),\qquad{\bf
x}\in\Omega,$ (13)

where the sign is positive if $\Omega=\Omega^{[i]}$ (interior problem) and
negative otherwise. We note that the numerical implementation of this formula
is not completely trivial because near singular integrals have to be
approximated, when ${\bf x}$ is very near to $\Gamma$. However, one can be
merely interested in $\phi$ or in evaluating $u$ just at points sufficiently
far from $\Gamma.$ For example, when an exterior problem is taken into
account, often the interest is just in the far field pattern of $u,$ that is
into the recovery of $u_{\infty}\,,$ where

$\displaystyle u_{\infty}({\bf w}):=\lim_{r\rightarrow\infty}re^{-i\kappa
r}u(r{\bf w})\,,\qquad{\bf w}\in S^{2},$

where $S^{2}$ is a unit sphere in ${\mathop{{\rm I}\kern-1.99997pt{\rm
R}}\nolimits}^{3}$ centered at the origin of the considered reference system.
In such case we can rely on the following formula whose computation just
requires rules for regular integrals [15],

$\displaystyle u_{\infty}({\bf w})=\frac{1}{4\pi}\,\int_{\Gamma}\left(i\kappa
u({\bf y})({\bf w}\cdot{\bf n}({\bf y}))\,+\,\frac{\partial u}{\partial{\bf
n}}({\bf y})\right)e^{-i\kappa({\bf w}\cdot{\bf y})}d\Gamma_{{\bf y}}.$

##  3 Multi-patch Isogeometric model

In this section we introduce the IgA-BEM model, with respect to the geometry
representation and the discretization space to express the numerical solution.
The flexibility to model complex geometries in the 3D case is facilitated by
adopting its multi-patch formulation. We assume that the boundary ${\Gamma}$
is a union of $M$ patches $\Gamma^{(\ell)},\ell=1,\ldots,M,$ and for $\ell\neq
k$ it holds $\Gamma^{(\ell)}\cap\Gamma^{(k)}=\emptyset$ and
$\partial\Gamma^{(\ell)}\cap\partial\Gamma^{(k)}$ is a boundary edge curve of
each patch, a corner point of each patch or an empty set. To each patch
$\Gamma^{(\ell)}$ we assign a geometry map ${\bf F}^{(\ell)}$, such that
$\overline{\Gamma}^{(\ell)}=\mbox{Image}({\bf F}^{(\ell)})$, ${\bf
F}^{(\ell)}~{}:~{}[0,1]^{2}\rightarrow\overline{\Gamma}^{(\ell)}\subset{\mathop{{\rm
I}\kern-1.99997pt{\rm R}}\nolimits}^{3}$ is regular and belongs to
$C^{2}([0,1]^{2})$. Each ${\bf F}^{(\ell)}$ is written in the following
standard NURBS form [2],

$\displaystyle{\bf F}^{(\ell)}({\bf t}):=\frac{\displaystyle{\sum_{{\bf
i}\in{\cal J_{\rm g}}^{(\ell)}}}w_{{\bf i}}^{(\ell)}\,{\bf Q}_{{\bf
i}}^{(\ell)}\,\hat{B}_{{\bf i},{\bf d}_{\rm g}}^{(\ell)}({\bf
t})}{\displaystyle{\sum_{{\bf i}\in{\cal J_{\rm g}}^{(\ell)}}}w_{{\bf
i}}^{(\ell)}\,\hat{B}_{{\bf i},{\bf d}_{\rm g}}^{(\ell)}({\bf t})},\quad{\bf
t}\in[0,1]^{2},$ (14)

where a set $\\{{\bf Q}_{{\bf i}}^{(\ell)}\in{\mathop{{\rm
I}\kern-1.99997pt{\rm R}}\nolimits}^{3},{\bf i}\in{\cal J_{\rm
g}}^{(\ell)}\\}$ defines a net of control points which, together with the
associated set of positive weights $\\{w_{{\bf i}}^{(\ell)}\in{\mathop{{\rm
I}\kern-1.99997pt{\rm R}}\nolimits},{\bf i}\in{\cal J_{\rm g}}^{(\ell)}\\},$
is the basic element typically used in the CAD environment to design free-form
surfaces. The tensor product B-spline basis function $\hat{B}_{{\bf i},{\bf
d}_{\rm g}}^{(\ell)}({\bf t}):=\hat{B}_{i_{1},d_{\rm g,1}}^{(\ell)}(t_{1})\
\hat{B}_{i_{2},d_{\rm g,2}}^{(\ell)}(t_{2})$, where ${\bf
i}:=(i_{1},i_{2})\in\cal J_{\rm g}^{(\ell)}$, is of bi-degree ${\bf d}_{\rm
g}=(d_{\rm g,1},d_{\rm g,2})$, defined on $[0,1]^{2}$ in variable ${\bf
t}:=(t_{1},t_{2})$ and with respect to the bidimensional knot array $T_{\rm
g}^{(\ell)}$ associated with the $\ell$-th patch, $T_{\rm g}^{(\ell)}:=T_{\rm
g,1}^{(\ell)}\times T_{\rm g,2}^{(\ell)}$, and $T_{\rm g,1}^{(\ell)}$ and
$T_{{\rm g},2}^{(\ell)}$ are both clamped knot vectors with entries in the
interval $[0,1]$. When two patches share an edge ${\cal
C}_{\ell,k}:=\partial{\Gamma}^{(\ell)}\cap\partial{\Gamma}^{(k)}$ for
$\ell\neq k$, the shared edge is a NURBS curve. For simplicity we usually
assume it to be parameterized with the same geometry map for both patches, up
to reversion in directions and direction swapping in the parametric space.

As commonly done also in the literature, in order to gain more flexibility, in
our generalized IgA-BEM setting the unknown Cauchy datum $\phi$ is
approximated in a functional space whose restriction to
$\bar{\Gamma}^{(\ell)}$ is the geometric counterpart of the space generated
from bivariate B-splines always defined in $[0,1]^{2}$ but not necessarily
obtained with a refinement procedure from $\hat{B}_{{\bf i},{\bf d}_{\rm
g}}^{(\ell)},\,{\bf i}\in{\cal J_{\rm g}}^{(\ell)}$. Such space is
characterized by its own degree ${\bf d}=(d_{1},d_{2})$ (common to all
patches) and its extended knot vectors, ${T}_{j}^{(\ell)},j=1,2.$ Then,
denoting with $h_{j}^{(\ell)}$ the maximal distance between successive knots
in ${T}_{j}^{(\ell)},$ setting
$h^{(\ell)}:=\max\\{h_{1}^{(\ell)},\,h_{2}^{(\ell)}\\}$ and introducing a set
of bivariate indices ${\cal J}^{(\ell)}$ for the basis functions, we can
introduce $\hat{\cal S}^{(\ell)}_{h^{(\ell)}}$ which denotes the spline space
generated by the tensor product B-spline basis $\\{\hat{B}_{{\bf j},{\bf
d}}^{(\ell)}:{\bf j}=(j_{1},j_{2})\in{\cal J}^{(\ell)}\\}$ of bi-degree ${\bf
d}$ defined on $[0,1]^{2}$ with respect to the bidimensional knot array
${T}^{(\ell)}:={T}_{1}^{(\ell)}\times{T}_{2}^{(\ell)}$.

###  3.1 Discontinuous basis

When no continuity constraint is imposed across patches for the basis
functions, the global space used to approximate the unknown Cauchy datum
$\phi$ is

$\displaystyle{\cal S}_{h,{\bf d}}\,:=\,\mathrm{span}\\{{B}_{{\bf j},{\bf
d}}^{(\ell)}:{\bf j}\in{\cal J}^{(\ell)},\quad\ell=1,\ldots,M\\},$

where $h:=\max\\{h^{(\ell)},\ell=1,\ldots,M\\}$ and

$\displaystyle{B}_{{\bf j},{\bf d}}^{(\ell)}({\bf
x})=\left\\{\begin{array}[]{cc}\hat{B}_{{\bf j},{\bf d}}^{(\ell)}\circ{{\bf
F}^{(\ell)}}^{-1}({\bf x}),&\mbox{if }{\bf x}\in{\Gamma}^{(\ell)},\cr
0,&\mbox{otherwise.}\end{array}\right.$

Clearly with this setting the global number $N_{\rm DOF}$ of degrees of
freedom used to approximate $\phi$ is

$\displaystyle{N}_{\rm DOF}=\sum_{\ell=1}^{M}|{\cal J}^{(\ell)}|,$

where $|{\cal J}^{(\ell)}|$ is the cardinality of the set. In this setting it
is convenient to select $|\cal J^{(\ell)}|$ distinct collocation points
belonging to the interior of each patch $\Gamma^{(\ell)}$ , $\ell=1,\ldots,M.$
In particular we fix ${\bf x}_{\bf j}^{(\ell)}={\bf F}^{(\ell)}({\bf s}_{\bf
j})$, ${\bf j}\in{\cal J}^{(\ell)},$ where $\\{{\bf s}_{\bf j}^{(\ell)}:{\bf
j}\in{\cal J}_{\ell}\\}$ denotes a set of points in $(0,1)^{2}$, defined as
the Cartesian product of two sets of abscissas, the so-called improved
Greville abscissas associated with ${T}_{1}^{(\ell)}$ and ${T}_{2}^{(\ell)}$
[23]. We recall that the improved Greville abscissae related to a set of
Greville abscissae $\bar{\xi}_{1},\bar{\xi}_{2},\ldots,\bar{\xi}_{n}$, are
Greville abscissae whose first and last elements are modified as follows:

$\bar{\xi}_{1}=\bar{\xi}_{1}+\omega(\bar{\xi}_{2}-\bar{\xi}_{1}),\;\;\;\bar{\xi}_{n}=\bar{\xi}_{n}-\omega(\bar{\xi}_{n}-\bar{\xi}_{n-1})\,,\quad\mbox{with}\,\,0<\omega<1\,.$
(15)

Thus, a discrete version of the integral formulation of the Dirichlet (9) or
Neumann (11) problem is obtained by approximating $\phi$ in the finite
dimensional composite space ${\cal S}_{h,{\bf d}}$. The applied collocation
method leads to a linear system

$\displaystyle A\,\mbox{\boldmath$\alpha$}=\mbox{\boldmath$\beta$}.$ (16)

The unknown entries in the vector $\alpha$ are the coefficients which allow us
to define our approximation $\phi_{h}$ of $\phi$ in the space ${\cal
S}_{h,{\bf d}}$ patch-wisely as follows,

$\displaystyle\phi_{h}({\bf x}):=\sum_{\ell=1}^{M}\sum_{{\bf j}\in{\cal
J}^{(\ell)}}\alpha^{(\ell)}_{\bf j}{B}_{{\bf j},{\bf d}}^{(\ell)}({\bf
x}),\qquad{\bf x}\in{\Gamma}.$

The square system matrix $A\in\mathbb{C}^{{N_{\rm DOF}}\times N_{\rm DOF}}$
and the right-hand side vector $\mbox{\boldmath$\beta$}\in\mathbb{C}^{N_{\rm
DOF}}$ have a block representation,
$A=(A^{(\ell,k)}),\,\,\mbox{\boldmath$\beta$}=(\mbox{\boldmath$\beta$}^{(\ell)})$,
for $\ell,k=1,\dots,M$, where $A^{(\ell,k)}\in\mathbb{C}^{|{\cal
J}_{\ell}|\times|{\cal J}_{k}|}$ and
$\mbox{\boldmath$\beta$}^{(\ell)}\in\mathbb{C}^{|{\cal J}_{\ell}|}$. For each
pair of patches $\Gamma^{\ell},\Gamma^{k}$, the rows (columns) of the related
submatrix are ordered by using a lexicographical ordering of the elements in
$\mathcal{J}^{(\ell)}$ ($\mathcal{J}^{(k)}$), namely we identify a single
index $i$ with ${\bf i}$ and $j$ with ${\bf j}$. After this simplification,
for the Dirichlet case such entries read as follows in both the physical and
the parametric domain,

$\displaystyle A^{(\ell,k)}_{i,j}:=$
$\displaystyle\int_{\Gamma^{(k)}}{{\mathcal{G}}}_{\kappa}({\bf x}_{\bf
i}^{(\ell)},{\bf y})B^{(k)}_{{\bf j},{\bf d}}({\bf y})\
d\Gamma\,=\int_{[0,1]^{2}}{{\mathcal{G}}}_{\kappa}({\bf x}_{\bf
i}^{(\ell)},{\bf F}^{(k)}({\bf t}))\hat{B}^{(k)}_{{\bf j},{\bf d}}({\bf t})\
J^{(k)}({\bf t})\ d{\bf t},$ (17)
$\displaystyle\mbox{\boldmath$\beta$}^{(\ell)}_{i}:=$
$\displaystyle\sum_{k=1}^{M}\int_{\Gamma^{(k)}}\frac{\partial{\mathcal{G}}_{\kappa}}{\partial{\bf
n}_{\bf y}}({\bf x}_{\bf i}^{(\ell)},{\bf y})\ u_{\rm D}({\bf y})\
d\Gamma+\frac{1}{2}u_{\rm D}({\bf x}_{\bf i}^{(\ell)})=$
$\displaystyle\sum_{k=1}^{M}\
\int_{[0,1]^{2}}\frac{\partial{\mathcal{G}}_{\kappa}}{\partial{\bf n}_{\bf
y}}({\bf x}_{\bf i}^{(\ell)},{\bf F}^{(k)}({\bf t}))\ u_{\rm D}({\bf
F}^{(k)}({\bf t}))\ J^{(k)}({\bf t})\ d{\bf t}+\frac{1}{2}u_{\rm D}({\bf
x}_{\bf i}^{(\ell)}),$ (18)

where $i=1,\ldots,|{\cal J}^{(\ell)}|,\quad j=1,\ldots,|{\cal
J}^{(k)}|,\,\,\hat{B}^{(k)}_{{\bf j},{\bf d}}\in{\hat{\cal
S}}^{(k)}_{h^{(k)}}$ and $J^{(k)}$ represents the infinitesimal surface area
element on the $k$-th patch,

$J^{(k)}(\cdot):=\left\|\frac{\partial{\bf F}^{(k)}}{\partial
t_{1}}(\cdot)\times\frac{\partial{\bf F}^{(k)}}{\partial
t_{2}}(\cdot)\right\|_{2}.$

Conversely, for the Neumann case we have

$\displaystyle A^{(\ell,k)}_{i,j}=$
$\displaystyle:=\int_{\Gamma^{(k)}}\frac{\partial{\mathcal{G}}_{\kappa}}{\partial{\bf
n}_{\bf y}}({\bf x}_{\bf i}^{(\ell)},{\bf y})B_{{\bf j},{\bf d}}^{(k)}({\bf
y})\ d\Gamma+\frac{1}{2}B_{{\bf j},{\bf d}}^{(k)}({\bf x}_{\bf i}^{(\ell)})=$
$\displaystyle\int_{[0,1]^{2}}\frac{\partial{\mathcal{G}}_{\kappa}}{\partial{\bf
n}_{\bf y}}({\bf x}_{\bf i}^{(\ell)},{\bf F}^{(k)}({\bf t}))\hat{B}_{{\bf
j},{\bf d}}^{(k)}({\bf t})J^{(k)}({\bf t})\ d{\bf t}+\frac{1}{2}\hat{B}_{{\bf
j},{\bf d}}^{(k)}({\bf x}_{\bf i}^{(\ell)}),$ (19)
$\displaystyle\mbox{\boldmath$\beta$}^{(\ell)}_{i}:=$
$\displaystyle\sum_{k=1}^{M}\ \int_{\Gamma^{(k)}}{\mathcal{G}}_{\kappa}({\bf
x}_{\bf i}^{(\ell)},{\bf y})\ u_{\rm N}({\bf y})\ d\Gamma\,=\sum_{k=1}^{M}\
\int_{[0,1]^{2}}{\mathcal{G}}_{\kappa}({\bf x}_{\bf i}^{(\ell)},{\bf
F}^{(k)}({\bf t}))\ u_{\rm N}({\bf F}^{(k)}({\bf t}))\ J^{(k)}({\bf t})\ d{\bf
t},$ (20)

Referring for brevity just to the Dirichlet case, we report below an
equivalent scaled expression of these entries, since it will be useful in
Section 4 to develop a reasonably scaled procedure to detect near singular
integrals to be computed during the assembly phase. Denoting with $\mu_{k}$ a
reference length for the $k$-th patch (e.g., the square root of its area or
its diameter), with $\varphi_{k}$ the associated uniform scaling such that
$\tilde{\bf z}=\varphi_{k}({\bf z})={\bf z}/\mu_{k}$ and with
$\tilde{\Gamma}^{(k)}$ the conformally scaled version of the $k$–th patch, we
have also

$\displaystyle A^{(\ell,k)}_{i,j}=$ $\displaystyle\mu_{k}\
\int_{\tilde{\Gamma}^{(k)}}{\mathcal{G}}_{\kappa_{k}}(\varphi_{k}({\bf x}_{\bf
i}^{(\ell)}),\tilde{{\bf y}})B^{(k)}_{{\bf j},{\bf d}}(\mu_{k}\tilde{{\bf
y}})\ d\tilde{\Gamma}=$ $\displaystyle\mu_{k}\
\int_{[0,1]^{2}}{{\mathcal{G}}}_{\kappa_{k}}(\varphi_{k}({\bf x}_{\bf
i}^{(\ell)}),\tilde{{\bf F}}^{(k)}({\bf t}))\hat{B}^{(k)}_{{\bf j},{\bf
d}}({\bf t})\ \tilde{J}^{(k)}({\bf t})\ d{\bf t},$ (21)
$\displaystyle\mbox{\boldmath$\beta$}^{(\ell)}_{i}=$
$\displaystyle\sum_{k=1}^{M}\
\int_{\tilde{\Gamma}^{(k)}}\frac{\partial{{\mathcal{G}}}_{\kappa^{(k)}}}{\partial{\bf
n}_{\tilde{y}}}(\varphi_{k}({\bf x}_{\bf i}^{(\ell)}),\tilde{\bf y})\ u_{\rm
D}(\mu_{k}\tilde{\bf y})\ d\tilde{\Gamma}+\frac{1}{2}u_{\rm D}({\bf x}_{\bf
i}^{(\ell)})=$ $\displaystyle\sum_{k=1}^{M}\
\int_{[0,1]^{2}}\frac{\partial{\mathcal{G}}_{\kappa_{k}}}{\partial{\bf
n}_{\tilde{{\bf y}}}}(\varphi_{k}({\bf x}_{\bf i}^{(\ell)}),\tilde{{\bf
F}}^{(k)}({\bf t}))\ u_{\rm D}(\mu_{k}\tilde{{\bf F}}^{(k)}({\bf t}))\
\tilde{J}^{(k)}({\bf t})\ d{\bf t}+\frac{1}{2}u_{\rm D}({\bf x}_{\bf
i}^{(\ell)}),$ (22)

where $\tilde{\bf F}^{(k)}:={\bf F}^{(k)}/\mu_{k},\,\tilde{J}^{(k)}$ is the
related infinitesimal scaled area element and $\kappa_{k}:=\kappa\mu_{k}.$

###  3.2 $C^{0}$ continuous basis

Let us denote the space of B-spline functions that join with $C^{0}$
regularity across patches by ${\cal S}^{0}_{h,{\bf d}}:={\cal S}_{h,{\bf
d}}\cap C^{0}({\Gamma})$. For simplicity let us assume that any two adjacent
patches $\Gamma^{(\ell)},\Gamma^{(k)}$ with a common boundary curve ${\cal
C}_{\ell,k}$ have the same knot vector on the edges of the patches. Namely, on
$\ell$-th patch we consider the knot vector $T_{1}^{(\ell)}$ if ${\bf
F}^{(\ell)}(t_{1},0)$ or ${\bf F}^{(\ell)}(t_{1},1)$ define the curve ${\cal
C}_{\ell,k}$, and $T_{2}^{(\ell)}$ in the other case. The same procedure is
done for $k$-th patch. Our assumption then states that the two vectors
$T_{i}^{(\ell)}$, $T_{j}^{(k)}$ need to coincide (up to a possible reversion
of one knot sequence, e.g., when ${\bf F}^{(\ell)}(1,t_{2})={\bf
F}^{(k)}(1,1-t_{2})$). For a more general setting to describe the geometry and
the discretization space without the matching of adjacent knot vectors we
refer to [1].

The so-called interior basis functions for ${\cal S}^{0}_{h,{\bf d}}$ includes
all basis functions ${B}_{{\bf j},{\bf d}}^{(\ell)}$ from ${\cal S}_{h,{\bf
d}}$, whose corresponding B-splines $\hat{B}_{{\bf j},{\bf d}}^{(\ell)}$
vanish on $\partial([0,1]^{2})$. The so-called edge or vertex basis functions
are obtained by identifying and summing together the remaining basis functions
from ${\cal S}_{h,{\bf d}}$ with shared knot vectors. More precisely, an edge
basis function $B\in{\cal S}^{0}_{h,{\bf d}}$ associated to ${\cal
C}_{\ell,k}$ is a sum of two remaining (not interior) functions $B_{{\bf
i},{\bf d}}^{(\ell)}$ and $B_{{\bf j},{\bf d}}^{(k)}$ from ${\cal S}_{h,{\bf
d}}$, that vanish at the endpoints of ${\cal C}_{\ell,k}$ and satisfy
$\hat{B}_{{\bf i},{\bf d}}^{(\ell)}\circ{{\bf F}^{(\ell)}}^{-1}({\bf
x})=\hat{B}_{{\bf j},{\bf d}}^{(k)}\circ{{\bf F}^{(k)}}^{-1}({\bf x})$ for
every ${\bf x}\in{\cal C}_{\ell,k}$. Similarly, a vertex function associated
to a patch corner vertex ${\bf x}_{v}$, is defined as a sum of all basis
functions $B_{{\bf i},{\bf d}}^{(\ell)}$ from ${\cal S}_{h,{\bf d}}$,
$\ell=1,\dots,M$ and ${\bf i}\in\cal J^{(\ell)}$ , that satisfy $\hat{B}_{{\bf
i},{\bf d}}^{(\ell)}\circ{{\bf F}^{(\ell)}}^{-1}({\bf x}_{v})\neq 0$.

In all our experiments with the $C^{0}$ continuous basis we have fixed the
collocation points as a set of images of Cartesian products of standard
Greville abscissas on $[0,1]^{2}$, mapped to $\overline{\Gamma}^{(\ell)}$ via
map ${\bf F}^{(\ell)}$ for all $\ell=1,\dots,M$, and by counting just once the
repeated points, which appear on the boundary curves of patches.

The approximate solution $\phi_{h}$ is constructed similarly to the
discontinuous case, bearing in mind a variation in the collocation points and
a sum of adequate entries in (17), (3.1) over all suitable $k$ to evaluate
integrals that span over the patches $\Gamma^{(k)}$ adjacent to
$\Gamma^{(\ell)}$.

##  4 Quadratures

###  4.1 Singularity extraction

Let us introduce our novel numerical integration method which makes possible a
function–by–function assembly phase. For the sake of simplicity, we refer just
to the one patch setting, i.e., we assume $M=1,$ thus avoiding any reference
to the patch index, when the source point and the integration domain belong to
the same patch. Some considerations about integration in the more general
multi-patch setting are also given.

It is convenient to write the two singular kernels ${\mathcal{G}}_{\kappa}$
and ${\partial{{\mathcal{G}}}_{\kappa}}/{\partial{\bf n}_{\bf y}}$ as in (7)
and (8) and analyze which terms are singular, hence a singular integration has
to be considered, and which terms define regular integrals. Referring to the
Appendix for some explicit insights, as already noted in Section 2, the real
part in (7) is weakly singular at $r=0$ with its non-singular part
$\cos(\kappa r)$ smooth under the assumption of a smooth geometry
parameterization. For the same reason, its imaginary part is a smooth
function. Similarly, the factor ${\bf r}\cdot{\bf n}_{\bf y}/r^{3}$ of the
real part in (8) is weakly singular, while the other factor, $(\cos(\kappa
r)+\kappa\sin(\kappa r)r),$ is smooth. Finally, the imaginary part in (8) is
in general only $C^{1}$ continuous at $r=0$ because the factor ${{\bf
r}\cdot{\bf n}_{\bf y}}/{r^{2}}$ is bounded and $(\sin(\kappa
r)/r-\kappa\cos(\kappa r))$ has a double zero at $r=0$. In our numerical
experiments such smoothness of the integrands was sufficient to numerically
integrate the corresponding integrals with the required accuracy.

###### Remark 1

To improve accuracy of the numerical integration for integrals involving the
imaginary part of (8), a more sophisticated approach could be considered. For
example, the integrals could be handled similarly as the singular ones,
treating the part ${{\bf r}\cdot{\bf n}_{\bf y}}/{r^{2}}$ as the “less regular
part” of the integrand, since the remaining part $\left({\sin(\kappa
r)}/r-\kappa\cos(\kappa r)\right)$ is an analytic function.

In order to give a compact introduction of our quadrature approach, first we
denote with ${\bf s}\in[0,1]^{2}$ the pre-image in the parametric domain of a
fixed generic collocation point ${\bf x}$. Note that for brevity in the
following we refer to the patch index only when it is essential. Conversely,
we refer to the patch scaled integrals appearing for example in formulas
(3.1)–(3.1), as well as to the related notation also introduced in Section 3.
So for example we denote with $\tilde{\bf x}$ the scaling of ${\bf x}$,
assuming that the considered scaling is that associated with the patch where
the integration is developed. Then, we can introduce the following general
form representing any integral to be dealt with,

$\displaystyle\int_{R_{\bf j}}U({\bf s},{\bf t})\hat{B}_{{\bf j},{\bf d}}({\bf
t})\ g({\bf t})\,d{\bf t},$ (23)

where $g$ is a given smooth function, $\hat{B}_{{\bf j},{\bf d}}$ is a
B-spline basis function, $R_{\bf j}:=\mathop{\mathrm{supp}}\\{\hat{B}_{{\bf
j},{\bf d}}\\}$ and $U$ refers to any of the following two weakly singular
kernels respectively associated with the single and double layer potential,

$\displaystyle U({\bf s},{\bf t})=U_{\rm SL}({\bf s},{\bf
t}):=\frac{1}{\|\tilde{\bf x}({\bf s})-\tilde{\bf F}({\bf
t})\|_{2}}\qquad\textrm{and}\qquad U({\bf s},{\bf t})=U_{\rm DL}({\bf s},{\bf
t}):={\frac{(\tilde{\bf x}({\bf s})-\tilde{\bf F}({\bf
t}))\cdot\tilde{}\mbox{\boldmath$\nu$}({\bf t})}{\|\tilde{\bf x}({\bf
s})-\tilde{\bf F}({\bf t})\|_{2}^{3}}},$ (24)

with $\tilde{}\mbox{\boldmath$\nu$}({\bf t}):=({\partial\tilde{\bf
F}}/{\partial t_{1}})({\bf t})\times({\partial\tilde{\bf F}}/{\partial
t_{2}})({\bf t})$ and $\tilde{\bf F}$ denoting the scaled normal and the
scaled mapping associated to the conformally scaled surface patch.

###### Remark 2

Integrals defining the entries of $\beta$ do not include the basis functions
$\hat{B}_{{\bf j},{\bf d}}$ directly. To split the integration domain
$[0,1]^{2}$ into smaller subdomains and arrive to the above general form, we
introduce additional knot vectors, $T_{{\rm r},1},T_{{\rm r},2}$, and inside
(3.1) and (20) we add a sum of all B-spline basis functions $\hat{B}_{{\bf
j},\boldsymbol{0}}$ of degree zero, defined on $T_{{\rm r},1}\times T_{{\rm
r},2}$. Therefore, the approach to evaluate integrals in $\beta$ is analogous
to the one for the system matrix $A$.

It is easy to determine when the integral in (23) is actually singular.
Indeed, since we are dealing with smooth surfaces without self-intersections,
this occurs if and only if $\tilde{\bf x}$ belongs to the scaled patch where
the integration is developed, and the associated ${\bf s}$ belongs to $R_{\bf
j}$. On the other hand, detection of nearly singular integrals that require
analogous treatment is more delicate if we want to preserve the accuracy of
the numerical scheme without deteriorating its efficiency. In such setting we
propose to consider nearly singular the integral in (23) when

$\displaystyle\tilde{r}_{{\rm min},{\bf j}}({\bf x}):=\min_{{\bf t}\in R_{\bf
j}}\|\tilde{\bf x}-\tilde{\bf F}({\bf t})\|_{2}\leq\,c\ \delta,$ (25)

where $0<c<1$ is a user-defined constant and $\delta$ is a suitable patch-
depended threshold. For further details concerning the definition of $\delta$
see the discussion in Subection 4.5.

When the condition (25) is not satisfied, the kernel $U$ is considered regular
and integration rules for regular integrals are applied. Conversely, when the
kernel $U$ in (23) is weakly singular or near singular, we preliminarily adopt
a regularization technique based on the subtraction of singularity,

$\displaystyle\int_{R_{\bf j}}U({\bf s},{\bf t})\hat{B}_{{\bf j},{\bf d}}({\bf
t})g({\bf t})\,d{\bf t}=\int_{R_{\bf j}}\big{(}U({\bf s},{\bf t})-{U}_{{\bf
s}}^{m}({\bf s}-{\bf t})\big{)}\hat{B}_{{\bf j},{\bf d}}({\bf t})g({\bf
t})\,d{\bf t}+\int_{R_{\bf j}}{U}_{{\bf s}}^{m}({\bf s}-{\bf t})\hat{B}_{{\bf
j},{\bf d}}({\bf t})g({\bf t})\,d{\bf t},$ (26)

where ${U}_{{\bf s}}^{m}$ is an approximation of $U({\bf s},\bullet)$,
obtained by truncating a particular series expansion of $U$ about the singular
point ${\bf t}={\bf s}$ after the $m$-th term; see [17] for the first
application of these expansions in the simplest form and [19] for detailed
construction and analysis. Note that the singularity subtraction technique
that splits an integral into two ones is applied only to a small portion of
integrals. In this case, the first term of the right hand side of (26) becomes
a _regularized integral_ while the second one is a singular integral with a
simplified singular kernel. The continuity of the integrand in the regularized
integral is controlled by the number $m$ of terms in $U_{{\bf s}}^{m}$ and is
equal to $C^{m-2}$ at ${\bf t}={\bf s}$ for the chosen $m\geq 1$ (with
$C^{-1}$ we denote functions that are integrable in the standard sense but
have discontinuities of the first kind at ${\bf t}={\bf s}$), if the geometry
map ${\bf F}$ is a smooth enough function [19].

By defining local coordinates ${\bf z}:={\bf s}-{\bf t}$, the kernel reads

$\displaystyle U^{m}_{\bf s}({\bf z})=\sum_{k=1}^{m}R_{\bf s}({\bf
z})^{-2(k+\gamma)+1}\;P_{{\bf s},3k+2\gamma-3}^{[k]}({\bf z}),$ (27)

where $R_{\bf s}({\bf z}):=\sqrt{{\bf z}^{\top}M_{\bf s}{\bf z}}$ is a square
root of a bivariate quadratic homogeneous polynomial and $M_{\bf s}$ is the
$2\times 2$ symmetric positive definite matrix associated with the first
fundamental form of a smooth surface (e.g., see formulas (45) and (47) in the
Appendix), $P_{{\bf s},3k+2\gamma-3}^{[k]}$ are suitable homogeneous
polynomials of degrees $3k+2\gamma-3$, and $\gamma=0$ and $\gamma=1$ for
${\mathcal{G}}_{\kappa}$ and ${\partial{{\mathcal{G}}}_{\kappa}}/{\partial{\bf
n}_{\bf y}}$, respectively. Note that coefficients inside $R_{\bf s}$ and
$P_{{\bf s},3k+2\gamma-3}^{[k]}$ depend on a chosen ${\bf s}$ and on
derivatives of the map $\tilde{\bf F}$. For example, for $m=1$ the kernels of
the single and double layer potentials are approximated as follows,

$\displaystyle U^{1}_{\bf s}({\bf
z})=\left\\{\begin{array}[]{ll}\frac{1}{4\pi}\left[R_{\bf s}({\bf
z})^{-1}\right]&\mbox{if }U=U_{\rm SL}\,,\cr\frac{1}{4\pi}\left[R_{\bf s}({\bf
z})^{-3}P^{[1]}_{{\bf s},2}({\bf z})\right]&\mbox{if }U=U_{\rm
DL}\,,\cr\end{array}\right.$

where $P^{[1]}_{{\bf s},2}({\bf z})=\frac{1}{2}\left(L\ z_{1}^{2}+2M\
z_{1}z_{2}+Nz_{2}^{2}\right)$, and $L,M,N$ are the coefficients of the second
fundamental form of the geometry, evaluated at ${\bf s}$, see (46) in the
Appendix.

Relating to the more general multi-patch setting, when source point ${\bf
x}={\bf F}^{(\ell)}({\bf s})$ and integration variable ${\bf y}={\bf
F}^{(k)}({\bf t})$ belong to different adjacent patches, $\Gamma^{(\ell)}$ and
$\Gamma^{(k)}$, the integral in (23) is nearly singular if condition (25) is
satisfied, and regular otherwise. In the latter case, the procedure is that
adopted for regular integrals on one patch while the former case clearly poses
an additional challenge. In order to explain our procedure in such setting, we
need first to make the notation more precise for $\ell\neq k$, adding the
superscript $(\ell,k)$ when referring to the kernel $U,\ U\rightarrow
U^{(\ell,k)},$ whenever $U_{\rm SL}$ or $U_{\rm DL}$ is referred to,

$\displaystyle U^{(\ell,k)}_{\rm SL}({\bf s},{\bf t}):=\frac{1}{\|\tilde{\bf
F}^{(\ell)}({\bf s})-\tilde{\bf F}^{(k)}({\bf t})\|_{2}},\quad
U^{(\ell,k)}_{\rm DL}({\bf s},{\bf t}):={\frac{(\tilde{\bf F}^{(\ell)}({\bf
s})-\tilde{\bf F}^{(k)}({\bf t}))\cdot\tilde{}\mbox{\boldmath$\nu$}^{(k)}({\bf
t})}{\|\tilde{\bf F}^{(\ell)}({\bf s})-\tilde{\bf F}^{(k)}({\bf
t})\|_{2}^{3}}},$ (28)

indicating that the source point $\tilde{\bf x}$ is connected with ${\bf s}$
in the parametric domain via the geometry map $\tilde{\bf F}^{(\ell)}$ but
integration has to be developed in the $k$-th patch. Since the
parameterization of the geometry is prescribed patch-wisely, the proposed
series expansion $U_{\bf s}^{(\ell,k),m}$ is limited to the $\ell$-th patch.
On the other hand, since ${\bf F}^{(k)}$ is regular, we can (slightly) enlarge
its definition outside the original domain $[0,1]^{2}$ and still have a
regular geometry map; for convenience let us use the same notation for the
extended function, ${\bf F}^{(k)}:[0,1]^{2}\subset
D_{k}\rightarrow{\mathop{{\rm I}\kern-1.99997pt{\rm R}}\nolimits}^{3}$.
Clearly on the enlarged part of the domain the map ${\bf F}^{(k)}$ describes
just an approximation of the actual geometry, which is exactly described by
${\bf F}^{(\ell)}$. By replacing $\tilde{\bf x}=\tilde{\bf F}^{(\ell)}({\bf
s})$ in (28) by $\tilde{\bf F}^{(k)}({\bf s}_{e})$ with

$\displaystyle{\bf s}_{e}:=\arg\min_{{\bf s}^{\prime}\in D_{k}}\|\tilde{\bf
x}-\tilde{\bf F}^{(k)}({\bf s}^{\prime})\|_{2}\,,$

we obtain the kernel $U_{{\bf s}_{e}}^{(k,k),m}$ that approximates $U_{{\bf
s}}^{(\ell,k),m}$ and naturally generalizes the case $\ell=k$. Since the
actual source point $\tilde{\bf x}$ is replaced by its nearby point
$\tilde{\bf F}^{(k)}({\bf s}_{e})$, the cancellation of the singularity is not
exact. This does not present any real problem because only the evaluation of
inter-patch nearly singular integrals is required.

As an example, let $\Gamma^{(\ell)}$ and $\Gamma^{(k)}$ be two adjacent NURBS
patches of a unit sphere (see the beginning of Section 5 for more details and
Figure 2(a) for visualization of the sphere representation). In Figure 1 the
function ${\bf t}\mapsto|U^{(\ell,k)}({\bf s},{\bf t})-U_{{\bf
s}_{e}}^{(k,k),2}({\bf s}_{e}-{\bf t})|$ is depicted for a chosen position of
${\bf x}={\bf F}^{(\ell)}({\bf s})\in\Gamma^{(\ell)}$ near the interface
between $\Gamma^{(\ell)}$ and $\Gamma^{(k)}$ (for simplicity, here we do not
consider any scaling, that is we assume $\mu_{k}=1$). Such function is the
absolute value of the regularized part of the integrand in (26) for the
kernels of the single and double layer potential. The figure shows its value
in the whole parametric domain $[0,1]^{2}$ of $\Gamma^{(k)}\,,$ together with
the position of ${\bf s}_{e}$ in $D_{k}.$ The distance between the computed
nearby point ${\bf F}^{(k)}({\bf s}_{e})$ and the actual source point ${\bf
x}={\bf F}^{(\ell)}({\bf s})$ is approximately $4.1\cdot 10^{-9}$ and
$9.5\cdot 10^{-9}$, when the source point is close to a vertex on the
interface (left figures) and to its central part (right figures),
respectively. The figure shows that the regularized function
$U^{(\ell,k)}({\bf s},{\bf t})-U_{{\bf s}_{e}}^{(k,k),2}({\bf s}_{e}-{\bf t})$
has small oscillations and is bounded–in comparison consider that the maximum
value of the function ${\bf t}\mapsto|U^{(\ell,k)}({\bf s},{\bf t})|$ in
$[0,1]^{2}$ is $41.9$ and $35.2$ for the two cases with the single layer
potential (top figures) and $25.0$, $31.6$ for the cases with the double layer
potential (bottom figures). Thus, it can be effectively integrated with the
quadrature rule for regular integrals presented in Subsection 4.3.

Figure 1:  An example of the absolute value of the regularized function of the
nearly singular kernel with the source point and integration domain belonging
to two different adjacent patches of a sphere. The top and bottom figures
correspond to regularization for the single and double layer potential,
respectively (with $m=2$). Two different configurations of points ${\bf
s}_{e}$, represented as a black dots, are considered on the left and on right.

###  4.2 B-spline quasi-interpolation scheme

To approximate integrals in 3D collocation BEM, we use the 2D Quasi
Interpolation (QI) scheme firstly introduced in [24] which is the tensor-
product extension of the univariate Hermite QI scheme introduced in [25]. With
such approach a bivariate function $f$ is approximated on a rectangle $R$ by
using a QI spline tensor product operator $f\rightarrow\sigma_{f}:=\sum_{{\bf
i}\in{\cal J}_{R}}\lambda_{\bf i}\hat{B}_{{\bf i},{\bf p}}$, where
$\\{\hat{B}_{{\bf i},{\bf p}},{\bf i}\in{\cal J}_{R}\\}$ is the tensor-product
B-spline basis generating the chosen finite dimensional space of splines
defined on $R,$ with ${\cal J}_{R}$ denoting a local set of multi–indices used
to identify the basis elements. In the considered QI scheme each coefficient
$\lambda_{\bf i},{\bf i}\in{\cal J}_{R},$ is computed as a linear combination
of the values assumed in a suitable local subset of breakpoints by $f,$ by
both its first partial derivatives and also by its second mixed derivative,
see [24] for details. However, in the present work we rely on a derivative-
free variant of such QI scheme, since we are interested to employ it for
numerical integration. Such variant is obtained approximating the required
derivative values by suitable finite difference formulas. Then, collecting
with lexicographical order all the coefficients $\lambda_{\bf i},{\bf
i}\in{\cal J}_{R}$ in a vector $\Lambda$, in compact matrix notation we have

$\displaystyle\mbox{\boldmath$\Lambda$}=\left({\hat{C}}^{(p)}\right)^{\top}{\bf
f}\,,$ (29)

where ${\bf f}$ is a vector collecting with the same ordering all the values
of $f$ at the breakpoints and

$\displaystyle{\hat{C}}^{(p)}:=\left(({\hat{A}_{2}}^{\top}\otimes{\hat{A}_{1}}^{\top})-({\hat{A}_{2}}^{\top}\otimes({\hat{B}_{1}}^{\top}D_{t_{1}}))-(({\hat{B}_{2}}^{\top}D_{t_{2}})\otimes{\hat{A}_{1}}^{\top})+(({\hat{B}_{2}}^{\top}D_{t_{2}})\otimes({\hat{B}_{1}}^{\top}D_{t_{1}}))\right)^{\top}\,.$
(30)

In the above formula $\hat{A_{i}},\,\hat{B_{i}}$ denote the two banded
matrices, explicitly reported for example in [26], which are associated to the
univariate original Hermite scheme formulated in the $i$-th direction,
$i=1,2;$ $D_{t_{i}}$ are the matrices associated with the two analogous
directional finite difference formulas used by the considered QI variant of
the original scheme for approximating the first partial derivatives. Clearly,
by incorporating a sufficiently accurate finite differences formula to
approximate the required derivatives, the resulting derivative-free QI scheme
shares the approximation power with the original scheme. On this concern,
referring to the considered case with uniform breakpoints on $R$, the
following error bound can be obtained for sufficiently smooth functions [24],

$\displaystyle\|f\,-\,\sigma_{f}\|_{\infty,R}\leq\,\|\partial^{p_{1}+1,0}f\|_{\infty,R}\left(\frac{H_{1}}{\nu_{1}-1}\right)^{p_{1}+1}+\|\partial^{0,p_{2}+1}f\|_{\infty,R}\left(\frac{H_{2}}{\nu_{2}-1}\right)^{p_{2}+1}$
(31)
$\displaystyle+\|\partial^{p_{1}+1,p_{2}+1}f\|_{\infty,R}\left(\frac{H_{1}}{\nu_{1}-1}\right)^{p_{1}+1}\left(\frac{H_{2}}{\nu_{2}-1}\right)^{p_{2}+1}\,,$

where $H_{i}$ and $\nu_{i}$ respectively denote the size of $R$ and the number
of uniform breakpoints in the $i$-th coordinate direction, $i=1,2,$ and
$\partial^{i_{1},i_{2}}f$ stands for ${\partial^{i_{1}+i_{2}}f}/{(\partial
t_{1}^{i_{1}}\partial t_{2}^{i_{2}})}.$ Thus the convergence order of the
scheme with respect to $\max\\{H_{1},H_{2}\\}\rightarrow 0$ but also to
$\max\\{(\nu_{1}-1)^{-1},(\nu_{2}-1)^{-1}\\}$ is equal to
$\min\\{p_{1},p_{2}\\}+1$.

###  4.3 B-spline quadrature rules for regular integrals

Since in the assayed BIEs the regular integrands consist of the product
between a bivariate regular function $f$ with a tensor product B-spline
$\hat{B}_{{\bf j},{\bf d}}$, we need to construct a quadrature rule to
evaluate numerically the following integral:

$\int_{R_{\bf j}}f({\bf t})\hat{B}_{{\bf j},{\bf d}}({\bf t})\,d{\bf t}.$ (32)

The first step consists in approximating $f$ with the QI scheme previously
described. In this way we can write:

$\displaystyle f\approx\sigma_{f}=\sum_{{\bf i}\in{\cal J}_{R_{\bf
i}}}\lambda_{\bf i}\hat{B}_{{\bf i},{\bf p}},$

where $\hat{B}_{{\bf i},{\bf p}}$ denotes a tensor product B-spline of bi-
degree ${\bf p}=(p_{1},p_{2})$, which is a trial function of a local spline
space defined in $R_{\bf j}$. Hence

$\displaystyle\int_{R_{\bf j}}f({\bf t})\hat{B}_{{\bf j},{\bf d}}({\bf
t})\,d{\bf t}\approx\int_{R_{\bf j}}\sigma_{f}({\bf t})\hat{B}_{{\bf j},{\bf
d}}({\bf t})\,d{\bf t}.$ (33)

The product $\sigma_{f}\hat{B}_{{\bf j},{\bf d}}$ is a new spline that can be
expressed in terms of another B-spline basis defined in a suitable local
“product” spline space of bi-degree ${\bf p}+{\bf d}$. Generalizing the idea
developed in [27, 28] for the one-dimensional case, the coefficients of
$\sigma_{f}\hat{B}_{{\bf j},{\bf d}}$ in the B-spline basis of the product
space can be written as $G_{{\bf j}}^{\top}\mbox{\boldmath$\Lambda$},$ with
$\Lambda$ defined by QI as specified in (29). In the current bivariate setting
the matrix $G_{\bf j}$ is defined as $G_{\bf
j}:=G_{j_{1}}^{(p_{1},d_{1})}\otimes G_{j_{2}}^{(p_{2},d_{2})},$ where
$G_{j_{1}}^{(p_{1},d_{1})}$ and $G_{j_{2}}^{(p_{2},d_{2})}$ are suitable
direction-wise coefficient matrices and ${\bf j}=(j_{1},j_{2})$. The dimension
of $G_{\bf j}$ is $|{\cal J}_{R_{\bf j}}|\times N_{\pi_{\bf j}}$, where
$N_{\pi_{\bf j}}$ denotes the dimension of the introduced local product spline
space. Thus, the integral in (33) can finally be evaluated as,

$\displaystyle\int_{R_{{\bf j}}}f({\bf t})\hat{B}_{{\bf j},{\bf d}}({\bf
t})d{\bf t}\approx{\bf v}^{\top}G_{\bf j}^{\top}\mbox{\boldmath$\Lambda$}={\bf
v}^{\top}G_{\bf j}^{\top}(\hat{C}^{(p)})^{\top}{\bf f},$

where the matrix $\hat{C}^{(p)}$ is defined in (30) and ${\bf v}$ denotes a
vector of suitable length collecting in lexicographical order the following
integrals of each B-spline $\hat{B}_{{\bf k},{\bf p}+{\bf d}}$ of the basis of
the local product spline space,

$\displaystyle\displaystyle{\int_{R_{\boldsymbol{\pi}}}\hat{B}_{{\bf k},{\bf
p}+{\bf d}}({\bf t})\,d{\bf t}\,=\,\frac{|\mbox{supp}(\hat{B}_{{\bf k},{\bf
p}+{\bf d}})|}{(p_{1}+d_{1}+1)(p_{2}+d_{2}+1)}}\,,\qquad{\bf k}\in{\cal
J}_{\pi_{\bf j}}\,,$

where ${\cal J}_{\pi_{\bf j}}$ is a set of $N_{\pi_{j}}$ distinct multi-
indices, as usually done to identify a tensor-product B-spline basis, and
${R_{\boldsymbol{\pi}}}$ is the support of the spline.

As already mentioned, the derived integration rule is the tensor product
extension of a variant specific for integrals including a B-spline weight of
the 1D quadrature formula introduced in [29]. Therefore, the error analysis
developed in [29] can be easily extended to such variant and lifted to the
tensor product setting. To be more concise, we assume that a uniform
distribution of quadrature nodes is used in each $R_{\bf j}$ (i.e., uniform
breakpoints in the associated local spline space are used by the QI scheme),
since this is our common setting in numerical experiments. Assuming $f$ is a
sufficiently smooth function, the quadrature error of the rule to approximate
the integral in (32) is upper bounded by a positive constant times the
following quantity,

$\displaystyle\frac{H_{1,{\bf j}}H_{2,{\bf
j}}}{(d_{1}+1)(d_{2}+1)}\,\left\\{\,\|\partial^{p_{1}+r_{1}+1,0}f\|_{\infty,R_{\bf
j}}\left(\frac{H_{1,{\bf
j}}}{\nu_{1}-1}\right)^{p_{1}+r_{1}+1}+\|\partial^{0,p_{2}+r_{2}+1}f\|_{\infty,R_{\bf
j}}\left(\frac{H_{2,{\bf j}}}{\nu_{2}-1}\right)^{p_{2}+r_{2}+1}\right.$ (34)
$\displaystyle+\left.\|\partial^{p_{1}+r_{1}+1,p_{2}+r_{2}+1}f\|_{\infty,R_{\bf
j}}\left(\frac{H_{1,{\bf
j}}}{\nu_{1}-1}\right)^{p_{1}+r_{1}+1}\left(\frac{H_{2,{\bf
j}}}{\nu_{2}-1}\right)^{p_{2}+r_{2}+1}\right\\}\,,$

where $H_{i,{\bf j}}$ and $\nu_{i}$ respectively denote the size of $R_{{\bf
j}}$ and the number of uniform quadrature nodes in the $i$-th coordinate
direction, $i=1,2$ and $r_{i}=1(0)$ for $p_{i}$ even (odd).

For fixed numbers $\nu_{i}$ of nodes, the convergence order of the quadrature
scheme with respect to $\max\\{H_{1,{\bf j}},H_{2,{\bf j}}\\}\rightarrow 0$ is

$\displaystyle C_{1}\min\\{H_{1,{\bf j}}^{p_{1}+r_{1}+3},H_{1,{\bf
j}}^{p_{2}+r_{2}+3}\\},$

if the quotient $H_{1,{\bf j}}/H_{2,{\bf j}}$ is bounded from below and above
for all nested integration subdomains $R_{\bf j}$. The constant $C_{1}$
contains norms of all the involved derivatives of $f$ in (34) for the largest
considered $R_{\bf j}$. For fixed $R_{\bf j}$, the convergence order with
respect to $\min\\{\nu_{1},\nu_{2}\\}$ is equal to

$\displaystyle\tilde{C}_{1}\min\\{(\nu_{1}-1)^{-p_{1}-r_{1}-1},(\nu_{2}-1)^{-p_{2}-r_{2}-1}\\},$

where the positive constant $\tilde{C}_{1}$ contains also constants $H_{i,{\bf
j}}$ but not varying $\nu_{i}$.

###  4.4 B-spline quadrature rules for singular integrals

Let us focus now on the rules developed for the numerical approximation of the
second addend on the right of (26), which involves the simplified singular
kernel $U_{\bf s}^{m}({\bf s}-{\bf t}).$ The first step is analogous to that
adopted for regular integrals, since $g$ is approximated with $\sigma_{g}$ by
using the same QI approach and then multiplied by $\hat{B}_{{\bf j},{\bf d}}.$
Thus we can write

$\displaystyle\int_{R_{\bf j}}{U}_{{\bf s}}^{m}({\bf s}-{\bf t})\hat{B}_{{\bf
j},{\bf d}}({\bf t})g({\bf t})\,d{\bf t}\approx\int_{R_{\bf j}}{U}_{{\bf
s}}^{m}({\bf s}-{\bf t})\hat{B}_{{\bf j},{\bf d}}({\bf t})\sigma_{g}({\bf
t})\,d{\bf t}={\bf v}^{\top}G_{\bf j}^{\top}(\hat{C}^{(p)})^{\top}{\bf g}\,,$
(35)

where ${\bf g}$ is the vector whose entries are the values of $g$ used by the
QI operator but now ${\bf v}$ collects in a suitable order all the following
integrals

$\displaystyle\int_{R_{\bf j}}{U}_{{\bf s}}^{m}({\bf s}-{\bf t})\hat{B}_{{\bf
k},{\bf p}+{\bf d}}({\bf t})\ d{\bf t}\,,\qquad{\bf k}\in{\cal J}_{\pi_{\bf
j}}.$ (36)

Thus, in order to complete the definition of the quadrature rule, we are
interested in deriving the exact expression of the above integrals (36). For
this aim we introduce the following notation,

$I^{q_{1},q_{2}}_{r_{1},r_{2}}(k_{1},k_{2}):=\int_{\tau^{(1)}_{k_{1}}}^{\tau^{(1)}_{k_{1}+r_{1}+1}}\int_{\tau^{(2)}_{k_{2}}}^{\tau^{(2)}_{k_{2}+r_{2}+1}}{U}_{{\bf
s}}^{m}({\bf s}-{\bf
t})(t_{1}-s_{1})^{q_{1}}(t_{2}-s_{2})^{q_{2}}\hat{B}_{k_{1},r_{1}}(t_{1})\hat{B}_{k_{2},r_{2}}(t_{2})dt_{1}dt_{2}\,,$

where ${\bf k}=(k_{1},k_{2})$ and
$\tau^{(\ell)}_{k_{\ell}+j},j=0,\ldots,r_{\ell},\ell=1,2$ define the two
univariate knot vectors active in the definition of $\hat{B}_{{\bf k},{\bf
p}+{\bf d}},$ where for brevity we have set $r_{\ell}:=p_{\ell}+d_{\ell}$,
which is the univariate degree in the $\ell$-th direction of the product
space. Then, using the tensor product extension of the Cox-de Boor recurrence
relation for B-splines [3], we can express
$I^{q_{1},q_{2}}_{r_{1},r_{2}}(k_{1},k_{2})$ as a linear combination of
$I^{q_{1}+m_{1},q_{2}+m_{2}}_{r_{1}-1,r_{2}-1}(k_{1}+w_{1},k_{2}+w_{2})$,
where $m_{1},m_{2},w_{1},w_{2}=0,1$. Starting with $q_{1}=q_{2}=0$ and
iterating this recurrence up to B-splines of degree 0 in both coordinate
directions, the procedure can be completed, provided that the following basic
moments have been preliminarily computed for $q_{\ell}$ and $i_{\ell}$,
respectively ranging in the set of indices $0,\ldots,p_{\ell}+r_{\ell}$ and
$k_{\ell},\ldots,k_{\ell}+r_{\ell},\ell=1,2,$

$\displaystyle I^{q_{1},q_{2}}_{0,0}(i_{1},i_{2})$
$\displaystyle=\int_{\tau^{(1)}_{i_{1}}}^{\tau^{(1)}_{i_{1}+1}}\int_{\tau^{(2)}_{i_{2}}}^{\tau^{(2)}_{i_{2}+1}}{U}_{{\bf
s}}^{m}({\bf s}-{\bf
t})(t_{1}-s_{1})^{q_{1}}(t_{2}-s_{2})^{q_{2}}dt_{1}dt_{2}$
$\displaystyle=\int_{s_{1}-\tau^{(1)}_{i_{1}+1}}^{s_{1}-\tau^{(1)}_{i_{1}}}\int_{s_{2}-\tau^{(2)}_{i_{2}+1}}^{s_{2}-\tau^{(2)}_{i_{2}}}{U}_{{\bf
s}}^{m}({\bf z})(-z_{1})^{q_{1}}(-z_{2})^{q_{2}}dz_{1}dz_{2}.$

For this aim we rely on the analytical expressions derived in [19] computed
with the help of Wolfram Mathematica software.

We are left to derive the asymptotic accuracy of these rules with respect to
the $\max\\{H_{1,{\bf j}}\,,H_{2,{\bf j}}\\}\rightarrow 0,$ where, as in the
previous subsection, $H_{i,{\bf j}}$ denotes the size of the integration
domain $R_{\bf j}$ in the $i$-th coordinate direction. Now from (35) we can
easily derive that the absolute value of the related quadrature error can be
upper bounded by

$\displaystyle\|g-\sigma_{g}\|_{\infty,R_{\bf j}}\ \int_{R_{\bf j}}|U_{{\bf
s}}^{m}({\bf s}-{\bf t})|\hat{B}_{{\bf j},{\bf d}}({\bf t})\,d{\bf t}.$ (37)

Using the triangle inequality $|U_{{\bf s}}^{m}({\bf s}-{\bf t})|\leq|U({\bf
s},{\bf t})|+|U_{{\bf s}}^{m}({\bf s}-{\bf t})-U({\bf s},{\bf t})|$, we can
bound the integral in (37) by a sum of two new integrals. Then, by applying
Proposition 2 reported in the Appendix, we can estimate the first integral as

$\displaystyle\ \int_{R_{\bf j}}|U({\bf s},{\bf t})|\hat{B}_{{\bf j},{\bf
d}}({\bf t})\,d{\bf t}\leq\ \int_{R_{\bf j}}|U({\bf s},{\bf t})|\,d{\bf t}\leq
C_{2}\max\\{H_{1},H_{2}\\}.$ (38)

The second integral involves a tail of the series expansion of $U({\bf
s},\bullet)$ about the source point ${\bf s}$. Referring to results in [19],
on this concern we can state that there exists a positive constant $C_{3}$
involving the infinity norm of derivatives of $\tilde{\bf F}$ of degree $\leq
m$ on $R_{\bf j}$, such that

$\displaystyle\|U_{{\bf s}}^{m}({\bf s}-{\bf t})-U({\bf s},{\bf
t})\|_{\infty,R_{\bf j}}\leq C_{3}\frac{\left(\max\\{H_{1,{\bf
j}}\,,\,H_{2,{\bf j}}\\}\right)^{m-1}}{m!},$

hence

$\displaystyle\int_{R_{\bf j}}|U_{{\bf s}}^{m}({\bf s}-{\bf t})-U({\bf s},{\bf
t})|\hat{B}_{{\bf j},{\bf d}}({\bf t})\,d{\bf t}\leq
C_{3}\frac{\left(\max\\{H_{1,{\bf j}}\,,\,H_{2,{\bf
j}}\\}\right)^{m}}{m!(d_{1}+1)(d_{2}+1)}.$ (39)

By combining these estimates with the bound in (31) for the QI error, we
finally derive the following error bound for the singular integral in (35)
with respect to size of $R_{\bf j}$,

$\displaystyle C\max\\{H_{1,{\bf j}}^{p_{1}+2},H_{2,{\bf j}}^{p_{2}+2}\\},$

where $C=C_{1}\max\\{C_{2},C_{3}\\}$ and $C_{1}$ contains norms of all the
involved derivatives of $f=g-\sigma_{g}$ in (34) for the largest considered
$R_{\bf j}$. Instead, by fixing domain $R_{\bf j}$, the error bound with
respect to $\min\\{\nu_{1}-1,\nu_{2}-1\\}$ is equal to

$\displaystyle\tilde{C}\max\\{(\nu_{1}-1)^{-p_{1}-1},(\nu_{2}-1)^{-p_{2}-1}\\},$

where $\tilde{C}$ contains a corresponding constant $\tilde{C}_{1}$ to the one
at the end of Section 4.3, involving the function $f=g-\sigma_{g}$, multiplied
by right-hand side estimates of (38) and (39).

###  4.5 The threshold for near singularity detection

The efficient and accurate evaluation of different types of integrals (23)
appearing in the considered BIEs is crucial, see also [30] on such concern.
For this aim in this subsection we are interested in determining a reasonable
threshold $\delta$ to be used in (25) to switch between nearly singular and
regular integrals. The goal is to apply the error bound (34) introduced in
Section 4.3 to $f=gU,$ in order to derive on each patch a mesh dependent
definition of the threshold $\delta$ to be used in (25).

By analyzing the nearly singular nature of the considered integral in (23), it
is important to observe that the function $f:=gU({\bf s},\bullet)$ for a fixed
${\bf s}=\tilde{\bf F}^{-1}(\tilde{\bf x})$ includes a factor proportional to
$1/\tilde{r}$, where $\tilde{r}=\|\tilde{\bf x}-\tilde{\bf F}(\bullet)\|_{2}$.
Hence, $\|\partial^{i_{1},i_{2}}f\|_{\infty,R_{\bf j}}$ includes a term
$\tilde{r}_{\rm{min},{\bf j}}({\bf x})^{-i_{1}-i_{2}-1},$ where
$\tilde{r}_{\rm{min},{\bf j}}({\bf x})$ is the minimal distance between
$\tilde{\bf x}$ and the physical integration domain image of $R_{\bf j}$ via
the geometry map $\tilde{\bf F},$ see its definition in formula (25). Thus,
considering the expression of the error bound in (34) for a fixed number of
uniformly spaced quadrature nodes, we can infer a remarkable reduction of the
convergence order with respect to $\max\\{H_{1,{\bf j}},H_{2,{\bf
j}}\\}\rightarrow 0$ if $\tilde{r}_{\rm{min},{\bf j}}({\bf x})$ would reduce
proportionally to it, due to increasing terms
“$\|\partial^{i_{1},i_{2}}f\|_{\infty,R_{\bf j}}$”. Instead, by setting

$\displaystyle\tilde{r}_{{\rm min},{\bf j}}({\bf x})\geq\eta_{\bf
j}:=\max\left\\{H_{1,{\bf j}}^{\frac{1}{2(p_{1}+r_{1}+2)}}\,,\,H_{2,{\bf
j}}^{\frac{1}{2(p_{2}+r_{2}+2)}}\right\\}$

this guarantees that for all the scaled integrals treated with the regular
rule in the governing BIEs, the quadrature convergence order is at least equal
to $5/2+\min\\{p_{1}+r_{1},p_{2}+r_{2}\\}.$ In order to have just one
threshold $\delta$ for all the integrals involving the same geometry map, for
the $k$-th patch we set

$\displaystyle\delta=\delta_{k}:=\max_{{\bf j}\in{\cal I}_{k}}\eta_{\bf j}.$

##  5 Numerical examples

In this section we test our model to numerically solve one interior and two
exterior Helmholtz problems, all with Neumann boundary conditions, which are
more commonly considered for applications.

The two well-known benchmarks for exterior problems are the pulsating sphere
and the rigid scattering on a sphere; see for instance [14, 15, 16]. The
sphere has radius $R=1$ and it is parameterized by $6$ quartic (${\bf
d}_{g}=(4,4)$) NURBS patches, see Figure 2(a) and refer to [31] for details.
As in [32], the considered quartic NURBS parameterization based on cube
topology is the preferred choice to the 8-patch tiling in order to avoid
singularities in the geometry description. The last example is a standard
acoustic problem interior to a torus; see for example [14, 15]. The toroidal
surface $\Gamma$ has inner and outer radii respectively equal to $1$ and $3$
and is exactly represented in a $16$-patch quadratic (${\bf d}_{g}=(2,2)$)
NURBS form, see Figure 2(b) and refer to [12] for details.

(a) The $6$-patch quartic NURBS representation of a sphere. (b) The $16$-patch
quadratic NURBS representation of a torus.

Figure 2: The two NURBS multipatch representations of a sphere and a torus.

For the discretization space $V$, on each patch we employ standard tensor
product B-splines of different bi-degrees $(d,d)$ and of highest continuity;
in the last two experiments we enforce $C^{0}$ continuity for the basis
functions across patch interfaces. The quality of the numerical results is
evaluated by considering the behavior of the on-surface relative error
measured in the $L^{2}(\Gamma)$ norm,

$e_{L^{2}}=\frac{\|u_{h}-u^{ex}\|_{L^{2}(\Gamma)}}{\|u^{ex}\|_{L^{2}(\Gamma)}}.$

###  5.1 Pulsating sphere

The following problem serves as a benchmark test for the developed quadrature
rules; it is also referred to as a patch test for IgA-BEM. The model can be
used for example to numerically compute the sound pressure at a distance $r$
from the center of the sphere for a constant speed of wave in the radial
direction. For the Neumann boundary conditions set to $u_{\rm N}\equiv
e^{i\kappa}(i\kappa-1)/(4\pi)$ for ${\bf x}\in\Gamma$, the analytic solution
is given as $u=e^{i\kappa r}/(4\pi r)$, hence the missing Cauchy data is
$\phi\equiv e^{i\kappa}/(4\pi)$ for ${\bf x}\in\Gamma$.

Since the geometry representation of the sphere is exact and $\phi$ is
constant over the entire boundary, the error can be attributed solely to
numerical integration. On each patch the knot vectors describing the local
part of the discretization space $\hat{\cal S}_{h^{(\ell)}}^{(\ell)}$ are

$\displaystyle
T_{1}^{(\ell)}=T_{2}^{(\ell)}=\begin{bmatrix}0&0.25&0.5&0.75&1\end{bmatrix}.$

The collocation points are midpoints of the knots, mapped to the physical
domain via geometry map. Fixing $\kappa=1,{\bf d}=(0,0)$ and setting $c=0.25$
for completing the definition of the threshold used for near singularity
detection, we measure the $L^{2}$ error of the approximate solution with
respect to different number of quadrature nodes on the support of B-splines.
The number $n_{q}$ of quadrature nodes increases in each coordinate direction
with $n_{q}=12\alpha+1$, for $\alpha=1,2,\dots,6$. The same quadrature nodes
are used for singular, regularized and regular integrals.

In Figure 3(a) we observe the error distribution (modulus of the difference
between the exact and approximate solution) for $\alpha=1$, $m=2$ (number of
terms in the singularity extraction) , $p=2$ (spline degree of the QI
operator). The error is well evenly spread on the boundary with a slight
increase around patch corners, where the kernel expansion is more critical.
The error is also slightly higher in the interior of the patches than in the
neighborhood of the edges since the area of interior mapped cells is slightly
larger than the boundary ones. In Figure 3(b) we observe that error
convergence order is impacted by the choice of $m$ and $p$. For $p=m=2$ the
obtained order of convergence is 3. If we increase both $m$ and $p$ to 3, the
order of convergence increases to 4.

(a) Error distribution on the boundary for $\alpha=6,m=2,p=2$. (b) Error
convergence plot with respect to multiplier $\alpha$ for quadrature number.

Figure 3: Pulsating sphere for $\kappa=1$.

The comparison can be done with Figure 6(a) of [15], where $R\kappa=1$ as we
have used in our experiments. In the referred picture the results obtained
with new approach (termed as New) are compared with those produced by the
method introduced in [14] (termed as Old). The relative $L^{2}$ error,
$e_{L^{2}}$, is shown against the total number of quadrature nodes $n_{qp}$
(where the -negligible- number of quadrature nodes in the elements with active
singularity is not taken into account). In our case the total number of
quadrature nodes, for $\alpha=6$, is $6\ \cdot 16\ \cdot 72^{2}\approx 5.0\
\cdot 10^{5}$ achieving an error (with $p=3$) of order $10^{-11}.$ This is a
much better result than that produced by the Old approach but also better than
that given by the New one.

###  5.2 Rigid Scattering on a sphere

In rigid scattering (exterior) problems reformulated in the frequency domain,
it is assumed that an acoustic pressure $p_{inc}$ of amplitude $A$ is produced
by a wave vector $\kappa{\bf v},$ with ${\bf v}$ denoting a unit vector
prescribing the direction of the wave,

$p_{inc}({\bf x}):=Ae^{i\kappa({\bf v}\cdot{\bf x})}.$

Such pressure is incident on a rigid body represented by a volume
$\Omega^{(i)}$ with boundary $\Gamma,$ where the rigidity assumption implies
that it reacts producing in $\Omega^{(e)}$ an additional scattered pressure
field $p$ verifying the Helmholtz equation in $\Omega^{(e)}\,,$ the
Sommerfield radiation condition at infinity and having a variation in the
normal direction on $\Gamma$ opposite to that of $p_{inc}$, see for example
[15]. The scattered pressure $p$ is the unknown to be determined in
$\Omega^{(e)}$ and it can be written as the difference between a total
pressure field $p_{tot}$ and $p_{inc},$ where $p_{tot}$ is such that

$\left\\{\begin{array}[]{ll}p_{tot}({\bf x})&=p_{inc}({\bf
x})+(V_{\kappa}\partial_{n}p_{tot})-(K_{\kappa}p_{tot})({\bf x})\,,\,\,{\bf
x}\in\Omega^{(e)}\cr\partial_{n}p_{tot}({\bf x})&=0\,,\,\,\forall{\bf
x}\in\Gamma\,.\cr\end{array}\right.$

Note that in this experiment we assume $\Gamma$ as a sphere with radius $R=1$
(again parameterized as previously described) and so, for symmetry, without
loss of generality we can assume ${\bf v}=(1,0,0).$ The problem can be used as
a benchmark because the analytic expression of $p$ is a priori known, see for
example [16],

$p({\bf x})=-A\ \sum_{n=0}^{\infty}\frac{i^{n}(2n+1)j_{n}^{\prime}(\kappa
R)}{h_{n}^{\prime}(\kappa R)}\ P_{n}(\cos(\theta))\ h_{n}(\kappa r)\,,$

where in our experiments we set $A=1$ and $n$ ranges from $1$ to $10$ because
this is sufficient to study the convergence of our numerical scheme for the
used frequencies. In this formula, assuming the sphere centered in ${\bf O}$,
$r=\|{\bf x}-{\bf O}\|_{2},$ $\theta$ is the angle between the vector ${\bf
x}-{\bf O}$ and ${\bf v},$ $h_{n}$ and $j_{n}$ are the spherical Hankel and
Bessel functions of order $n$ (the derivative with respect to its argument is
denoted by ′) and $P_{n}$ is the Legendre polynomial of order $n.$

(a) The real part of the total pressure on the sphere (b) The imaginary part
of the total pressure on the sphere (c) The real part of the total pressure on
the equatorial plane (d) The imaginary part of the total pressure on the
equatorial plane

Figure 4: Rigid scattering on a sphere with $\kappa=2$: on and off surface
total pressure distributions.

(a) On surface error distribution ($d=2\,,n=5$) (b) On surface error
distribution ($d=4\,,n=3$)

Figure 5: Rigid scattering on a sphere with $\kappa=2$: on surface
distribution of the modulus of the absolute error ($N_{\rm DOF}=218$).

The values considered for $\kappa$ in the experiments are $\kappa=1,2,3$ which
in the air corresponds to a frequency of about $50,100,150$ Hz. In the
implementation we aim to approximate the unknown Dirichlet Cauchy datum
$\phi=p_{tot}|\Gamma$ (alternatively also denoted below just with $u$ which
denotes also its extension to $\Omega^{(e)}$). Thus, considering that the
corresponding given Neumann datum is homogeneous, the boundary integral
equation to be considered simplifies to the following one,

$(K_{\kappa}+\frac{1}{2}I)\phi({\bf x})\,=\,p_{inc}({\bf x})\,,\quad{\bf
x}\in\Gamma\,.$ (40)

The discretization of this BIE has been done by using a globally $C^{0}$
multi-patch IgA space with the same size of uniform elements on each patch.
For this experiment we consider the following simplification

$\displaystyle\frac{{\bf r}\cdot{\bf n}_{\bf y}}{r^{2}}=-\frac{1}{2R},$

where ${\bf n}$ points toward the center of the sphere since the considered
problems are exterior. Then by taking into account (8), for the imaginary part
of the double layer kernel $K_{\kappa}$ we use regular rules while, for its
real part we can rely on the singular rules for the kernel $1/r.$ The constant
$c$ involved in the near singularity detection phase in this experiment is
chosen to be equal to $0.1$. In all the tests reported for this and also in
the next experiment for all the three kinds of quadrature rules the QI degree
is $2$ for the singular and regularized rules and it is $4$ for the regular
ones. We observe that a low degree is reasonable for regularized integrals,
since higher degrees of the quasi interpolating spline are profitable only if
the function it approximates is sufficiently regular. For singular integrals a
low degree has been preferred to reduce the cost of the related weight
computation. For all the involved quadrature rules, in this experiment the
chosen number of QI nodes (uniformly distributed in the support of each trial
B-spline function) is $7,9,11,$ respectively for degrees
$d=d_{1}=d_{2}=2,3,4,$ on both coordinate directions. Note that for any degree
this corresponds to select the nodes at the knots of the trial B-spline and at
their midpoints in each coordinate.

(a) Modulus of the exact and of the numerical solution on ${\cal C}$ for
$\kappa=1.$ (b) Modulus of the exact and of the numerical solution on ${\cal
C}$ for $\kappa=3.$ (c) Distribution of the error $e_{P}$ on ${\cal C}$ for
$\kappa=1.$ (d) Distribution of the error $e_{P}$ on ${\cal C}$ for
$\kappa=3.$

Figure 6: Rigid scattering on a sphere: results on the exterior equatorial
circumference ${\cal C}$ ($d=2,p=2,2,4$).

In order to underline the oscillating nature of the solution of (40) and of
the related domain solution, Figure 4 shows the real and the imaginary part of
the exact Cauchy datum $\phi=p_{tot}|\Gamma$ for $\kappa=2$, (frequency of
about $100$ Hz) together with the corresponding real and imaginary parts of
the total pressure in an area of the equatorial plane exterior to the sphere.
Figure 5 shows the distribution on the sphere of the modulus of the absolute
error obtained with both the choices $d=2$ and $d=4.$ Note that the number of
degrees of freedom ($N_{\rm DOF}=218$) is the same on the left and on the
right of this figure, since the number $n^{2}$ of elements in each patch is
chosen in order to ensure the fulfillment of such goal, that is $n=5$ for
$d=2$ and $n=3$ for $d=4$. This means that when $d=4$, the elements are larger
than in the other case and so an analogous maximal error on the surface is
acceptable. Note also that the error distribution on the sphere implies a
slightly lower $L^{2}(\Gamma)$ relative error on surface for $d=4$ (it is
$3.20e-03$ for $d=2$ and $2.43e-03$ for $d=4$). Our results can be compared
with those, less accurate, shown in Figure 13 of [16], which refers to
$\kappa=2$, $R=1$, $d=4$ as well, and are obtained with $N_{\rm DOF}=200$.

The off surface quality of our approach is clarified by the plots reported in
Figure 6, as it is done in [14]. Focusing on the approximation of the modulus
of the solution, the figure shows the results obtained on the exterior
equatorial circumference ${\cal C}$ with radius $10$ setting $d=2$ (the other
parameters are chosen as already described), where the reported error is the
following one

$e_{P}({\bf x}):=\frac{|\,|u_{h}({\bf x})|-|u^{ex}({\bf x})|\,|}{|u^{ex}({\bf
x})|}.$ (41)

The figure outlines a good conformability between the moduli of the exact and
of the numerical solution obtained with very few degrees of freedom for both
$\kappa=1$ and $\kappa=3.$ The good quality of our exterior reconstructions
can be further checked looking at the distribution on ${\cal C}$ of the error
$e_{P}.$ Figures 6(c) and 6(d) can be compared respectively with Figure 17
(where $\kappa=2$, but $R=0.5$) and Figure 20 (where $\kappa=6$, but $R=0.5$)
of [14], confirming the better performance of our approach. We should remark
however, as underlined in [15], that the accuracy of the obtained results in
[14] is probably also influenced by the adopted singular parameterization of
the sphere.

(a) $\kappa=1$ (b) $\kappa=2$ (c) $\kappa=3$

Figure 7: Rigid scattering on a sphere: behavior of the $L^{2}(\Gamma)$ error,
$e_{L^{2}}$, versus the number $n$ of elements in each coordinate direction
for each patch.

In order to complete the analysis of the results for this experiment, we
present also Figure 7, which shows the dependency from $n$ of the
$L^{2}(\Gamma)$ relative error, $e_{L^{2}}$ for $d$ ranging between $2$ and
$4,$ confirming its nice convergence behavior for all the considered values of
$\kappa$. We remark that these results confirm a better performance than the
approach used in [15], compare the plot corresponding to $d=4$ in our Figure
7(a) with the plot related to Parameterization 2 (our 6 patch parameterization
of the sphere), IGABEM CCBIE (Conventional BIE, Collocation), shown in Figure
12 of [15]. For a precise comparison we note that $n=1,2,\ldots,7$ in our case
corresponds to $N_{\rm DOF}=98,152,218,296,386,488,602$.

###  5.3 Acoustic problem interior to a torus

This last example is taken from [14]. The Helmholtz equation is considered
with exact acoustic potential chosen to be

$\displaystyle\phi({\bf x})=\sin(\kappa x/\sqrt{3})\sin(\kappa
y/\sqrt{3})\sin(\kappa z/\sqrt{3}).$ (42)

The Neumann datum is prescribed by computing the acoustic velocity field as
$({{\partial\phi}/{\partial{\bf n}}})({\bf x})$ for ${\bf x}\in\Gamma$, where
${\bf n}$ denotes the outward unit normal vector to the surface. In the
experiment we consider the wavenumber $\kappa=2$ and global $C^{0}$
discretization space with spline bi-degree $(d,d)$ for $d=2,3$. The constant
$c$ used in (25) for near singularity detection is set to $0.1$. Regarding the
adopted quadrature scheme, we use QI degree $4$ for regular integrals and
degree $2$ for singular and regularized integrands. For $d=2$ we set $m=2$
(number of terms in the singularity extraction) and the number of uniformly
spaced quadrature nodes in each parametric direction on the support of a
B-spline is set to $7$; hence nodes are just knots and knot midpoints of
B-splines. When $d=m=3,$ we use $9$ quadrature nodes for singular and regular
integrals; for regularized integrals we use $13$ nodes. To obtain a more even
knot span in both parametric directions of the physical space, for all 16
patches the starting mesh is described by the following knot vectors,

$\displaystyle
T_{1}^{(\ell)}=\begin{bmatrix}0&1/9&2/9&3/9&4/9&5/9&6/9&7/9&8/9&1\end{bmatrix},\qquad
T_{2}^{(\ell)}=\begin{bmatrix}0&1/3&2/3&1\end{bmatrix}.$

While the imaginary part of the acoustic potential $\phi$ is zero on the
entire boundary $\Gamma$, the chosen wavenumber produces several oscillation
in the real part of $\phi$, see Figure 8(a). The error $e_{L^{2}}$ is plotted
with respect to size $h$ in Figure 8(b); a convergence rate equal to $d+1$ can
be observed. In Figure 8(c) and (d) the absolute error for the real part and
the imaginary part of the unknown acoustic potential $\phi$ is shown,
respectively.

(a) Acoustic potential. (b) Error convergence against mesh size $h$.

(c) On surface distribution of the real part of the absolute error for $d=3$
and $h=1/9$ ($N_{\rm DOF}=640$). (d) On surface distribution of the imaginary
part of the absolute error for $d=3$ and $h=1/9$ ($N_{\rm DOF}=640$).

Figure 8: Acoustic problem with $\kappa=2$, interior to a torus.

##  6 Conclusion

In this paper we address 3D Helmholtz problems using isogeometric BEMs and
relying on conformal multi-patch smooth geometries and spline discretization
spaces. The governing boundary integral equations are numerically approximated
by a collocation scheme. The adopted quadrature rules for both regular and
singular integrals are tailored for B-splines and they allow a profitable
function-by-function matrix assembly. A procedure is proposed for the
automatic detection of near singular integrals. Besides studying the
approximation power of the considered rules, an arrangement with a singularity
extraction technique is proposed, extending it to the multi-patch setting for
the first time. Numerical results confirm the effectiveness of the method by
achieving sufficient accuracy of the numerical solution with a small number of
uniformly distributed quadrature nodes.

A possible future work includes a generalization of the proposed simulation
model to more flexible non-conformal discretization spaces and
parameterizations. Further testing of the proposed methodology in both the
conformal and non-conformal setting for more general types of contacts between
patches could be also taken into account.

##  7 Appendix

To provide a self-contained analysis of the developed integration procedures,
a more technical part related to the behavior of the considered singular
kernels is given in this appendix.

Let $\Gamma$ be a regular surface admitting a $C^{2}$ regular parameterization
on $D:=[0\,,\,1]^{2}\,,\,\,{\bf F}:D\rightarrow\Gamma,\,\,{\bf
t}=(t_{1},t_{2})\mapsto{\bf F}({\bf t})$. Furthermore, for any fixed ${\bf
s}\in D,$ let us introduce the following notation,

${\bf r}_{\bf s}({\bf t}):={\bf F}({\bf t})-{\bf F}({\bf s})\,,\,\,r_{\bf
s}({\bf t}):=\|{\bf r}_{\bf s}({\bf t})\|_{2}\,,\,\,\mbox{\boldmath$\nu$}({\bf
t}):=\frac{\partial{\bf F}}{\partial t_{1}}({\bf t})\times\frac{\partial{\bf
F}}{\partial t_{2}}({\bf t})\,.$

###### Proposition 1

For each fixed ${\bf s}\in D$ and for any unit vector ${\bf
u}=(\cos\theta,\sin\theta)$ such that ${\bf s}+\rho{\bf u}$ is included in $D$
for $\rho$ positive and sufficiently small, it is

$\lim_{\rho\rightarrow 0^{+}}\frac{\rho}{r_{\bf s}({\bf s}+\rho{\bf
u})}=\frac{1}{\sqrt{E\cos^{2}\theta+2F\sin\theta\cos\theta+G\sin^{2}\theta}},$
(43) $\lim_{\rho\rightarrow 0^{+}}\frac{{\bf r}_{\bf s}({\bf s}+\rho{\bf
u})\cdot\mbox{\boldmath$\nu$}({\bf s}+\rho{\bf u})}{r^{2}({\bf s}+\rho{\bf
u})}=-\frac{1}{2}\frac{L\cos^{2}\theta+2M\sin\theta\cos\theta+N\sin^{2}\theta}{E\cos^{2}\theta+2F\sin\theta\cos\theta+G\sin^{2}\theta},$
(44)

which are both finite limits, with $E,F,G$ and $L,M,N$ respectively denoting
the coefficients of the first and second fundamental forms (with respect to
the non-normalized normal $\nu$) of ${\bf F}$ evaluated at ${\bf s}$,

$\displaystyle E:=\frac{\partial{\bf F}}{\partial t_{1}}({\bf
s})\cdot\frac{\partial{\bf F}}{\partial t_{1}}({\bf s}),$ $\displaystyle
F:=\frac{\partial{\bf F}}{\partial t_{1}}({\bf s})\cdot\frac{\partial{\bf
F}}{\partial t_{2}}({\bf s}),$ $\displaystyle G:=\frac{\partial{\bf
F}}{\partial t_{2}}({\bf s})\cdot\frac{\partial{\bf F}}{\partial t_{2}}({\bf
s}),$ (45) $\displaystyle L:=\mbox{\boldmath$\nu$}({\bf
s})\cdot\frac{\partial^{2}{\bf F}}{\partial t_{1}^{2}}({\bf s}),$
$\displaystyle M:=\mbox{\boldmath$\nu$}({\bf s})\cdot\frac{\partial^{2}{\bf
F}}{\partial t_{1}\partial t_{2}}({\bf s}),$ $\displaystyle
N:=\mbox{\boldmath$\nu$}({\bf s})\cdot\frac{\partial^{2}{\bf F}}{\partial
t_{2}^{2}}({\bf s}).$ (46)

Proof. Setting for brevity ${\bf F}_{i}:=\frac{\partial{\bf F}}{\partial
t_{i}}({\bf s}),\,\,i=1,2$ and ${\bf F}_{i,j}:=\frac{\partial^{2}{\bf
F}}{\partial t_{i}\partial t_{j}}({\bf s}),\,i,j=1,2,$ and using the short
notation ${\bf r}={\bf r}({\bf s}+\rho{\bf u}),r=r({\bf s}+\rho{\bf
u}),\mbox{\boldmath$\nu$}=\mbox{\boldmath$\nu$}({\bf s}+\rho{\bf u}),$ we can
write,

${\bf r}=\rho({\bf F}_{1}\cos\theta+{\bf
F}_{2}\sin\theta)+\rho^{2}(\frac{1}{2}{\bf F}_{1,1}\cos^{2}\theta+{\bf
F}_{1,2}\cos\theta\sin\theta+\frac{1}{2}{\bf
F}_{2,2}\sin^{2}\theta)+O(\rho^{3})$

and

$\begin{array}[]{ll}\mbox{\boldmath$\nu$}&=[{\bf F}_{1}+\rho\big{(}{\bf
F}_{1,1}\cos\theta+{\bf F}_{1,2}\sin\theta\big{)}+O(\rho^{2})]\times[{\bf
F}_{2}+\rho({\bf F}_{1,2}\cos\theta+{\bf
F}_{2,2}\sin\theta)+O(\rho^{2})]=\cr&=({\bf F}_{1}\times{\bf
F}_{2})+\rho[({\bf F}_{1}\times{\bf F}_{1,2}-{\bf F}_{2}\times{\bf
F}_{1,1})\cos\theta+({\bf F}_{1}\times{\bf F}_{2,2}-{\bf F}_{2}\times{\bf
F}_{1,2})\sin\theta]+O(\rho^{2})\,.\end{array}$

Since
$r^{2}=\rho^{2}(E\cos^{2}\theta+2F\sin\theta\cos\theta+G\sin^{2}\theta)+O(\rho^{3})\,,$
we have,

$\frac{\rho}{r}=\frac{1}{\sqrt{E\cos^{2}\theta+2F\sin\theta\cos\theta+G\sin^{2}\theta+O(\rho)}}\,,\quad{\bf
r}\cdot\mbox{\boldmath$\nu$}=-\frac{1}{2}\rho^{2}\big{(}L\cos^{2}\theta+2M\sin\theta\cos\theta+N\sin^{2}\theta\big{)}+O(\rho^{3})\,.$

In order to confirm that the two limits in (43) and (44) are always finite, we
observe that their denominator can be respectively written as $\sqrt{{\bf
u}^{T}M_{\bf s}{\bf u}}$ and ${\bf u}^{T}M_{\bf s}{\bf u}\,,$ with $M_{\bf s}$
denoting the following symmetric $2\times 2$ positive definite matrix
associated to the first fundamental form of ${\bf F}$ evaluated in ${\bf s},$

$M_{\bf s}:=\left(\begin{array}[]{ll}E&F\cr F&G\end{array}\right)\,.$ (47)

Now the proof is complete.

Note that the previous result implies that both the functions ${\bf
t}\rightarrow\|{\bf t}-{\bf s}\|_{2}/r_{\bf s}({\bf t})$ and ${\bf
t}\rightarrow({\bf r}_{\bf s}({\bf t})\cdot\nu({\bf t}))/r_{\bf s}({\bf
t})^{2}$ are discontinuous at ${\bf t}={\bf s}.$ Despite this, the following
corollary allows one to verify their uniform boundedness in $D$.

###### Corollary 1

Under the assumptions of the previous proposition, both the functions
$g_{1,{\bf s}}({\bf t}):=\|{\bf t}-{\bf s}\|_{2}/r_{\bf s}({\bf t})$ and
$g_{2,{\bf s}}({\bf t}):=({\bf r}_{\bf s}({\bf t})\cdot\nu({\bf t}))/r_{\bf
s}^{2}({\bf t})$ are bounded in their domain $D\setminus\\{{\bf s}\\},$
uniformly with respect to ${\bf s}\in D.$

Proof. Since the proof is analogous for $g_{1,{\bf s}}$ and $g_{2,{\bf s}},$
let us refer only to $g_{1,{\bf s}}$.

Using proof by contradiction, we can assume that there exists a convergent
sequence $\\{({\bf s}_{k_{1}},{\bf t}_{k_{1}})\\}$ in $D^{2}$ such that
$g_{1,{\bf s}_{k_{1}}}({\bf t}_{k_{1}})\rightarrow+\infty$. Then, considering
that ${\bf F}$ is continuous and injective, it holds $\rho_{k_{1}}:=\|{\bf
t}_{k_{1}}-{\bf s}_{k_{1}}\|_{2}\rightarrow 0$ and the limit of the convergent
sequence is equal to $({\bf s},{\bf s})$ for some ${\bf s}\in D$. Now, setting
${\bf t}_{k_{1}}={\bf
s}_{k_{1}}+\rho_{k_{1}}(\cos\theta_{k_{1}},\sin\theta_{k_{1}}),$ and by
considering that $\theta_{k_{1}}$ varies in the periodic interval $[0,2\pi]$,
we can find a subsequence $\\{({\bf s}_{k_{2}},{\bf t}_{k_{2}})\\}$ of
$\\{({\bf s}_{k_{1}},{\bf t}_{k_{1}})\\}$ such that $\\{\theta_{k_{2}}\\}$
tends to a certain angle $\theta$. Then, using arguments analogous to the ones
used in the proof of the previous proposition, it can be proven that

$\displaystyle\lim_{k_{2}\rightarrow\infty}g_{1,{\bf s}_{k_{2}}}({\bf
t}_{k_{2}})=\frac{1}{\sqrt{E\cos^{2}\theta+2F\sin\theta\cos\theta+G\sin^{2}\theta}},$
(48)

where $E,F,G$ are defined in (45).

Since the limit in (48) exists, we obtain a contradiction, since $g_{1,{\bf
s}_{k_{3}}}({\bf t}_{k_{3}})\rightarrow+\infty$ for every subsequence
$\\{({\bf s}_{k_{3}},{\bf t}_{k_{3}})\\}$ of sequence $\\{({\bf
s}_{k_{1}},{\bf t}_{k_{1}})\\}$. Since ${\bf F}$ is regular, we can bound
$E,F,G$ for all ${\bf s}\in D$ and obtain a common upper bound for the limit
(48) for all ${\bf s}\in D$.

The following last proposition is useful to estimate an upper bound for our
quadrature rules for singular integrals.

###### Proposition 2

Let $R\subset D$ be a rectangular integration domain with edges of dimensions
$H_{1}$ and $H_{2}.$ Then there exists a constant $C$ such that for each ${\bf
s}\in R$ it is

$\int_{R}|U({\bf s},{\bf t})|d{\bf t}\,\,\leq C\max\\{H_{1},H_{2}\\}\,,$ (49)

where $U$ is one of the two kernels $U_{\rm SL}$ and $U_{\rm DL}$ defined in
(24).

Proof. First of all, setting $\rho_{\bf s}({\bf t}):=\|{\bf t}-{\bf s}\|_{2},$
we observe that we can write $U_{\rm SL}({\bf s},{\bf t})=g_{1,{\bf s}}({\bf
t})/\rho_{\bf s}({\bf t})$ and $U_{\rm DL}({\bf s},{\bf t})=g_{1,{\bf s}}({\bf
t})g_{2,{\bf s}}({\bf t})/\rho_{\bf s}({\bf t})\,,$ where $g_{1,{\bf s}}$ and
$g_{2,{\bf s}}$ are the two functions defined in the previous corollary. Thus
we can derive from the corollary itself that, for both the considered
definitions of $U,$ there exists a positive constant $Q$ such that

$|U({\bf s},{\bf t})|\leq Q/\rho_{\bf s}({\bf t})\,,\,\,\forall({\bf s},{\bf
t})\in D^{2},\,\,{\bf s}\neq{\bf t}\,.$

Then, denoting with ${\cal A}_{\bf s}$ the circle centered at ${\bf s}$ with
radius $\sqrt{2}\max\\{H_{1},H_{2}\\}$ we have that $R\subset D\cap A_{\bf s}$
which implies the following inequalities,

$\int_{R}|U({\bf s},{\bf t})|d{\bf t}\,\,\leq\int_{D\cap{\cal A}_{\bf
s}}|U({\bf s},{\bf t})|d{\bf t}\leq Q\int_{{\cal A}_{\bf s}}\frac{1}{\rho_{\bf
s}({\bf t})}d{\bf t}\,.$

Thus, using polar coordinates centered at ${\bf s}$ the thesis easily
descends.

The constant $C$ in (49) depends on ${\bf F}$. By prescribing common bounds
for coefficients in (45) and (46) for all “well-behaved” functions ${\bf F}$
considered in practical applications, a common constant $C$ can be
sufficiently well estimated.

## Acknowledgements

Tadej Kanduč was partially funded for developing this research by SIMAI-ACRI,
which is gratefully acknowledged. The research of Antonella Falini was founded
by PON Project AIM 1852414 CUP H95G18000120006 ATT1. The Italian authors are
members of the INdAM Research group GNCS. The INdAM support through GNCS and
Finanziamenti Premiali SUNRISE is gratefully acknowledged.

## References

  * [1] J. Cottrell, T. Hughes, Y. Bazilevs, Isogeometric analysis: toward integration of CAD and FEA, John Wiley & Sons, 2009.
  * [2] G. Farin, J. Hosheck, M. S. Kim, Handbook of Computer Aided Geometric Design, Elsevier, Amsterdam, 2002.
  * [3] C. de Boor, A practical guide to splines, revised Edition, Vol. 27 of Applied Mathematical Sciences, Springer-Verlag, New York, 2001.
  * [4] L. Schumaker, Spline functions: basic theory, Cambridge University Press, 2007.
  * [5] E. Sande, C. Manni, H. Speleers, Explicit error estimates for spline approximation of arbitrary smoothness in isogeometric analysis, Numer Math 144 (2020) 889–929.
  * [6] C. Politis, A. I. Ginnis, P. D. Kaklis, K. Belibassakis, C. Feurer, An isogeometric bem for exterior potential-flow problems in the plane, in: 2009 SIAM/ACM Joint Conference on Geometric and Physical Modeling, 2009, pp. 349–354.
  * [7] R. Simpson, S. Bordas, J. Trevelyan, T. Rabczuk, A two-dimensional isogeometric boundary element method for elastostatic analysis, Comput Methods Appl Mech Engrg 209–212 (2012) 87–100.
  * [8] A. Aimi, M. Diligenti, M. Sampoli, A. Sestini, Isogeometric Analysis and Symmetric Galerkin BEM: a 2D numerical study, Appl Math Comp 272 (2016) 173–186.
  * [9] G. Beer, B. Marussig, C. Duenser, The Isogeometric Boundary Element Method, Springer, 2020.
  * [10] Q. Zang, J. Liu, W. Ye, G. Lin, Isogeometric boundary element for analyzing steady-state heat conduction problems under spatially varying conductivity and internal heat source, Comput Math Appl 80 (7) (2020) 1767–1792.
  * [11] G. Gantner, D. Praetorius, Adaptive BEM for elliptic PDE systems, part II: Isogeometric analysis with hierarchical B-splines for weakly-singular integral equations, Comput Math Appl 117 (2022) 74–96.
  * [12] A. Falini, C. Giannelli, T. Kanduč, M. L. Sampoli, A. Sestini, A collocation IGA-BEM for 3D potential problems on unbounded domains, in: Geometric Challenges in Isogeometric Analysis, Springer, 2022, pp. 31–47.
  * [13] S. Kirkup, The boundary element method in acoustics: a survey, Appl Sci 9 (2019) 1642.
  * [14] R. Simpson, M. Scott, M. Taus, D. Thomas, H. Lian, Acoustic isogeometric boundary element analysis, Comput Methods Appl Mech Engrg 269 (2014) 265–290.
  * [15] J. Venas, T. Kvamsdal, Isogeometric boundary element method for acustic scattering by a submarine, Comput Methods Appl Mech Engrg 359 (2020) 112670.
  * [16] A. Shaaban, C. Anitescu, E. Atroshchenko, T. Rabczuk, 3D isogeometric boundary element analysis and structural shape optimization for Helmholtz acoustic scattering problems, Comput Methods Appl Mech Engrg 384 (2021) 113950.
  * [17] A. Falini, T. Kanduč, M. Sampoli, A. Sestini, Cubature rules based on bivariate spline quasi-interpolation for weakly singular integrals, in: G. Fasshauer, M. Neamtu, L. Schumaker (Eds.), Approximation Theory XVI: Nashville 2019, Springer, 2021, pp. 73–86.
  * [18] A. Aimi, F. Calabrò, M. Diligenti, M. Sampoli, G. Sangalli, A. Sestini, Efficient assembly based on B-spline tailored quadrature rules for the IgA-SGBEM, Comput Methods Appl Mech Engrg 331 (2018) 327–342.
  * [19] T. Kanduč, Isoparametric singularity extraction technique for 3D potential problems in BEM, Comput Methods Appl Mech Engrg 398 (2022) 115271.
  * [20] D. Grebenkov, D. Nguyen, Geometrical structure of Laplacian eigenfunctions, SIAM Review 55 (4) (2013) 601–667.
  * [21] J. Venas, T. Jenserud, Exact 3D scattering solutions for spherical symmetric scatterers, J Sound Vib 440 (2019) 439–479.
  * [22] K. Atkinson, The numerical solution of integral equations of the second kind, Cambridge Monographs on Applied and Computational Mathematics, Cambridge University Press, New York, 1997.
  * [23] Y. Wang, D. Benson, Multi-patch nonsingular isogeometric boundary element analysis in 3D, Comput Methods Appl Mech Engrg 293 (2015) 71–91.
  * [24] C. Bracco, C. Giannelli, F. Mazzia, A. Sestini, Bivariate hierarchical Hermite spline quasi-interpolation, BIT Numer Math 56 (4) (2016) 1165–1188.
  * [25] F. Mazzia, A. Sestini, The BS class of Hermite spline quasi-interpolants on nonuniform knot distributions, BIT 49 (3) (2009) 611–628.
  * [26] F. Calabrò, A. Falini, M. Sampoli, A. Sestini, Efficient quadrature rules based on spline quasi-interpolation for application to IgA-BEMs, J Comput Appl Math 338 (2018) 153–167.
  * [27] K. Mørken, Some identities for products and degree raising of splines, Contr Approx 7 (1991) 195–208.
  * [28] A. Falini, C. Giannelli, T. Kanduč, M. L. Sampoli, A. Sestini, An adaptive IgA-BEM with hierarchical B-splines based on quasi-interpolation quadrature schemes, Int J Numer Methods Eng 117 (10) (2019) 1038–1058.
  * [29] F. Mazzia, A. Sestini, Quadrature formulas descending from BS Hermite spline quasi-interpolation, J Comput Appl Math 236 (2012) 4105–4118.
  * [30] L. af Klinteberg, C. Sorgentone, A. Tornberg, Quadrature error estimates for layer potentials evaluated near curved surfaces in three dimensions, Comput Math Appl 111 (2022) 1–19.
  * [31] J. Cobb, Tiling the sphere with rational Bézier patches, Technical report UUCS-88-009, Computer Science, University of Utah (1988).
  * [32] J. Dölz, H. Harbrecht, S. Kurz, S. Schöps, F. Wolf, A fast isogeometric BEM for the three dimensional Laplace-and Helmholtz problems, Comput Methods Appl Mech Engrg 330 (2018) 83–101.