verdeRGB46,139,87

11institutetext: Antonella Falini, Dipartimento di Informatica, Università
degli Studi di Bari Aldo Moro, Italy 11email<EMAIL_ADDRESS>
22institutetext: Carlotta Giannelli, Dipartimento di Matematica e Informatica
Ulisse Dini, Università degli Studi di Firenze, Italy 22email:
<EMAIL_ADDRESS>
33institutetext: Tadej Kanduč, Faculty of Mathematics and Physics, University
of Ljubljana, Ljubljana, Slovenia 33email<EMAIL_ADDRESS>
44institutetext: M. Lucia Sampoli, Dipartimento di Ingegneria
dell’Informazione e Scienze Matematiche, Università degli Studi di Siena,
Italy44email<EMAIL_ADDRESS>
55institutetext: Alessandra Sestini, Dipartimento di Matematica e Informatica
Ulisse Dini, Università degli Studi di Firenze, Italy 55email:
<EMAIL_ADDRESS>
# A collocation IGA-BEM for 3D potential problems on unbounded domains

Antonella Falini  Carlotta Giannelli  Tadej Kanduč  M. Lucia Sampoli and
Alessandra Sestini

###### Abstract

In this paper the numerical solution of potential problems defined on 3D
unbounded domains is addressed with Boundary Element Methods (BEMs), since in
this way the problem is studied only on the boundary, and thus any finite
approximation of the infinite domain can be avoided. The isogeometric analysis
(IGA) setting is considered and in particular B-splines and NURBS functions
are taken into account. In order to exploit all the possible benefits from
using spline spaces, an important point is the development of specific
cubature formulas for weakly and nearly singular integrals. Our proposal for
this aim is based on spline quasi-interpolation and on the use of a spline
product formula. Besides that, a robust singularity extraction procedure is
introduced as a preliminary step and an efficient function-by-function
assembly phase is adopted. A selection of numerical examples confirms that the
numerical solutions reach the expected convergence orders.

##  1 Introduction

The framework of Isogeometric Analysis (IGA), introduced in 2005 by Thomas J.
R. Hughes and coauthors, has emerged as a very attractive approach to solve
differential problems LibroHughes ; Hughes_2005 . IGA bridges the gap between
Computational Mechanics and Computer Aided Design (CAD) by employing a unified
representation of geometric objects and physical quantities. In particular,
the spline basis functions commonly used for CAD geometry parameterization are
also used to span the discretization space adopted for the numerical solution.
This not only eliminates the need for slow and error-prone conversion
processes between the geometry and simulation models, but also enables the
possibility of a widespread use of design optimization tools in simulation-
based engineering.

While IGA has achieved great success within the finite element paradigm by
attracting most of the IGA researchers, a remaining critical challenge is how
to efficiently and exactly parameterize a volumetric domain from its boundary
al2016Brep ; cohen_etal2010 ; pan2020volumetric . Three dimensional problems
need trivariate NURBS solids for analysis, but CAD systems use a boundary
representation where only bivariate NURBS are involved. A possible solution to
this issue consists in adopting boundary element methods (BEMs) beer2020LIBRO
; BEMbook , whose isogeometric formulation is briefly referred to as IGA-BEM
and has recently been succesfully applied to several classical problems, e.g.,
Laplace and Helmholtz equation TauRodHug ; Dolz18 , elastostatics and
elasticity Simp2012 ; Wang15 ; beer_etal2017 ; nguyen2017 , potential flow
gong2017  and wave-resistance problems ginnis . BEMs reformulate the original
differential problem through the fundamental solution into boundary integral
equations. Thus, adopting a boundary representation, they can be easily
coupled with standard CAD techniques and ensure dimension reduction of the
computational domain where the numerical simulation is developed. Besides
these key points, they are attractive because offering an easy treatment of
problems on unbounded domains, see for example the potential or acoustic
problems considered in Coox ; Nguyen16 . Within this kind of problems, screen
or crack problems are particularly important for applications, being in this
case the domain given by the whole Euclidean space external to an open limited
arc in 2D or surface in 3D.

Relying on the general NURBS parametric description of the domain boundary and
on (refinements of) the related spline space for the approximation of the
unknowns, in this paper we confirm the effectiveness of the IGA-BEM in the 3D
setting, previously experimented in 2D ADSS1 . In particular, we focus on a
collocation discretization strategy for 3D exterior potential problems. The
combination of the superior smoothness of spline basis functions with the low
computational cost and simplicity of collocation techniques seems to
constitute an optimal basis for accurately modeling complex and
computationally demanding problems, at least when smooth geometries are dealt
with, as also outlined in TauRodHug . However it is worth mentioning that
considering a Galerkin based variant of the developed method could be of great
interest in the future, either because the method could become more robust for
treating problems with low regular solutions sutradhar2008SGBEM  or even from
a more theoretical point of view, since most of the convergence results
available in the literature refer to the Galerkin approach.
Since an effective implementation of BEMs surely requires suitable rules for
the approximation of arising singular integrals, in this paper we focus on
such aspect. In particular, the cubature rules for bivariate weakly singular
integrals developed in Nash20  are here firstly applied to the 3D IGA-BEM
context. These rules are based on spline Quasi-Interpolation (QI) and on the
usage of a spline product formula Morken91 . Specifically formulated for
weakly singular integrals containing a tensor-product B-spline factor in the
integrand, their peculiarity consists in being applicable on its whole
support. Since tensor product B-splines are here adopted to span the
discretization space, this feature implies that the assembly phase of the
system matrix can be done function-by function. This is a clear advantage of
our approach, as already outlined in the 2D setting ACDS3 . Furthermore, as a
preliminary step the multiplicative singularity extraction procedure
considered in Nash20  is here improved with the aim to apply the QI operator
adopted by our rules to sufficiently smooth functions. Note that in the IGA-
BEM setting, according to the distance between the source point and the
integration domain, not all the integrals to be dealt with are singular, some
being near singular and others regular. The final accuracy of our IGA-BEM
implementation gained from treating near singular integrals with the same
strategy adopted for singular ones. For regular integrals, a simplified
formulation of the mentioned cubature rules has been developed, relying on the
same QI operator, but possibly selecting for it a different spline degree
and/or knot spacing.

The structure of the paper is as follows. Section 2 recalls the indirect
integral formulation of a potential problem on a 3D unbounded domain, and then
it introduces the related discretization based on an isogeometric collocation
BEM. Section 3 presents cubature rules based on quasi-interpolation and on the
spline product formula. Section 4 reports the results obtained with a Matlab
implementation of the scheme for a screen problem and for two different
problems on the unbounded domain external to a torus. Finally, some
conclusions are given in Section 5.

##  2 Isogeometric indirect BEMs on 3D unbounded domains

As well as for interior problems, even for general exterior problems different
integral formulations are available in the literature, coming from different
simple or double layer representation formulas. BEMs derived from an integral
equation which directly involves the unknown Cauchy datum are called (simple
or double layer) “direct” methods. As an alternative, (simple or double layer)
“indirect” approaches can also be adopted. In such case the unknown appearing
in the integral equation and in the associated representation formula is a
different function, usually referred to as “density” function. In this paper
we adopt a BEM indirect formulation, since it is the only one applicable also
to screen/crack problems, costabel1986principles . In particular we focus on
such formulation for 3D Dirichlet exterior potential problems,

$\left\\{\begin{array}[]{l}\Delta u=0\,,\,\,{\bf x}\in{\mathop{{\rm
I}\kern-1.99997pt{\rm R}}\nolimits}^{3}\setminus{\overline{\Omega}}\,,\cr
u=g\,,\,\,{\bf x}\in\Gamma\,,\cr\left.\begin{array}[]{l}|u({\bf
x})|\rightarrow 0\cr\|\nabla u({\bf x})\|=O(\|{\bf
x}\|^{-2})\end{array}\right\\}\,\mbox{for }\|{\bf
x}\|\rightarrow\,+\infty\,,\par\end{array}\right.$ (1)

where $g$ is the given Dirichlet datum and $\Gamma:=\partial\Omega,$ with
$\Omega$ denoting a bounded domain in ${\mathop{{\rm I}\kern-1.99997pt{\rm
R}}\nolimits}^{3}$. Note that for screen problems $\Omega$ has to be replaced
with $\Gamma,$ where $\Gamma$ is a given open and limited surface. Observe
also that, since we deal with unbounded domains, it is necessary to specify
the required behavior of the harmonic solution at infinity to ensure
uniqueness BEMbook . The boundary integral formulation of (1) is given by the
following BIE,

$V\psi({\bf x}):=\int_{\Gamma}U({\bf x},{\bf y})\psi({\bf
y})\,\mathrm{d}\Gamma_{{\bf y}}\,=\,g({\bf x})\,,\qquad{\bf x}\in\Gamma\,,$
(2)

where the unknown density function $\psi$ represents the jump across the
boundary of the normal derivative of $u$ (in the outward direction from
${\mathop{{\rm I}\kern-1.99997pt{\rm
R}}\nolimits}^{3}\setminus\overline{\Omega}$ to $\Omega$) and the kernel $U$
is the fundamental solution of the Laplace equation in 3D,

$U({\bf x}\,,\,{\bf y}):=\frac{1}{4\pi}\,\frac{1}{\|{\bf x}-{\bf y}\|_{2}}.$
(3)

The operator $V$ defined in (2) is an elliptic isomorphism between its
functional domain $H^{-\frac{1}{2}}(\Gamma)$ and its codomain
$H^{\frac{1}{2}}(\Gamma),$ while the BIE also defined in (2) is a Fredholm
integral equation of the first kind with weakly singular kernel.

Once the density function $\psi$ has been (approximately) computed, the
(approximated) expression of $u$ at any point in ${\mathop{{\rm
I}\kern-1.99997pt{\rm R}}\nolimits}^{3}\setminus\overline{\Omega}$ can be
derived by relying on its simple layer representation formula,

$u({\bf x})=\int_{\Gamma}U({\bf x}\ ,{\bf y})\,\psi({\bf
y})\,\mathrm{d}\Gamma_{{\bf y}}\,,\qquad{\bf x}\in{\mathop{{\rm
I}\kern-1.99997pt{\rm R}}\nolimits}^{3}\setminus\overline{\Omega}\,.$ (4)

Note that this formula automatically ensures the required decay behavior of
$u$ at infinity if $\psi$ is continuous in $\Gamma$ BEMbook , and also that
its numerical approximation can be challenging only when ${\bf x}$ is very
close to $\Gamma.$

To compute a numerical solution of the BIE in (2) we adopt collocation, since
this is probably the easiest discretization method suited for smooth problems
(i.e., problems with a smooth $\Gamma$ and sufficiently smooth Dirichlet
datum). However, we observe that the related theory is not yet fully developed
in particular for the 3D case; for example in Sections 9.2 and 9.3 of
Atkinson_2009  the available theoretical results can be applied just for the
direct approach which is not usable for screen problems, as we already
remarked. In general, after having numerically solved the considered BIE with
a BEM approach, the representation formula associated with it has to be
considered in order to compute values of the approximated solution inside the
domain. However there are several applications, e.g. the simulation for screen
or crack problems, which do not require such post-processing step, since their
physically relevant quantity is the unknown Cauchy datum on the boundary
itself, see for example costabel1986principles .

We assume that the geometry is CAD-generated, that is $\overline{\Gamma}$ has
a parametric representation, $\overline{\Gamma}=\mbox{Image}({\bf F})\,,$
where ${\bf F}:R\rightarrow{\mathop{{\rm I}\kern-1.99997pt{\rm
R}}\nolimits}^{3}$ is a one-to-one mapping111With the exception of the
boundary of the parametric domain in case of a closed $\Gamma$ or with
cylinder-like shape. defined on the parametric rectangular domain
$R:=[a_{1}\,,\,b_{1}]\times[a_{2}\,,\,b_{2}]$ in the following standard NURBS
form,

${\bf F}({\bf
t}):=\frac{\displaystyle{\sum_{i=0}^{n_{1}}\sum_{j=0}^{n_{2}}}w_{i,j}\,{\bf
Q}_{i,j}\,B^{d_{1}}_{i,T_{1}}(t_{1})B^{d_{2}}_{j,T_{2}}(t_{2})}{\displaystyle{\sum_{i=0}^{n_{1}}\sum_{j=0}^{n_{2}}}w_{i,j}\,B^{d_{1}}_{i,T_{1}}(t_{1})B^{d_{2}}_{j,T_{2}}(t_{2})}\,,\quad{\bf
t}\in R\,,$ (5)

where ${\bf t}:=(t_{1},t_{2})$. Each set of functions
$B^{d_{j}}_{i,T_{j}},i=0,\ldots,n_{j},j=1,2$ is made up of univariate B-spline
functions defined in $[a_{j}\,,\,b_{j}],$ of degree $d_{j}$ and with extended
knot vector $T_{j}$. We assume that the mapping ${\bf F}$ is at least
differentiable. The set $\\{{\bf Q}_{i,j}\in{\mathop{{\rm
I}\kern-1.99997pt{\rm R}}\nolimits}^{3},i=0,\ldots,n_{1},j=0,\ldots,n_{2}\\}$
defines the net of control points222In order to avoid degenerate nets we
assume the control net composed by distinct points. and each $w_{i,j}$ is a
positive weight used to obtain additional shape control (obviously, if the
weights are all equal to a common constant value, ${\bf F}$ becomes just a
vector tensor-product polynomial spline represented in B-spline form). This
surface representation is often adopted in CAD environments to design free-
form surfaces and it has the facility of involving just one patch, that is a
unique continuous vector function ${\bf F}$ is used to define the whole
considered surface. On this concern we note that clearly this representation
can be not flexible enough for complex geometries characterized for example by
linkage curves with just geometric continuity. Note also that, even if the
NURBS form was originally introduced for guaranteeing exact representation of
conic sections and quadric patches, it cannot be used for example for
representing closed surfaces with zero genus without singularities in the
parameterization. For example in reference TauRodHug  the one-patch NURBS
parameterization of the sphere is singular at the poles and this has required
the development of special cubature formulas to obtain accurate results. In
order to deal with surface representations more flexible and more suitable for
the analysis, multi-patch NURBS representation could be considered, as well as
different kinds of spline spaces, e.g., multi-degree polar splines TSHH20 .

Focusing on the basic standard single patch NURBS representation and referring
for brevity only to the case of clamped knot vectors $T_{j},j=1,2$, let us
define ${\cal T}_{j}$ as a possible refinement of $T_{j},j=1,2,$ obtained by
inserting new breakpoints in $[a_{j}\,,\,b_{j}],$ or by increasing the
multiplicity of a knot already included in $T_{j}.$ Denoting with $h_{j}$ the
maximal distance between successive knots in ${\cal T}_{j},j=1,2$ and with
$S_{h_{j}}$ the univariate spline space of dimension $m_{j}+1\,,m_{j}\geq
n_{j},$ spanned by the B-splines of degree $d_{j}$ associated with ${\cal
T}_{j},$ we can define the tensor-product spline space ${\cal
S}^{h}:=S_{h_{1}}\otimes S_{h_{2}}$ of dimension $D:=(m_{1}+1)(m_{2}+1),$
where $h:=\max\\{h_{1},\,h_{2}\\}.$ Note that in the following we simplify the
notation by denoting as $B_{j}$ the bivariate B-spline $B^{d_{1}}_{i_{1},{\cal
T}_{1}}B^{d_{2}}_{i_{2},{\cal T}_{2}}$ generating ${\cal S}^{h},$ where
$j=i_{2}(m_{1}+1)+i_{1}.$ Thus, the density function $\psi$ is approximated in
a finite dimensional space $\hat{\cal S}^{h}$ derived from ${\cal S}^{h}$ by
lifting all the B-splines to the physical space

$\hat{\cal S}^{h}\,:=\,\mathrm{span}\\{B_{j}\circ{\bf
F}^{-1},j=0,\ldots,D-1\\}\,.$

Using collocation, the coefficient vector
$\mbox{\boldmath$\alpha$}_{h}=(\alpha_{0},\cdots,\alpha_{D-1})^{T}\in{\mathop{{\rm
I}\kern-1.99997pt{\rm R}}\nolimits}^{D}$ associated with the IGA-BEM discrete
solution $\psi_{h}\in\hat{\cal S}^{h},\ \psi_{h}({\bf
x})=\sum_{j=0}^{D-1}\alpha_{j}(B_{j}\circ{\bf F}^{-1})({\bf x})\,,$ is
determined by solving the following linear system with $D$ equations,

$A_{h}\mbox{\boldmath$\alpha$}_{h}={\bf g}_{h}\,,$ (6)

where

$(A_{h})_{i,j}:=\int_{\Gamma}U({\bf x}_{i},{\bf y})\ (B_{j}\circ{\bf
F}^{-1})({\bf y})\ d\Gamma_{{\bf y}}\,,$ (7)

and $({\bf g}_{h})_{i}:=g({\bf x}_{i})\,.$

The points ${\bf x}_{i},i=0,\ldots,D-1,$ are $D$ distinct collocation points
belonging to $\Gamma$ defined as the image through ${\bf F}$ of the Cartesian
product of two sets of distinct abscissas, defined in $[a_{1},\,b_{1}]$ and in
$[a_{2},\,b_{2}],$ respectively, as for example Greville abscissas (or
suitable variants, see for example Wang15 , to avoid collocation on the
boundary of $\Gamma$). Thus, for each collocation point ${\bf
x}_{i}\in\Gamma,i=1,\ldots,D,$ there exists ${\bf s}_{i}\in R$ such that ${\bf
x}_{i}={\bf F}({\bf s}_{i}),$ with ${\bf s}_{i}\neq{\bf s}_{j},i\neq j.$

Then, moving to the parametric domain, the entries of the matrix $A_{h}$ can
be written as follows,

$(A_{h})_{i,j}\,=\,\int_{R_{j}}U({\bf F}({\bf s}_{i}),{\bf F}({\bf t}))\
B_{j}({\bf t})\,J({\bf t})\ dt_{1}dt_{2}\,,$ (8)

where $R_{j}:=\mathrm{supp}(B_{j})$ and where $J$ represents the infinitesimal
surface area element,

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

Regarding the matrix $A_{h},$ we observe that it is a full matrix whose
characteristics also depend on the distribution of the collocation points in
$\Gamma$. On this concern we mention that the nonsingularity of the
collocation matrix descending from a second kind Fredholm BIE has been
theoretically studied in Atkinson_2009 . However our $A_{h}$ comes from a
Fredholm integral equation of the first kind and, up to our knowledge,
sufficient conditions ensuring such property are not available in this case,
unless for very special situations Costabel92screen3D . We can only add that
in all our experiments we have always dealt with positive definite matrices,
see also the comment about the condition number in the numerical experiments
presented in Section 4.

In order to be able to compute the discrete solution efficiently and
accurately, it is important to have an efficient and stable method for solving
the linear system in (6) but even more to consider suitable cubature rules for
weakly singular integrals to be used in the assembly phase. In this paper we
focus on this second issue, extending to cubature the quadrature rules
previously developed for analogous univariate weakly singular integrals. As
already mentioned in the introduction, the use of such rules (which are
specific for integrals containing an explicit B-spline factor in the
integrand) allows us to avoid the element-by-element assembly strategy, since
they are directly applied on the whole support of each $B_{j}.$

##  3 Cubature rules

In this section we present the cubature rules introduced in Nash20  and here
firstly applied in the 3D IGA-BEM context for the numerical approximation of
all the integrals on the right hand side of (8). Such rules rely on spline
quasi-interpolation and are an extension to the bivariate case of those
firstly studied in the univariate setting CFSS18  and successfully tested also
working in adaptive spline spaces FGKSS_2019 .

We first observe that the integrals defining the entries of $A_{h}$ are weakly
singular only if the source point ${\bf s}_{i}$ belongs to
${\overline{R}}_{j}.$ If this is not the case, they are nearly singular or
regular, depending on the distance between ${\bf s}_{i}$ and the integration
domain $R_{j},$ clearly taking also into account if $\Gamma$ is a closed
surface. This initial classification is important because nearly singular
integrals need as much care as weakly singular ones. In our implementation we
adopt for them the same strategy used for weakly singular integrals, while a
simplified version of the developed rules is used to approximate regular
integrals, possibly using also higher QI spline degree in such a case.

Referring first of all to singular and nearly singular integrals, we observe
that, as well as in the 2D case ACDS3 ; FGKSS_2019 , a preliminary singularity
extraction procedure is necessary in order to increase the performance of our
cubature rule. In Nash20  an analysis has been developed to show that two
singularity extraction procedures, based on a Taylor expansion of ${\bf
F}({\bf t})$ about ${\bf s}$, are possible, one using an additive splitting of
$U$ and the other of multiplicative type. Considering that the first splitting
requires the numerical computation of two integrals instead of one, here we
have chosen to rely on the multiplicative technique. In this case the idea is
to rewrite the entries of the matrix $A_{h}$ given in (8) in the following
equivalent form,

$(A_{h})_{i,j}:=\int_{R_{j}}K_{\bf F}({\bf s}_{i},{\bf t})\ \rho_{{\bf
s}_{i},{\bf F}}({\bf t})B_{j}({\bf t})\,J({\bf t})\
\mathrm{d}t_{1}\mathrm{d}t_{2}\,,$ (9)

where

$\rho_{{\bf s},{\bf F}}({\bf t}):=\frac{U({\bf F}({\bf s}),{\bf F}({\bf
t}))}{K_{\bf F}({\bf s},{\bf t})}\,,$

and $K_{\bf F}$ is a suitably locally defined weakly singular kernel such that
the function $\rho_{{\bf s},{\bf F}}$ is at least continuous. In particular in
Nash20  we considered the following definition of such reference kernel,

$K_{\bf F}:=K_{\bf F}({\bf s},{\bf t}):=\frac{1}{4\pi}\ \frac{1}{\sqrt{P_{{\bf
s},{\bf F}}({\bf t})}}\,,$ (10)

with $P_{{\bf s},{\bf F}}$ denoting the bivariate homogeneous quadratic
polynomial vanishing with its first derivatives at ${\bf t}={\bf s}$ and
obtained by replacing ${\bf F}({\bf t})$ with its first degree Taylor
expansion about ${\bf s}$ in $\|{\bf F}({\bf t})-{\bf F}({\bf s})\|_{2}^{2}$.
Note anyway that for a general surface this definition just ensures continuity
to $\rho_{{\bf s},{\bf F}}$ at ${\bf t}={\bf s}$. However, if ${\bf F}$ is
sufficiently smooth at the singular point ${\bf s}$, it is possible to use
higher order expansions for ${\bf F}$, in order to ensure higher regularity to
$\rho_{{\bf s},{\bf F}}$. In particular here we consider an expansion with up
to $3$ terms, which ensures up to $C^{2}$ smoothness to $\rho_{{\bf s},{\bf
F}}$. Details about this improved singularity extraction procedure will be
presented in a forthcoming paper333T. Kanduč. Isoparametric singularity
extraction technique for 3D potential problems in BEM, arXiv:2203.11538.,
since they form a technical part necessary for an exhaustive explanation but
too long for being here reported. We just mention that in this case $K_{\bf
F}$ is defined with the following expression which reduces to (10) for $n=1,$

$K_{\bf F}({\bf s},{\bf t}):=\frac{1}{4\pi}\ \sum_{\ell=1}^{n}\sqrt{P_{{\bf
s},{\bf F}}({\bf t})}^{-2\ell+1}P_{3\ell-3}^{({\bf s},{\bf F})}({\bf t}),$
(11)

where $P_{3\ell-3}^{({\bf s},{\bf F})}$ are appropriate homogeneous
polynomials of total degree $3\ell-3\,,\ell=1,\ldots,n.$ Concerning this
improved singularity extraction procedure, finally we add that, when $n>1$ and
the computational mesh is not fine enough, a correction term $\eta\|{\bf
s}-{\bf t}\|_{2}^{2}$ for some suitably small $\eta>0$ has to be added on the
right hand side of (11) to avoid the introduction of new singularities in
$\rho_{{\bf s},{\bf F}}$.

After this preliminary rewriting of the entries of the matrix $A_{h},$ each of
them is approximated as follows,

$(A_{h})_{i,j}\approx\int_{R_{j}}K_{\bf F}({\bf s}_{i},{\bf t})\ \left(Q_{{\bf
p},\Delta_{j}}\left(J\rho_{{\bf s}_{i},{\bf F}}\right)\right)({\bf t})\
B_{j}({\bf t})\ \mathrm{d}t_{1}\mathrm{d}t_{2}\,,$ (12)

where $\Delta_{j}$ is a uniform partition of $R_{j}$ and $Q_{{\bf
p},\Delta_{j}}$ denotes a tensor-product QI operator which approximates a
bivariate function with a tensor-product spline of bi-degree ${\bf
p}=(p_{1},p_{2})$ with respect to the partition $\Delta_{j}.$

In order to compute the integral on the right hand side of (12), the spline
$Q_{{\bf p},\Delta_{j}}\left(J\ \rho_{{\bf s}_{i},{\bf F}}\right)$ is
multiplied by the B-spline $B_{j},$ using for this aim the extension to the
tensor-product setting of the Mørken formula Morken91  for spline product. As
a result, the whole product $\left(J\ B_{j}\ \rho_{{\bf s}_{i},{\bf
F}}\right)$ is approximated with a local tensor-product spline function
defined in $R_{j},$ having bi-degree ${\bf p}+{\bf d},$ and breakpoints in
$R_{j}.$ We observe that, in order to implement the so derived cubature rules,
the evaluation of the following modified moments is firstly necessary,

$\mu_{K_{\bf F},B}({\bf s}):=\int_{\mathrm{supp}(B)}K_{\bf F}({\bf s},{\bf
t})\ B({\bf t})\ \mathrm{d}t_{1}\mathrm{d}t_{2}\,,$ (13)

with $B$ denoting any tensor-product B-spline. This computation is done by
relying on the recursion formula for B-splines and on the preliminary analytic
computation of element integrals whose integrand is the product of a monomial
times $K_{\bf F},$ procedure possible for a kernel $K_{\bf F}$ defined as in
(10) or more generally as in (11).

When the integral defining the $(i,j)$–th entry of $A_{h}$ is classified as a
regular one, we do not need the preliminary computation of the modified
moments and the strategy is simplified by directly setting

$(A_{h})_{i,j}\approx\int_{R_{j}}\left(Q_{{\bf q},\Lambda_{j}}\left(U({\bf
F}({\bf s}_{i}),\cdot)J\right)\right)({\bf t})\ B_{j}({\bf t})\
\mathrm{d}t_{1}\mathrm{d}t_{2}\,,$ (14)

where ${\bf q}=(q_{1}\,,\,q_{2})$ denotes a bi-degree for the local splines
possibly different from ${\bf p}$ and $\Lambda_{j}$ is a uniform partition of
$R_{j}$, also possibly different from $\Delta_{j}.$

Concerning the QI spline operator used by the proposed cubature rules, for our
experiments we have always considered the uniform tensor-product formulation
of the derivative free variant of the Hermite QI operator introduced in
MSbit09 , see for example BGMS16  for an introduction to its formulation in
the bivariate setting and for an application in the context of adaptive
bivariate spline approximation. This univariate spline QI operator has
approximation order $p+1$ and in the context of numerical integration it has
special interest because its application to both regular and weakly singular
quadrature is superconvergent when even spline degrees are considered on
uniform knot distributions MSJcam12 ; CFSS18 . Thus for smooth functions
$Q_{{\bf p},\Delta_{j}}$ has approximation order $\min\\{p_{1},p_{2}\\}+1$ and
it inherits the superconvergence feature when analogously applied for
cubature. Observe however that, in our approach, choosing values of
$p_{k},k=1,2,$ greater than $2$ is not useful, because of the limited
regularity of the function $\rho_{{\bf s}_{i},{\bf F}}.$ It can be profitable
instead to select $q_{k}>2,k=1,2$, if ${\bf F}$ is smooth in $R_{j}.$
Considering the assumed analytic expression of ${\bf F},$ this is surely true
when the interior of $R_{j}$ does not contain any straight line
$t_{j}=\tau_{j},$ with $\tau_{j}$ geometric knot in the $j$–th parametric
direction, that is $\tau_{j}\in T_{j}\,,j=1,2,$ where $T_{j},j=1,2$ are the
extended knot vectors used to define ${\bf F}.$

##  4 Numerical simulations

In this section we test the uniform formulation of the proposed IGA boundary
element model for the numerical solution of two Laplace problems on two
different 3D unbounded domains. Furthermore a nonuniform implementation of the
scheme is also tested for solving a more difficult Laplace problem exterior to
the first considered domain. In all the examples the improved singularity
extraction procedure is used within the developed singular cubature rules,
always setting the parameter $\eta$ (see Section 3) equal to $1.$ Concerning
the discretization, in all the presented experiments approximate densities are
constructed on a given initial mesh and also on successive finer ones obtained
from a dyadic $h$-refinement procedure which inserts a simple knot at the
midpoint between every pair of consecutive breakpoints in each parametric
direction.
For the first two considered examples the density $\psi$ is known, while the
Dirichlet datum $g$ is not available in closed form. Thus, the entries of the
right-hand side in (6) are computed numerically,

$({\bf g}_{h})_{i}=\int_{R}U({\bf F}({\bf s}_{i}),{\bf F}({\bf t}))\
(\psi\circ{\bf F})({\bf t})\,J({\bf t})\ \mathrm{d}t_{1}\mathrm{d}t_{2}\,,$
(15)

by applying the same approach adopted to compute the system matrix. For these
examples the accuracy of the scheme is measured by the $L^{2}(\Gamma)$ norm of
the absolute density error $(\psi_{h}-\psi),$ in particular considering its
dependency on the mesh size $h.$ The third exterior Laplace problem is based
on Test $5.3$ of Dolz18  where the more general Helmholtz equation is solved.
Since in this case the analytic expression of the exact potential $u$ is
known, the Dirichlet datum is given and the accuracy in the potential
approximation $u_{h}\approx u$ can be measured, provided that the
representation formula given in (4) is previously applied to $\psi_{h}$ to
define $u_{h}$ in the considered unbounded domain ${\mathop{{\rm
I}\kern-1.99997pt{\rm R}}\nolimits}^{3}\setminus\overline{\Omega}.$

.

###  4.1 Problem on an unbounded domain exterior to a torus

We consider a problem exterior to a toroidal surface $\Gamma$ with a major
radius $\rho_{M}=3$ and a minor radius $\rho_{m}=1$. The geometry is
represented in NURBS form by defining as in (5) the vector function ${\bf F}$
on the parametric domain $R=[-3\,,\,3]\times[-1\,,\,1]$ with ${\bf d}=(2,2),$
fixing the following knot vectors,

$\displaystyle
T_{1}=(-3,\;-3,\;-3,\;-1.5,\;-1.5,\;-1.5,\;0,\;0,\;0,\;1.5,\;1.5,\;1.5,\;3,\;3,\;3),$
$\displaystyle
T_{2}=(-1,\;-1,\;-1,\;-0.5,\;-0.5,\;-0.5,\;0,\;0,\;0,\;0.5,\;0.5,\;0.5,\;1,\;1,\;1),$

see solid lines in Fig. 1(a). In this form the torus is composed by the $16$
rational biquadratic sub-patches shown in Fig. 1(c) (separated by solid
lines), which are locally represented in rational Bernstein form, since all
the breakpoints have maximal multiplicity (equal to $3$) in both the knot
vectors $T_{j},j=1,2.$ The associated control net has the size $12\times 12;$
for clarity in Fig. 1(b) we show a sub-patch and the related $3\times 3$ sub-
net of control points with the corresponding weights. Dashed lines in Fig.
1(a) and Fig. 1(c) represent additional inserted knots for the initial
computational mesh. Note that, despite the maximal multiplicity of the
breakpoints in $T_{j},j=1,2,$ the vector function ${\bf F}$ is $C^{1}$
continuous because of the symmetry of the considered control net and of the
related weights. We remark that, even if not strictly necessary, maximal
multiplicity is used for the geometric knots to ensure to ${\bf F}$ maximal
regularity when restricted to the support of a trial function $B_{j}.$ As
explained at the end of the previous section, this is useful for the accuracy
of our cubature rules when regular integrals are dealt with, since the support
of a trial function can never cross a line $t_{j}=\tau_{j}\in T_{j},j=1,2\,,$
i.e. solid lines in Fig. 1(a).
The exact function $\psi\circ{\bf F}$ is set to

$\displaystyle(\psi\circ{\bf F})({\bf t})=\sin(4/3\;t_{1}\,t_{2})$

and its shape in the parametric domain is depicted in Fig. 1(d). Each rational
sub-patch of the torus contains the same number of collocation points obtained
from the Cartesian product of the improved Greville abscissas Wang15 , defined
on the corresponding sub-rectangles of $R$, and mapped to $\Gamma$ by ${\bf
F}.$
Two different experiments have been performed for this example. The first one
is a standard test where initially ${\cal T}_{2}=T_{2}$ and ${\cal T}_{1}$ is
obtained by merging $T_{1}$ with
$\\{-2.5\,,-2\,,-1\,,-0.5\,,0.5\,,1\,,2\,,2.5\\}.$ Since in all the
simulations the geometric knots keep maximal multiplicity, at any refinement
level also discontinuous functions belong to the space ${\cal S}_{h}.$ In the
second experiment initially we define ${\cal T}_{2}$ still from $T_{2}$ but by
decreasing the multiplicity of each inner knot by one. ${\cal T}_{1}$ is
initialised as in the other case but again by reducing by one the multiplicity
of all the geometric knots. Thus $3\times 3$ degrees of freedom are eliminated
with respect to the previous experiment and as a consequence the simulation
space is always included in $C^{0}(R).$ Note that anyway we keep unchanged the
set of collocation points so we end up with an overdetermined linear system
which is solved in the least squares sense.

In both experiments the results shown are obtained by deriving the reference
kernel $K_{\bf F}$ by using the first two terms of the Taylor expansion of
${\bf F}.$ For the weakly or near singular cubature the QI bi-degree ${\bf p}$
is set to ${\bf p}=(2,2),$ using in each parametric direction 13 uniform nodes
on the support of any B-spline trial function. Due to higher global smoothness
of integrands, for regular integrals we can profitably set ${\bf q}=(4,4)$ and
just use 7 uniform nodes on the support of any B-spline in each parametric
direction.

Convergence plots of the behaviour of approximation error
$\|\psi-\psi_{h}\|_{L^{2}(\Gamma)}$ using uniform $h$-refinement strategy
together with the theoretical convergence order 3 is shown in Fig. 3(a).
Looking at the figure, we see that both the plots have the expected behaviour
and also that the $C^{0}$ one is very near to the other. The system matrix $A$
is well conditioned for all tests (in the last step its condition number for
the first experiment is about $5.1\cdot 10^{3}$).

(a) Initial mesh in the parametric domain $R$ \begin{overpic}[trim=28.45274pt
14.22636pt 14.22636pt 56.9055pt,clip={true},height=99.58464pt]{torus_net.eps}
\put(12.0,41.0){\scriptsize 1} \put(6.0,77.0){\scriptsize$\sqrt{2}/2$}
\put(30.0,80.0){\scriptsize 1} \par\put(65.0,18.0){\scriptsize$\sqrt{2}/2$}
\put(75.0,48.0){\scriptsize$1/2$} \put(72.0,67.0){\scriptsize$\sqrt{2}/2$}
\par\put(140.0,37.0){\scriptsize$1$} \put(127.0,61.0){\scriptsize$\sqrt{2}/2$}
\put(123.0,76.0){\scriptsize$1$} \end{overpic} (b) Control net and weights for
one patch (c) Surface $\Gamma=\partial\Omega={\bf F}(R)$ and initial mesh (d)
Graph of composition with ${\bf F}$ of the exact BIE solution $\psi$

Figure 1: Example 4.1: problem exterior to a torus.

###  4.2 A screen problem for a curved obstacle

As a second example we consider a screen problem (1) with $\Gamma$ taken as
the limited portion of a saddle open surface parametrized by ${\bf
F}:R\rightarrow\overline{\Gamma}$, with $R=[-1\,,\,1]^{2}$ and ${\bf
F}(t_{1},t_{2})=(t_{1}\,,\,t_{2}\,,\,t_{1}^{2}-t_{2}^{2})$. The Dirichlet
datum $g$ is numerically computed such that the composition with ${\bf F}$ of
the exact density function, solution of Eq. (2), is given by

$(\psi\circ{\bf F})({\bf t})=\sqrt{1+4t_{1}^{2}-0.5t_{2}^{2}}.$

In Fig. 2 we show the geometry and $\psi\circ{\bf F}.$ Note that the
considered ${\bf F}$ can be written as in (5) with ${\bf d}=(2,2)$ and
$T_{j}=(-1,\,-1,\,-1,\,1,\,1,\,1)\,,\,j=1,2\,,$ since polynomials (restricted
to $R$) are a subset of NURBS. Thus the problem is discretised by using
quadratic B-splines and setting ${\bf p}=(2,2)$ for the QI approximation of
singular and nearly singular integrals, with $7$ quadrature uniform nodes on
the support of any B-spline trial function in each parametric direction. The
reference kernel $K_{\bf F}$ is obtained in this case using the first three
terms in the Taylor expansion of ${\bf F},$ while the improved Greville
abscissas Wang15  are employed as collocation points also in this case.

In Fig. 3(b) we show the convergence plots for ${\bf q}$ varying from $(2,2)$
to $(4,4)$, using from $5$ to $7$ uniform quadrature nodes per B-spline
support in each parametric direction for regular integrals. Better accuracy of
the numerical solution corresponds to higher ${\bf q}$ but the order $3$ is
observed in all the cases. Finally, the condition number for the system matrix
$A$ has order $10^{2}$ at the last level.

(a) Surface $\Gamma=\partial\Omega={\bf F}(R)$ (b) Graph of composition with
${\bf F}$ of the exact BIE solution $\psi$

Figure 2: Example 4.2: screen problem.

(a) Problem exterior to a torus $\big{(}{\bf d}=(2,2)\big{)}$ (b) Screen
problem $\big{(}{\bf d}=(2,2)\big{)}$

Figure 3: Convergence plots in $L^{2}(\Gamma)$ norm of the density error
$\psi-\psi_{h}$ for Examples 4.1 and 4.2.

###  4.3 A nearly singular potential problem exterior to a torus

We consider the same toroidal surface presented in Example 4.1, solving the
problem (1) with the Dirichlet datum $g$ chosen as,

$\displaystyle g({\bf x}):=\frac{1}{\|{\bf
x}-\mbox{\boldmath$\nu$}\|_{2}}\quad\mbox{for }\,\,{\bf x}\in\Gamma\,,$ (16)

which is also the known analytic expression of the exact potential $u$ when
${\bf x}\in{\mathop{{\rm I}\kern-1.99997pt{\rm
R}}\nolimits}^{3}\setminus{\overline{\Omega}}.$ In particular in the reported
experiments we fix $\mbox{\boldmath$\nu$}=(-3,0,0)\in\Omega\,,$ see Fig. 4(a)
for the graph of $g\circ{\bf F}$ with this setting. This type of Dirichlet
datum gives rise to sharp gradients of the corresponding function
$\psi\circ{\bf F}$ near two opposite boundary edges of the parametric domain,
see Fig. 4(b) where we can only show the graph of $\psi_{h}\circ{\bf F},$
since the exact expression of $\psi$ is unknown.

A suitable mesh, which is nonuniform in the horizontal parametric direction,
has been manually constructed in order to get better accuracy. Clearly it
would be advisable to rely on a fully automatic adaptive mesh generation
strategy but this is out of the scope of the presented paper. The initial mesh
used for the reported experiments is shown in Fig. 4(c); as in Example 4.1,
solid lines correspond to the knots of the vectors $T_{1}$ and $T_{2}$ which
are defined as before, while dashed lines correspond to additional inserted
knots. As done in the second subtest described in Subsection 4.1, the
multiplicity of each inner geometry knot is decreased by one in order to deal
with continuous spline spaces. The collocation points are obtained from the
improved Greville abscissas and their number matches the global number of
degrees of freedom characterising the discretization space.

(a) Composition with ${\bf F}$ of $g$ (b) Composition with ${\bf F}$ of
$\psi_{h}$ (mesh from three refinements of the initial one) (c) Initial mesh
in the parameter domain (d) Distribution on $S^{2}$ of the potential error
$|u-u_{h}|$ (mesh from three refinements of the initial one)

Figure 4: Example 4.3: another problem exterior to the torus considered in
Example 4.1.

Using the representation formula (4), for this test we also evaluate the
approximated potential at some point ${\bf x}$ of the unbounded domain,

$\displaystyle u_{h}({\bf x}):=V\psi_{h}({\bf x})\quad\mbox{for }{\bf
x}\in{\mathop{{\rm I}\kern-1.99997pt{\rm
R}}\nolimits}^{3}\setminus{\overline{\Omega}}\,.$ (17)

In particular, in order to evaluate a posteriori the accuracy of our scheme,
we compute the potential error $u-u_{h}$ at $2^{m+2}$ points uniformly sampled
on a spherical surface $S^{2}$ with center located at the origin and with
radius equal to $10$ and so exterior to the considered toroidal surface. Note
that here $m$ denotes the dyadical refinement step, where $m=1$ for the
initial mesh and in particular $m=4$ was adopted to get the error distribution
in $S^{2}$ shown in Fig. 4(d). Looking at the figure, it can be observed that
the potential error has a nice constant order of magnitude on the sphere equal
to $10^{-8}$ and that its distribution is symmetric with respect to the $x$
axis where the point $\nu$ is located (it is not shown in the figure since it
is inside the sphere).

The accuracy achieved at different refinement levels $m$ is reported in Table
1 together with the corresponding global (Ndof) and directional (Ndof1, Ndof2)
number of degrees of freedom. It is measured in both the maximum and the
$L^{2}(S^{2})$ norms, both approximated using the available discrete
information on the sphere. The test has been performed by using quadratic
discretization splines and selecting ${\bf p}={\bf q}=(2,2)$. Concerning the
cubature rules, $7$ and $13$ uniform nodes on the support of every basis
function are used in each parametric direction, respectively for the
evaluation of singular/nearly singular and regular integrals. It is worth
noticing that, since the evaluation of $u_{h}$ is done at points located far
enough from the surface $\Gamma$, the integral in (17) is non singular. Hence,
the same cubature rule constructed for the evaluation of the regular integrals
has been employed for the representation formula. Finally for completeness we
mention that also for this example the coefficient matrix in (6) is well
conditioned in all the developed numerical experiments, having a condition
number about equal to $3.3\cdot 10^{3}$ at the last considered level.

$m$ | Ndof1 | Ndof2 | Ndof | $\|u-u_{h}\|_{L^{2}(S^{2})}$ | $\|u-u_{h}\|_{\infty}$
---|---|---|---|---|---
$1$ | $11$ | $9$ | $99$ | $1.24\cdot 10^{-2}$ | $5.88\cdot 10^{-4}$
$2$ | $17$ | $13$ | $221$ | $3.58\cdot 10^{-3}$ | $2.21\cdot 10^{-4}$
$3$ | $29$ | $21$ | $609$ | $1.51\cdot 10^{-4}$ | $1.00\cdot 10^{-5}$
$4$ | $53$ | $37$ | $1961$ | $4.68\cdot 10^{-7}$ | $1.82\cdot 10^{-8}$
Table 1: Example 4.3: degrees of freedom and norms of the potential error
$u-u_{h}$ obtained at different refinement level $m$ setting ${\bf d}={\bf
p}={\bf q}=(2,2).$

##  5 Conclusion

In this paper a collocation IGA-BEM scheme relying on the indirect integral
formulation of Laplace problems is proposed and tested for 3D exterior
problems, obtaining the expected convergence order in the density
approximation. Cubature rules based on spline quasi-interpolation and on a
spline product formula are applied for the numerical approximation of both
weakly and nearly singular bivariate integrals appearing in the 3D IGA-BEM
scheme. They are combined with a robust multiplicative singularity extraction
procedure to increase the accuracy of the rules. A variant of such rules is
proposed for the regular integrals to be also dealt with in the assembly
phase. This can be efficiently implemented in a function-by-function form
thanks to the peculiarity of the considered cubature.

###### Acknowledgements.

A. Falini, C. Giannelli, M. L. Sampoli, and A. Sestini are members of the
INdAM Research group GNCS. The INdAM support through GNCS and Finanziamenti
Premiali SUNRISE is gratefully acknowledged. A. Falini was also supported by
the INdAM-GNCS project “Finanziamenti Giovani Ricercatori 2019-2020”.

## References

  * [1] A. Aimi, F. Calabrò, M. Diligenti, M. L. Sampoli, G. Sangalli, and A. Sestini.  Efficient assembly based on B-spline tailored quadrature rules for the IgA-SGBEM.  Comput. Methods Appl. Mech. Engrg., 331:327–342, 2018.
  * [2] A. Aimi, M. Diligenti, M. L. Sampoli, and A. Sestini.  Isogeometric Analysis and Symmetric Galerkin BEM: a 2D numerical study.  Appl. Math. Comp., 272:173–186, 2016.
  * [3] H. Al Akhras, T. Elguedj, A. Gravouil, and M. Rochette.  Isogeometric analysis-suitable trivariate nurbs models from standard b-rep models.  Comput. Methods Appl. Mech. Engrg., 307:256–274, 2016.
  * [4] K. E. Atkinson.  The Numerical Solution of Integral Equations of the Second Kind.  Cambridge University Press, 2009.
  * [5] G. Beer, V. Mallardo, E. Ruocco, B. Marussig, J. Zechner, C. Dünser, and T.-P. Fries.  Isogeometric Boundary Element analysis with elasto-plastic inclusions. Part 2: 3-D problems.  Comput. Methods Appl. Mech. Engrg., 315:418–433, 2017.
  * [6] G. Beer, B. Marussig, and C. Duenser.  The Isogeometric Boundary Element Method.  Springer, 2020.
  * [7] C. Bracco, C. Giannelli, F. Mazzia, and A. Sestini.  Bivariate hierarchical Hermite spline quasi–interpolation.  BIT, 56:1165–1188, 2016.
  * [8] F. Calabrò, A. Falini, M. L. Sampoli, and A. Sestini.  Efficient quadrature rules based on spline quasi-interpolation for application to IgA-BEMs.  J Comput Appl Math, 338:153–167, 2018.
  * [9] E. Cohen, T. Martin, R. M. Kirby, T. Lyche, and R. F. Riesenfeld.  Analysis-aware modeling: Understanding quality considerations in modeling for isogeometric analysis.  Comput. Methods Appl. Mech. Engrg., 199(5-8):334–356, 2010.
  * [10] L. Coox, O. Atak, D. Vandepitte, and W. Desmet.  An isogeometric indirect boundary element method for solving acoustic problems in open-boundary domains.  Comput. Methods Appl. Mech. Engrg., 316:186–208, 2017.
  * [11] M. Costabel.  Principles of boundary element methods.  Techn. Hochsch., Fachbereich Mathematik, 1986.
  * [12] M. Costabel, F. Penzel, and R. Schneider.  Error analysis of a boundary element collocation method for a screen problem in $r^{3}$.  Math. Comp., 58(198):575–586, 1992.
  * [13] J. A. Cottrell, T. J. R. Hughes, and Y. Bazilevs.  Isogeometric analysis: toward integration of CAD and FEA.  John Wiley & Sons, 2009.
  * [14] J. Dölz, H. Harbrecht, S. Kurz, S. Schöps, and F. Wolf.  A fast isogeometric BEM for the three dimensional Laplace and Helmholtz problems.  Comput. Methods Appl. Mech. Engrg., 330:83–101, 2018.
  * [15] A. Falini, C. Giannelli, T. Kanduč, M. L. Sampoli, and A. Sestini.  An adaptive IgA-BEM with hierarchical B-splines based on quasi-interpolation quadrature schemes.  Int. J. Numer. Methods Engrg., 117(10):1038–1058, 2019.
  * [16] A. Falini, T. Kanduč, M. L. Sampoli, and A. Sestini.  Cubature rules based on bivariate spline quasi-interpolation for weakly singular integrals.  In G. E. Fasshauer, M. Neamtu, and L. L. Schumaker, editors, Approximation Theory XVI. AT 2019, volume 336 of Springer in Mathematics & Statistics, pages 73–86. Springer Cham., 2020.
  * [17] A. I. Ginnis, K. V. Kostas, C. G. Politis, P. D. Kaklis, K. A. Belibassakis, Th. P. Gerostathis, M. A. Scott, and T. J. R. Hughes.  Isogeometric boundary-element analysis for the wave-resistance problem using T-splines.  Comput. Methods Appl. Mech. Engrg., 279:425–439, 2014.
  * [18] Y. P. Gong, C. Y. Dong, and X. C. Qin.  An isogeometric boundary element method for three dimensional potential problems.  J. Comput. Appl. Math., 313:454–468, 2017.
  * [19] T. J. R. Hughes, J. A. Cottrell, and Y. Bazilevs.  Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement.  Comput. Methods Appl. Mech. Engrg., 194(39–41):4135–4195, 2005\.
  * [20] F. Mazzia and A. Sestini.  The BS class of Hermite spline quasi-interpolants on nonuniform knot distributions.  BIT, 49(3):611–628, 2009.
  * [21] F. Mazzia and A. Sestini.  Quadrature formulas descending from BS Hermite spline quasi-interpolation.  J. Comput. Appl. Math., 236:4105–4118, 2012.
  * [22] K. Mørken.  Some identities for products and degree raising of splines.  Constr. Approx., 7:195–208, 1991.
  * [23] B. H. Nguyen, H. D. Tran, C. Anitescu, X. Zhuang, and T. Rabczuk.  An isogeometric symmetric Galerkin boundary element method for two–dimensional crack problems.  Comput. Methods Appl. Mech. Engrg., 306:252–275, 2016.
  * [24] B. H. Nguyen, X. Zhuang, P. Wriggers, T. Rabczuk, M. E. Mear, and H. D. Tran.  Isogeometric symmetric galerkin boundary element method for three-dimensional elasticity problems.  Comput. Methods Appl. Mech. Engrg., 323:132–150, 2017.
  * [25] M. Pan, F. Chen, and W. Tong.  Volumetric spline parameterization for isogeometric analysis.  Comput. Methods Appl. Mech. Engrg., 359:112769, 2020.
  * [26] S. A. Sauter and C. Schwab.  Boundary element methods, volume 39 of Springer Series in Computational Mathematics.  Springer-Verlag, Berlin, Heidelberg, 2011.
  * [27] R. N. Simpson, S. P. A. Bordas, J. Trevelyan, and T. Rabczuk.  A two-dimensional Isogeometric Boundary Element Method for elastostatic analysis.  Comput. Methods Appl. Mech. Engrg., 209–212:87–100, 2012.
  * [28] A. Sutradhar, G. Paulino, and L. J. Gray.  Symmetric Galerkin boundary element method.  Springer Science & Business Media, 2008.
  * [29] M. Taus, G. J. Rodin, and T. J. R. Hughes.  Isogeometric analysis of boundary integral equations: High-order collocation methods for the singular and hyper-singular equations.  Math. Models and Methods in Appl. Sci., 26(8):1447–1480, 2016.
  * [30] D. Toshniwal, H. Speleers, R. R. Hiemstra, and T. J. R. Hughes.  Multi-degree smooth polar splines: A framework for geometric modeling and isogeometric analysis.  Comput. Methods Appl. Mech. Engrg., 316:1005–10061, 2017.
  * [31] Y. J. Wang and D. J. Benson.  Multi-patch nonsingular isogeometric boundary element analysis in 3d.  Comput. Methods Appl. Mech. Engrg., 293:71–91, 2015.