# Isoparametric singularity extraction technique
for 3D potential problems in BEM

Tadej Kanduč<EMAIL_ADDRESS>Faculty of Mathematics and Physics
and Faculty of Mechanical Engineering, University of Ljubljana, Ljubljana,
Slovenia

###### Abstract

To solve boundary integral equations for potential problems using collocation
Boundary Element Method (BEM) on smooth curved 3D geometries, an analytical
singularity extraction technique is employed. By adopting the isoparametric
approach, curved geometries that are represented by mapped rectangles or
triangles from the parametric domain are considered. The singularity
extraction on the governing singular integrals can be performed either as an
operation of subtraction or division, each having some advantages.

A particular series expansion of a singular kernel about a source point is
investigated. The series in the intrinsic coordinates consists of functions of
a type $R^{p}x^{q}y^{r}$, where $R$ is a square root of a quadratic bivariate
homogeneous polynomial, corresponding to the first fundamental form of a
smooth surface, and $p,q,r$ are integers, satisfying $p\leq-1$ and $q,r\geq
0$. By extracting more terms from the series expansion of the singular kernel,
the smoothness of the regularized kernel at the source point can be increased.
Analytical formulae for integrals of such terms are obtained from
antiderivatives of $R^{p}x^{q}y^{r}$, using recurrence formulae, and by
evaluating them at the edges of rectangular or triangular parametric domains.

Numerical tests demonstrate that the singularity extraction technique can be a
useful prerequisite for a numerical quadrature scheme to obtain accurate
evaluations of the governing singular integrals in 3D collocation BEM.

###### keywords:

Analytical integration , Singular integral , Singularity extraction , Boundary
element method, Isogeometric analysis

##  1 Introduction

Isogeometric Analysis (IGA) [32, 12] has in recent years brought a renewed
interest in developing highly accurate simulation models in Boundary Element
Method (BEM); e.g., see advances in local and adaptive mesh refinements [27,
23, 24], trimming [5], fast matrix formation [14], fluid–structure interaction
[31], fracture [43], shape optimization [37], acoustics and linear elasticity
problems [49, 53, 50], and references therein. A predominant setting in IGA
combines exact representation of the geometry using CAD standards such as
polynomial and rational splines, and the same type of basis functions to
represent the approximate solutions of PDEs. To fully profit by high
convergence orders of approximate solutions of isogeometric methods, the
construction of efficient numerical rules for (singular) integrals is one of
the most crucial challenges to be addressed in IGA-BEM [1, 2, 17, 30, 51].
Although the Gaussian quadratures are usually considered as the optimal choice
for smooth integrands with polynomial bases, there are better alternatives for
integrals with B-splines bases, which exploit the interelement continuity of
splines; see the review paper [8].

A special care is needed to treat singular integrals with numerical schemes –
a continuous research in the last decades (mostly outside isogeometric
context) to design accurate, efficient and easy to implement integration rules
produced a vast number of approaches, including explicit formulae for flat
domains [22, 33], radial integration [25, 26], (near) singularity smoothing
[34, 40, 48, 52, 54], true desingularisation [36] and other techniques [9];
see also references therein. Coordinate transformation, such as polar [29, 35,
45, 51] and Duffy transformation [15, 6, 42, 44], can reduce and simplify the
nature of the singularity. An established technique, used in various forms, is
the singularity extraction, in which the singular part of the integral is
evaluated analytically [3, 22, 26, 29, 33]. In the preceding research to this
paper, a promising approach for curved isoparametric boundaries for 2D
problems combines an elegant singularity extraction and a local spline quasi-
interpolation operators [2, 7, 17]. For this type of quadrature schemes the
optimal convergence orders of the approximate solution can be recovered with a
small numbers of quadrature nodes [19].

In this paper we study a singularity extraction for singular integrals for 3D
potential problems. Focusing on isogeometric collocation BEM, curved boundary
domains are represented as a set of rectangles or triangles in the parametric
domain, mapped to the physical space using B-spline (or NURBS) functions. As a
proof of concept, the main focus is on weakly and nearly singular integrals,
appearing for Laplace problems, although the ideas could be applied to other
types of integrals as well. The singularity extraction utilizes a particular
series expansion of the singular kernel about a source point in the parametric
domain – to the best of author’s knowledge this type of series expansions have
not been studied before in this form in BEM. No distortion of the integration
domain occurs since the extraction is applied directly in the initial
parametric domain; this is an important accuracy and efficiency feature since
the same quadrature nodes can be used for regular and singular integrals and
for several basis functions and source points. The expansion depends on the
(higher) derivatives of the geometry parameterization, which are considered to
be easily computable in IGA setting. For each summand in the series,
multiplied with a polynomial basis function, recursive formulae to evaluate
its double integral analytically are provided. The smoothness of the
regularized kernel at the source point is controlled by the number of
considered terms in the series. An additional constant cost to evaluate more
terms can be outweighed by the higher accuracy of a quadrature for regularized
integrals – this is especially a preferred trade-off for isogeometric methods
with high convergence orders. The developed singularity extraction has already
been applied in BEM for Laplace problems [18, 16] and Helmholtz equation can
be dealt similarly [20], since the appearing singularities are essentially the
same.

The structure of the paper is as follows. In Section 2.1 we recall a basic
setting of integral formulation for potential problems – it serves as a short
introduction and motivation to Section 2.2, where series expansions of kernels
for the single and double layer potentials are derived. For completeness, we
briefly outline the considered numerical integration schemes for the governing
integrals in BIE in Section 2.3 and we provide some directions for
implementation speedups for the analytical evaluation of the studied singular
integrals in 2.4. In 2.5, smoothness of a more general regularized kernel is
studied with respect to the two regularization techniques: subtraction and
division. In Section 3 we perform tests of numerical integrations for the two
common singular kernels for Laplace problems and for the two regularization
techniques. In appendix we collect basic properties of functions
$R^{p}x^{q}y^{r}$ and derive the recursive formulae to compute their
antiderivatives.

##  2 Singularity extraction

In this section we demonstrate how to apply a singularity extraction technique
on weakly singular kernels for the single and double layer potential for the
3D Laplace problems (see [20] how to rewrite the kernels in Helmholtz
equation). As a first step, a series expansion of the singular kernels about
the source point in intrinsic coordinates is derived. A truncated series is an
approximation of the singular kernel and it is used either in a singularity
subtraction or in a singularity division to obtain regularized kernels. The
remaining regular parts of the integrands can be locally approximated with a
suitable polynomial. Therefore, the singular integrals in the underlying
Boundary Integral Equations (BIE) are approximated as a finite sum of
integrals, described in Section A. At the end of the section we also analyze
the smoothness of the derived regularized kernels.

###  2.1 Integral formulation of the potential problem

To be concise, let us consider 3D potential problems described by the Laplace
equation with Dirichlet boundary conditions on finite volumes $\Omega$ with
closed smooth boundary surfaces $\Gamma$,

$\displaystyle\left\\{\begin{array}[]{ll}\Delta u=0&{\rm in}\;\Omega,\\\
u=u_{D}&{\rm on}\;\Gamma,\end{array}\right.$ (3)

where the solution $u$ belongs to the Sobolev space $H^{1}(\Omega)$ and the
Dirichlet boundary datum $u_{D}$ is in $H^{1/2}(\Gamma)$, the trace space of
$H^{1}(\Omega)$. For more details we refer to [4, 11, 46].

To solve (3) using the so-called direct approch we rewrite the problem into
the following BIE (the Symm’s integral equation)

$\displaystyle\int_{\Gamma}{\mathcal{G}}({\boldsymbol{x}},{\boldsymbol{y}})\phi({\boldsymbol{y}})\mathop{d\Gamma_{{\boldsymbol{y}}}}=\frac{1}{2}u_{D}({\boldsymbol{x}})+\int_{\Gamma}{\mathcal{H}}({\boldsymbol{x}},{\boldsymbol{y}})u_{D}({\boldsymbol{y}})\mathop{d\Gamma_{{\boldsymbol{y}}}},\qquad{\boldsymbol{x}}\in\Gamma,$
(4)

where the function $\phi$ denotes the unknown flux of $u$ and belongs to
$H^{-1/2}(\Gamma)$, the dual space of $H^{1/2}(\Gamma)$ (duality is defined
with respect to the usual $L^{2}(\Gamma)$-scalar product). The weakly singular
kernels for the single and the double layer potentials are

$\displaystyle{\mathcal{G}}({\boldsymbol{x}},{\boldsymbol{y}})=\frac{1}{4\pi}\frac{1}{\|{\boldsymbol{x}}-{\boldsymbol{y}}\|_{2}},\qquad{\mathcal{H}}({\boldsymbol{x}},{\boldsymbol{y}})=\frac{\partial{\mathcal{G}}({\boldsymbol{x}},{\boldsymbol{y}})}{\partial{\boldsymbol{n}}_{{\boldsymbol{y}}}}=\frac{1}{4\pi}\frac{({\boldsymbol{x}}-{\boldsymbol{y}})^{\top}}{\|{\boldsymbol{x}}-{\boldsymbol{y}}\|_{2}^{3}}{\boldsymbol{n}}({\boldsymbol{y}}),$
(5)

and ${\boldsymbol{n}}$ denotes the outward unit normal vector.

Following the isoparametric approach, the boundary $\Gamma$ is represented as
a mapped domain $\hat{\Gamma}$ in ${\mathop{{\rm I}\kern-1.99997pt{\rm
R}}\nolimits}^{2}$,

$\displaystyle{\boldsymbol{F}}:\hat{\Gamma}\to{\mathop{{\rm
I}\kern-1.99997pt{\rm
R}}\nolimits}^{3},\;{\boldsymbol{t}}\mapsto{\boldsymbol{F}}({\boldsymbol{t}}),$

where the geometry mapping ${\boldsymbol{F}}$ is described in terms of
bivariate splines or NURBS [13, 21, 47] in parametric coordinates
${\boldsymbol{t}}:=(t_{1},t_{2})$. By applying the collocation discretization,
we search for the approximate solution $\phi_{h}$ of (4) in the discretization
space, spanned by basis functions $B_{j}:\hat{\Gamma}\to{\mathop{{\rm
I}\kern-1.99997pt{\rm R}}\nolimits}$, $j=1,2,\dots,N_{\rm DOF}$,

$\displaystyle\phi_{h}({\boldsymbol{x}}):=\sum_{j=1}^{N_{\rm
DOF}}\alpha_{j}(B_{j}\circ{\boldsymbol{F}}^{-1})({\boldsymbol{x}}).$

By fixing suitable collocation sites ${\boldsymbol{x}}_{i}$, for
$i=1,2,\dots,N_{\rm DOF}$, the unknown entries in the vector
$\boldsymbol{\alpha}=(\alpha_{1},\dots,\alpha_{N_{\rm DOF}})^{\top}$ are
obtained by solving a square linear system,

$\displaystyle A\boldsymbol{\alpha}=\boldsymbol{\beta}.$

The system matrix $A$ is fully populated and its entries $a_{ij}$ are

$\displaystyle a_{ij}$
$\displaystyle=\frac{1}{4\pi}\int_{\Gamma}\frac{1}{\|{\boldsymbol{x}}_{i}-{\boldsymbol{y}}\|_{2}}(B_{j}\circ{\boldsymbol{F}}^{-1})({\boldsymbol{y}})\mathop{d\Gamma_{{\boldsymbol{y}}}}=\frac{1}{4\pi}\int_{\hat{\Gamma}}\frac{1}{\|{\boldsymbol{F}}({\boldsymbol{s}}_{i})-{\boldsymbol{F}}({\boldsymbol{t}})\|_{2}}B_{j}({\boldsymbol{t}})J({\boldsymbol{t}})\mathop{d\boldsymbol{t}},$

where $J:=\|D_{1}{\boldsymbol{F}}\times D_{2}{\boldsymbol{F}}\|_{2}$ is the
area of the infinitesimal surface element,
$D_{1}:=\partial{\boldsymbol{F}}/{\partial
t_{1}},D_{2}:=\partial{\boldsymbol{F}}/{\partial t_{2}}$, and
${\boldsymbol{F}}({\boldsymbol{s}}_{i})={\boldsymbol{x}}_{i}$,
${\boldsymbol{F}}({\boldsymbol{t}})={\boldsymbol{y}}$. Entries of the right-
hand side vector $\boldsymbol{\beta}$ are

$\displaystyle\beta_{i}$
$\displaystyle=\frac{1}{2}(u_{D}\circ{\boldsymbol{F}})({\boldsymbol{s}}_{i})+\frac{1}{4\pi}\int_{\hat{\Gamma}}{\mathcal{H}}({\boldsymbol{F}}({\boldsymbol{s}}_{i}),{\boldsymbol{F}}({\boldsymbol{t}}))\,J({\boldsymbol{t}})\,(u_{D}\circ{\boldsymbol{F}})({\boldsymbol{t}})\mathop{d\boldsymbol{t}}.$

Integral in $a_{ij}$ is classified singular if the source point
${\boldsymbol{s}}_{i}$ is inside the support of $B_{j}$. If
${\boldsymbol{s}}_{i}$ is outside the support but close to it, we declare the
integral nearly singular. In the remaining case we say the integral is
regular.

###  2.2 Series expansion of a singular kernel

In this subsection we derive a series expansion of kernels ${\mathcal{G}}$ and
${\mathcal{H}}$ for fixed source point ${\boldsymbol{s}}$. By restricting the
kernel to a line segment, the singular part of the kernel is easily decoupled
from the regular part – for the latter part we then apply the Taylor series
expansion. The procedure shares some similarity to polar and Duffy
transformation and to radial integration approach by Gao [25], with the
difference that in our case the final expression of the series is written
again in the starting intrinsic coordinates. Each summand in the series is a
function of a type, described in Section A.2, thus a closed form expression
for its integral exists.

####  2.2.1 Series expansion of ${\mathcal{G}}$

As a proof of concept we focus on the simplest kernel, ${\mathcal{G}}$ in (5).
Let us write ${\boldsymbol{x}}={\boldsymbol{F}}({\boldsymbol{s}})$,
${\boldsymbol{y}}={\boldsymbol{F}}({\boldsymbol{t}})$ and let us assume
${\boldsymbol{F}}$ to be analytical (or sufficiently smooth) near fixed
${\boldsymbol{x}}={\boldsymbol{F}}({\boldsymbol{s}})$. For a fixed
${\boldsymbol{s}}$ we can express kernel ${\mathcal{G}}$,
${\mathcal{G}}_{\boldsymbol{s}}({\boldsymbol{z}}):={\mathcal{G}}({\boldsymbol{F}}({\boldsymbol{s}}),{\boldsymbol{F}}({\boldsymbol{s}}-{\boldsymbol{z}}))$,
in local coordinates
${\boldsymbol{z}}:=(z_{1},z_{2}):={\boldsymbol{s}}-{\boldsymbol{t}}$. Taylor
expansion of ${\boldsymbol{F}}$ about ${\boldsymbol{s}}$ can be compactly
written as

$\displaystyle{\boldsymbol{F}}({\boldsymbol{t}})={\boldsymbol{F}}({\boldsymbol{s}})+\sum_{i=1}^{\infty}{\boldsymbol{D}}_{i}{\boldsymbol{F}}({\boldsymbol{t}}),\qquad{\boldsymbol{D}}_{i}{\boldsymbol{F}}({\boldsymbol{t}}):=\sum_{|\alpha|=i}{\boldsymbol{a}}_{\alpha}{\boldsymbol{z}}^{\alpha},\qquad{\boldsymbol{a}}_{\alpha}:=\frac{(-1)^{|\alpha|}}{\alpha!}D^{\alpha}{\boldsymbol{F}}({\boldsymbol{s}}).$

Here we adopt the multi-index notation. Specifically,
$\alpha=(\alpha_{1},\alpha_{2})$, $\alpha_{1},\alpha_{2}\in\mathbb{N}_{0}$ and
$|\alpha|=\alpha_{1}+\alpha_{2}$, $\alpha!=\alpha_{1}!\alpha_{2}!$. The
partial derivative operator $D^{\alpha}$ acts on each component of a vector
function ${\boldsymbol{F}}$ separately and on both variables,

$\displaystyle
D^{\alpha}=D_{1}^{\alpha_{1}}D_{2}^{\alpha_{2}}=\frac{\partial^{|\alpha|}}{\partial
t_{1}^{\alpha_{1}}\partial t_{2}^{\alpha_{2}}}.$

Also we compactly write the scalar
${\boldsymbol{z}}^{\alpha}=z_{1}^{\alpha_{1}}z_{2}^{\alpha_{2}}$.

From the expression

$\displaystyle{\boldsymbol{F}}({\boldsymbol{s}})-{\boldsymbol{F}}({\boldsymbol{t}})=-\sum_{|\beta|\geq
1}{\boldsymbol{a}}_{\beta}{\boldsymbol{z}}^{\beta}$

we can easily express the kernel ${\mathcal{G}}_{{\boldsymbol{s}}}$ as a power
series of ${\boldsymbol{z}}$ about $(0,0)$,

$\displaystyle
4\pi{\mathcal{G}}({\boldsymbol{F}}({\boldsymbol{s}}),{\boldsymbol{F}}({\boldsymbol{t}}))$
$\displaystyle=\frac{1}{\|{\boldsymbol{F}}({\boldsymbol{s}})-{\boldsymbol{F}}({\boldsymbol{t}})\|_{2}}=\left(\sum_{{|\alpha|\geq
2}}c_{\alpha}{\boldsymbol{z}}^{\alpha}\right)^{-\frac{1}{2}},$ (6)

where

$\displaystyle c_{\alpha}=\sum_{\begin{subarray}{c}\beta+\gamma=\alpha\\\
|\beta|\geq 1,|\gamma|\geq
1\end{subarray}}{\boldsymbol{a}}_{\beta}^{\top}{\boldsymbol{a}}_{\gamma}$ (7)

are the coefficients in the series expansion of
$\|{\boldsymbol{F}}({\boldsymbol{s}})-{\boldsymbol{F}}({\boldsymbol{t}})\|_{2}^{2}$
about ${\boldsymbol{s}}$ in local variable ${\boldsymbol{z}}$. From (7) it is
easy to check that coefficient $c_{\alpha}$ is a linear combination of
products of derivatives of ${\boldsymbol{F}}$, where the total degree of the
derivatives is at most $|\alpha|-1$.

In an ideal world one would probably just truncate the series in (6) and
integrate the obtained function, involving the Taylor polynomial.
Unfortunately, closed forms for these types of integrals are not known in
general, therefore a further transformation of the expression is needed. In
the next step we derive a Taylor expansion of the regular part of (6),
restricted to an arbitrary line $(z_{1},vz_{1})$, $v\in\mathbb{R}$. From the
parameter $v$ depended linear expansion it is relatively straightforward to
derive a truncated series expansion of ${\mathcal{G}}_{\boldsymbol{s}}$ that
can be integrated analytically.

Let $L_{v}:=\\{(z_{1},vz_{1}),z_{1}\in(-\delta,\delta)\\}$ be a line segment
passing the origin $(0,0)$, where $v\in\mathbb{R}$ is a fixed parameter and
$\delta$ is sufficiently small. With ${\mathcal{G}}_{v}$ let us write the
kernel ${\mathcal{G}}$, restricted on $L_{v}$ and written as a function of
$z_{1}$. By multiplying the expression (6) with $|z_{1}|$ we obtain

$\displaystyle
r_{v}(z_{1}):=|z_{1}|\,4\pi{\mathcal{G}}_{v}(z_{1})=|z_{1}|\left(\sum_{|\alpha|\geq
2}c_{\alpha}z_{1}^{\alpha_{1}}(vz_{1})^{\alpha_{2}}\right)^{-\frac{1}{2}}=\left(\sum_{|\alpha|\geq
2}c_{\alpha}v^{\alpha_{2}}z_{1}^{|\alpha|-2}\right)^{-\frac{1}{2}},$ (8)

and the newly defined function $r_{v}$ is regular (it holds
$\sum_{|\alpha|=2}c_{\alpha}v^{\alpha_{2}}>0$ due to the first fundamental
form at a regular point on the surface). That way we separate the singular
part $|z_{1}|^{-1}$ of the kernel ${\mathcal{G}}_{v}$ from the regular part
$r_{v}(z_{1})$ ,

$\displaystyle 4\pi{\mathcal{G}}_{v}(z_{1})=|z_{1}|^{-1}r_{v}(z_{1}).$ (9)

Since $r_{v}$ is a regular function in $z_{1}$, we can replace it with its
Taylor series expansion in $z_{1}$ about origin 0 using expression (8),

$\displaystyle
r_{v}(z_{1})=\sum_{\ell=0}^{\infty}\frac{r_{v}^{(\ell)}(0)}{\ell!}z_{1}^{\ell}=\big{(}P_{v,2}\big{)}^{-\frac{1}{2}}-\frac{1}{2}\big{(}P_{v,2}\big{)}^{-\frac{3}{2}}P_{v,3}\,z_{1}+\bigg{(}-\big{(}P_{v,2}\big{)}^{-\frac{3}{2}}P_{v,4}+\frac{3}{4}\big{(}P_{v,2}\big{)}^{-\frac{5}{2}}\big{(}P_{v,3}\big{)}^{2}\bigg{)}\frac{z_{1}^{2}}{2}+\dots$
(10)

where $P_{v,\ell}:=\sum_{|\alpha|=\ell}c_{\alpha}v^{\alpha_{2}}$ is a
polynomial of degree $\ell$ in $v$.

The obtained expression for $4\pi{\mathcal{G}}_{v}$ is valid for any $v$, we
insert back $\displaystyle{v=\frac{z_{2}}{z_{1}}}$ into (9) using the
expansion (10) and get the sought series expansion of
${\mathcal{G}}_{\boldsymbol{s}}$,

$\displaystyle{\mathcal{G}}_{\boldsymbol{s}}({\boldsymbol{z}})=\frac{1}{4\pi}\left[R({\boldsymbol{z}})^{-1}-\frac{1}{2}R({\boldsymbol{z}})^{-3}P_{3}({\boldsymbol{z}})+\left(-\frac{1}{2}R({\boldsymbol{z}})^{-3}P_{4}({\boldsymbol{z}})+\frac{3}{8}R({\boldsymbol{z}})^{-{5}}P_{3}({\boldsymbol{z}})^{2}\right)+\dots\right],$
(11)

where
$R({\boldsymbol{z}}):=R(c_{(2,0)},c_{(1,1)},c_{(0,2)},{\boldsymbol{z}})=(\sum_{|\alpha|=2}c_{\alpha}{\boldsymbol{z}}^{\alpha})^{1/2}$,
$P_{\ell}({\boldsymbol{z}}):=\sum_{|\alpha|=\ell}c_{\alpha}{\boldsymbol{z}}^{\alpha}$
is the bivariate homogeneous polynomial of degree $\ell$, and coefficients
$c_{\alpha}$ are defined in (7). Observe that in (11) there is no problem in
the definition when $z_{1}=0$ and $z_{2}\neq 0$; it is only a limitation of
the definition of ${\mathcal{G}}_{v}$.

For the sake of easier notation and vocabulary, we consider an $\ell$-th term
of (11) a weighted sum of all entries $R^{p}z_{1}^{q}z_{2}^{r}$, that satisfy
$\zeta(R^{p}z_{1}^{q}z_{2}^{r})=\ell-2$. Indeed, compactly we could express
(11) as

$\displaystyle{\mathcal{G}}_{\boldsymbol{s}}=\frac{1}{4\pi}\sum_{\ell=1,2,\dots}R^{-2\ell+1}P_{3\ell-3}^{[\ell]},$
(12)

where $P_{3\ell-3}^{[\ell]}$ are appropriate homogeneous polynomials of degree
$3\ell-3$ and clearly for the $\ell$-th term it holds
$\zeta(R^{-2\ell+1}P_{3\ell-3}^{[\ell]})=\ell-2$. In this notation we assume
that there are no common factors between polynomials $R^{2}$ and
$P_{3\ell-3}^{[\ell]}$.

With ${\mathcal{G}}_{{\boldsymbol{s}},n}$ we denote a kernel, obtained from
${\mathcal{G}}_{\boldsymbol{s}}$ by truncating the series in (11) after the
$n$-th term. The kernel ${\mathcal{G}}_{{\boldsymbol{s}},n}$ is a sum of $n$
addends and double integral of each of the terms is a sum of functions of a
type $I_{p,q,r}(z_{1},z_{2})$, defined in (33). As we show in Section A.2,
these type of integrals have closed form expressions.

######  Remark 1.

If the geometry is locally flat and parameterized by a linear map, then the
cancellation of the singularity with the first term in the series expansion is
exact,
$\rho:={\mathcal{G}}_{\boldsymbol{s}}-{\mathcal{G}}_{{\boldsymbol{s}},1}\equiv
0$ and
$\rho:={\mathcal{G}}_{\boldsymbol{s}}/{\mathcal{G}}_{{\boldsymbol{s}},1}\equiv
1$ , for singularity subtraction and division, respectively (see Section 2.3
for more details about utilization of the singularity extraction in numerical
integration). Thus, the procedure generalizes the common singularity
extraction formulae for flat surfaces.

####  2.2.2 Series expansion of ${\mathcal{H}}$

We proceed similarly as for the kernel ${\mathcal{G}}$. Since we want to
express the kernel ${\mathcal{H}}$ in the intrinsic coordinates from the
information of the geometry mapping ${\boldsymbol{F}}$, it is more convenient
to consider the kernel $\bar{\mathcal{H}}:={\mathcal{H}}J$, where
$J({\boldsymbol{t}})=\|D_{1}{\boldsymbol{F}}({\boldsymbol{t}})\times
D_{2}{\boldsymbol{F}}({\boldsymbol{t}})\|_{2}$,

$\displaystyle\bar{\mathcal{H}}({\boldsymbol{F}}({\boldsymbol{s}}),{\boldsymbol{F}}({\boldsymbol{t}}))={\mathcal{H}}({\boldsymbol{F}}({\boldsymbol{s}}),{\boldsymbol{F}}({\boldsymbol{t}}))\,J({\boldsymbol{t}})=\left(\frac{{\boldsymbol{F}}({\boldsymbol{s}})-{\boldsymbol{F}}({\boldsymbol{t}})}{\|{\boldsymbol{F}}({\boldsymbol{s}})-{\boldsymbol{F}}({\boldsymbol{t}})\|^{3}}\right)^{\top}(D_{1}{\boldsymbol{F}}({\boldsymbol{t}})\times
D_{2}{\boldsymbol{F}}({\boldsymbol{t}})).$ (13)

The inner product of the two vectors in (13) can be compactly written as a
linear combination of its components. Therefore

$\displaystyle\bar{\mathcal{H}}({\boldsymbol{F}}({\boldsymbol{s}}),{\boldsymbol{F}}({\boldsymbol{t}}))=\frac{\varepsilon_{\ell_{1}\ell_{2}\ell_{3}}\big{(}{\boldsymbol{F}}({\boldsymbol{s}})-{\boldsymbol{F}}({\boldsymbol{t}})\big{)}_{\ell_{1}}\big{(}D_{1}{\boldsymbol{F}}({\boldsymbol{t}})\big{)}_{\ell_{2}}\big{(}D_{2}{\boldsymbol{F}}({\boldsymbol{t}})\big{)}_{\ell_{3}}}{\|{\boldsymbol{F}}({\boldsymbol{s}})-{\boldsymbol{F}}({\boldsymbol{t}})\|^{3}},$
(14)

where we use Levi-Civita symbol $\varepsilon$, components of vectors
$(\bullet)_{\ell_{i}}$, for ${\ell_{i}}=1,2,3$, and Einstein summation
convention.

By fixing the source points ${\boldsymbol{s}}$ we define
$\bar{\mathcal{H}}_{\boldsymbol{s}}({\boldsymbol{z}}):=\bar{\mathcal{H}}({\boldsymbol{F}}({\boldsymbol{s}}),{\boldsymbol{F}}({\boldsymbol{s}}-{\boldsymbol{z}}))$
and the geometric quantities in the kernel can be written in the Taylor series
about ${\boldsymbol{s}}$ in the local variable
${\boldsymbol{z}}:={\boldsymbol{s}}-{\boldsymbol{t}}$,

$\displaystyle{\boldsymbol{F}}({\boldsymbol{s}})-{\boldsymbol{F}}({\boldsymbol{t}})$
$\displaystyle=-\sum_{|\beta|\geq
1}{\boldsymbol{a}}_{\beta}{\boldsymbol{z}}^{\beta},$ $\displaystyle
D_{1}{\boldsymbol{F}}({\boldsymbol{t}})$
$\displaystyle=D_{1}\bigg{(}{\boldsymbol{F}}({\boldsymbol{s}})+\sum_{|\gamma|\geq
1}{\boldsymbol{a}}_{\gamma}{\boldsymbol{z}}^{\gamma}\bigg{)}=-\sum_{|\gamma|\geq
0}(\gamma_{1}+1){\boldsymbol{a}}_{(\gamma_{1}+1,\gamma_{2})}{\boldsymbol{z}}^{\gamma},$
$\displaystyle D_{2}{\boldsymbol{F}}({\boldsymbol{t}})$
$\displaystyle=D_{2}\bigg{(}{\boldsymbol{F}}({\boldsymbol{s}})+\sum_{|\delta|\geq
1}{\boldsymbol{a}}_{\delta}{\boldsymbol{z}}^{\delta}\bigg{)}=-\sum_{|\delta|\geq
0}(\delta_{2}+1){\boldsymbol{a}}_{(\delta_{1},\delta_{2}+1)}{\boldsymbol{z}}^{\delta},$

where $\gamma=(\gamma_{1},\gamma_{2})$, $\delta=(\delta_{1},\delta_{2})$.
Therefore the kernel $\bar{\mathcal{H}}_{\boldsymbol{s}}$ can be compactly
written as

$\displaystyle\bar{\mathcal{H}}_{\boldsymbol{s}}=\frac{\displaystyle\sum_{{|\alpha|\geq
2}}d_{\alpha}{\boldsymbol{z}}^{\alpha}}{\displaystyle\left(\sum_{{|\alpha|\geq
2}}c_{\alpha}{\boldsymbol{z}}^{\alpha}\right)^{\frac{3}{2}}},$

where $c_{\alpha}$ are defined as in (7) and

$\displaystyle
d_{\alpha}=-\sum_{\begin{subarray}{c}\beta+\gamma+\delta=\alpha\\\ |\beta|\geq
1,|\gamma|\geq 0,|\delta|\geq
0\end{subarray}}\varepsilon_{\ell_{1}\ell_{2}\ell_{3}}\big{(}{\boldsymbol{a}}_{\beta}\big{)}_{\ell_{1}}\big{(}(\gamma_{1}+1){\boldsymbol{a}}_{(\gamma_{1}+1,\gamma_{2})}\big{)}_{\ell_{2}}\big{(}(\delta_{2}+1){\boldsymbol{a}}_{(\delta_{1},\delta_{2}+1)}\big{)}_{\ell_{3}}.$

Note that $d_{\alpha}=0$ for $|\alpha|=1$; namely

$\displaystyle d_{(1,0)}$
$\displaystyle=\varepsilon_{\ell_{1}\ell_{2}\ell_{3}}\big{(}{\boldsymbol{a}}_{(1,0)}\big{)}_{\ell_{1}}\big{(}{\boldsymbol{a}}_{(1,0)}\big{)}_{\ell_{2}}\big{(}{\boldsymbol{a}}_{(0,1)}\big{)}_{\ell_{3}}=\left|\begin{matrix}{\boldsymbol{a}}_{(1,0)}&{\boldsymbol{a}}_{(1,0)}&{\boldsymbol{a}}_{(0,1)}\end{matrix}\right|=0,$
$\displaystyle d_{(0,1)}$
$\displaystyle=\varepsilon_{\ell_{1}\ell_{2}\ell_{3}}\big{(}{\boldsymbol{a}}_{(0,1)}\big{)}_{\ell_{1}}\big{(}{\boldsymbol{a}}_{(1,0)}\big{)}_{\ell_{2}}\big{(}{\boldsymbol{a}}_{(0,1)}\big{)}_{\ell_{3}}=\left|\begin{matrix}{\boldsymbol{a}}_{(0,1)}&{\boldsymbol{a}}_{(1,0)}&{\boldsymbol{a}}_{(0,1)}\end{matrix}\right|=0,$

since in both cases the involved three vectors are linearly dependent.

Again, let $L_{v}:=\\{(z_{1},vz_{1}),z_{1}\in(-\delta,\delta)\\}$ for fixed
$v\in\mathbb{R}$, $\delta$ is sufficiently small and let
$\bar{\mathcal{H}}_{v}$ be the kernel $\bar{\mathcal{H}}$, restricted on
$L_{v}$ and written as a function of $z_{1}$. Then

$\displaystyle
r_{v}(z_{1}):=|z_{1}|\,4\pi\bar{\mathcal{H}}_{v}(z_{1})=|z_{1}|\frac{\displaystyle\sum_{{|\alpha|\geq
2}}d_{\alpha}z_{1}^{\alpha_{1}}(vz_{1})^{\alpha_{2}}}{\displaystyle\left(\sum_{{|\alpha|\geq
2}}c_{\alpha}z_{1}^{\alpha_{1}}(vz_{1})^{\alpha_{2}}\right)^{\frac{3}{2}}}=\frac{\displaystyle\sum_{{|\alpha|\geq
2}}d_{\alpha}v^{\alpha_{2}}z_{1}^{|\alpha|-2}}{\displaystyle\left(\sum_{{|\alpha|\geq
2}}c_{\alpha}v^{\alpha_{2}}z_{1}^{|\alpha|-2}\right)^{\frac{3}{2}}}.$ (15)

Its Taylor expansion reads

$\displaystyle r_{v}(z_{1})$
$\displaystyle=\sum_{\ell=0}^{\infty}\frac{r_{v}^{(\ell)}(0)}{\ell!}z_{1}^{\ell}$
$\displaystyle=\big{(}P_{v,2}\big{)}^{-\frac{3}{2}}Q_{v,2}+\left(\big{(}P_{v,2}\big{)}^{-\frac{3}{2}}Q_{v,3}-\frac{3}{2}\big{(}P_{v,2}\big{)}^{-\frac{5}{2}}Q_{v,2}P_{v,3}\right)\,z_{1}$
$\displaystyle+\left(2(P_{v,2})^{-\frac{3}{2}}Q_{v,4}-3(P_{v,2})^{-\frac{5}{2}}Q_{v,3}P_{v,3}-3\big{(}P_{v,2}\big{)}^{-\frac{5}{2}}Q_{v,2}P_{v,4}+\frac{15}{4}\big{(}P_{v,2}\big{)}^{-\frac{7}{2}}Q_{v,2}(P_{v,3})^{2}\right)\frac{1}{2}z_{1}^{2}+\dots$
(16)

where $P_{v,\ell}=\sum_{|\alpha|=\ell}c_{\alpha}v^{\alpha_{2}}$,
$Q_{v,\ell}:=\sum_{|\alpha|=\ell}d_{\alpha}v^{\alpha_{2}}$ are polynomials of
degree $\ell$ in $v$.

Since the obtained expression for $4\pi\bar{\mathcal{H}}_{v}$ is valid for any
$v$,

$\displaystyle 4\pi\bar{\mathcal{H}}_{v}(z_{1})=|z_{1}|^{-1}r_{v}(z_{1}),$
(17)

we insert back $\displaystyle{v=\frac{z_{2}}{z_{1}}}$ into (17) using the
expansion (2.2.2) and get the sought series expansion of
$\bar{\mathcal{H}}_{\boldsymbol{s}}$,

$\displaystyle\bar{\mathcal{H}}_{\boldsymbol{s}}({\boldsymbol{z}})=$
$\displaystyle\frac{1}{4\pi}\Bigg{[}R({\boldsymbol{z}})^{-3}Q_{2}({\boldsymbol{z}})+\left(R({\boldsymbol{z}})^{-3}Q_{3}({\boldsymbol{z}})-\frac{3}{2}R({\boldsymbol{z}})^{-5}Q_{2}({\boldsymbol{z}})P_{3}({\boldsymbol{z}})\right)$
$\displaystyle+\left(R({\boldsymbol{z}})^{-3}Q_{4}({\boldsymbol{z}})-\frac{3}{2}R({\boldsymbol{z}})^{-5}Q_{3}({\boldsymbol{z}})P_{3}({\boldsymbol{z}})-\frac{3}{2}R({\boldsymbol{z}})^{-5}Q_{2}({\boldsymbol{z}})P_{4}({\boldsymbol{z}})+\frac{15}{8}R({\boldsymbol{z}})^{-7}Q_{2}({\boldsymbol{z}})P_{3}({\boldsymbol{z}})^{2}\right)+\dots\Bigg{]},$
(18)

where
$P_{\ell}({\boldsymbol{z}})=\sum_{|\alpha|=\ell}c_{\alpha}{\boldsymbol{z}}^{\alpha}$,
$Q_{\ell}({\boldsymbol{z}}):=\sum_{|\alpha|=\ell}d_{\alpha}{\boldsymbol{z}}^{\alpha}$
are two homogeneous polynomials. Compactly we can write (2.2.2) as

$\displaystyle\bar{\mathcal{H}}_{\boldsymbol{s}}=\frac{1}{4\pi}\sum_{\ell=1,2,\dots}R^{-2\ell-1}P_{3\ell-1}^{[\ell]},$
(19)

where $P_{3\ell-1}^{[\ell]}$ are appropriate homogeneous polynomials of degree
$3\ell-1$, there are no common factors between $R^{2}$ and
$P_{3\ell-1}^{[\ell]}$ and the continuity of the $\ell$-th term is
characterized by $\zeta(R^{-2\ell-1}P_{3\ell-1}^{[\ell]})=\ell-2$.

###  2.3 Numerical integration scheme

By utilizing the kernel expansion from the previous subsection we can now
focus on numerical integration of integrals in the Symm’s integral equation
(4). In this section we consider a more general singular kernel
${\mathcal{K}}_{\boldsymbol{s}}$ that can be written as a similar series
expansion and two regularization approaches are introduced.

Let
${\mathcal{K}}_{\boldsymbol{s}}={\mathcal{K}}_{\boldsymbol{s}}({\boldsymbol{z}})$
be a singular kernel, expressed about a source point ${\boldsymbol{s}}$ in
local intrinsic coordinates ${\boldsymbol{z}}$,

$\displaystyle{\mathcal{K}}_{\boldsymbol{s}}=\sum_{\ell=1,2,\dots}R^{-2(\ell+m_{1})+1}P_{3\ell+m_{2}}^{[\ell]},$
(20)

with regularity $m:=\zeta({\mathcal{K}}_{\boldsymbol{s}})=-2m_{1}+m_{2}+2$ and
integers $m_{2}\geq-3$, $m_{1}\geq 0$. We assume $m<0$ and smaller $m$
corresponds to a stronger type of singularity. Here $P_{3\ell+m_{2}}^{[\ell]}$
are again suitable homogeneous polynomials of degree $3\ell+m_{2}$.
Additionally, let ${\mathcal{K}}_{\boldsymbol{s}}$ be sufficiently smooth for
${\boldsymbol{z}}\in D_{\varepsilon}$, where

$\displaystyle
D_{\varepsilon}:=\\{{\boldsymbol{z}}:0<\|{\boldsymbol{z}}\|_{2}\leq\varepsilon\\}$

is a punctured disk for some fixed $\varepsilon>0$. Let
${\mathcal{K}}_{{\boldsymbol{s}},n}$ be an approximation of
${\mathcal{K}}_{\boldsymbol{s}}$ by truncating the infinite series (20) after
the $n$-th term,

$\displaystyle{\mathcal{K}}_{{\boldsymbol{s}},n}=\sum_{\ell=1,2,\dots,n}R^{-2(\ell+m_{1})+1}P_{3\ell+m_{2}}^{[\ell]}.$
(21)

The singular kernel ${\mathcal{K}}_{\boldsymbol{s}}$ is either regularized by
the subtraction, i.e.,
$\rho:={\mathcal{K}}_{\boldsymbol{s}}-{\mathcal{K}}_{{\boldsymbol{s}},n}$, or
the division, i.e.,
$\rho:={\mathcal{K}}_{\boldsymbol{s}}/{\mathcal{K}}_{{\boldsymbol{s}},n}$. Our
goal is to numerically compute integrals

$\displaystyle\frac{1}{4\pi}\int_{I_{1}}\int_{I_{2}}{\mathcal{K}}_{\boldsymbol{s}}({\boldsymbol{z}})B_{j}({\boldsymbol{s}}-{\boldsymbol{z}})v({\boldsymbol{s}}-{\boldsymbol{z}})\mathop{d\boldsymbol{z}},$
(22)

where $B_{j}$ is a polynomial basis function and $v$ is an auxiliary smooth
function that can include for example the area of the infinitesimal surface
element $J$ or a boundary datum.

  1. Case 1:

If the integral (22) is regular, we can apply one of the common numerical
integration schemes since the integrand is a smooth function.

  2. Case 2:

If the integral (22) is singular or nearly singular, the following
transformation needs to be applied first.

    1. (a)

The first strategy, which is more common in literature, is the singularity
subtraction. By writing
$\rho:={\mathcal{K}}_{\boldsymbol{s}}-{\mathcal{K}}_{{\boldsymbol{s}},n}$, the
integral (22) is expressed as a sum of two integrals

$\displaystyle\frac{1}{4\pi}\int_{I_{1}}\int_{I_{2}}\rho({\boldsymbol{z}})B_{j}({\boldsymbol{s}}-{\boldsymbol{z}})v({\boldsymbol{s}}-{\boldsymbol{z}})+\frac{1}{4\pi}\int_{I_{1}}\int_{I_{2}}{\mathcal{K}}_{{\boldsymbol{s}},n}({\boldsymbol{z}})B_{j}({\boldsymbol{s}}-{\boldsymbol{z}})v({\boldsymbol{s}}-{\boldsymbol{z}})\mathop{d\boldsymbol{z}}.$
(23)

The integrand in the first integral in (23) is sufficiently regular if the
geometry parameterization (hence $v$) and $\rho$ are sufficiently regular. The
smoothness of $\rho$ is studied in the next subsection. Thus, for this type of
integrals we can apply a quadrature scheme for regular integrals. In the
second integral in (23) the regular part
$B_{j}({\boldsymbol{s}}-{\boldsymbol{z}})v({\boldsymbol{s}}-{\boldsymbol{z}})$
can be approximated sufficiently well with an appropriate polynomial $P$ via
least square fitting, quasi-interpolation etc.

    2. (b)

In the regularization process via division step we consider
$\rho:={\mathcal{K}}_{\boldsymbol{s}}/{\mathcal{K}}_{{\boldsymbol{s}},n}$ and
express the integral (22) as

$\displaystyle\frac{1}{4\pi}\int_{I_{1}}\int_{I_{2}}{\mathcal{K}}_{{\boldsymbol{s}},n}({\boldsymbol{z}})B_{j}({\boldsymbol{s}}-{\boldsymbol{z}})v({\boldsymbol{s}}-{\boldsymbol{z}})\rho({\boldsymbol{z}})\mathop{d\boldsymbol{z}}.$
(24)

The regular part of the integrand,
$B_{j}({\boldsymbol{s}}-{\boldsymbol{z}})v({\boldsymbol{s}}-{\boldsymbol{z}})\rho({\boldsymbol{z}})$,
is again approximated with a suitable polynomial $P$.

In both strategies (a), (b), the singular part of the integral (22) is
approximated with

$\displaystyle\frac{1}{4\pi}\int_{I_{1}}\int_{I_{2}}{\mathcal{K}}_{{\boldsymbol{s}},n}({\boldsymbol{z}})P({\boldsymbol{z}})\mathop{d\boldsymbol{z}}$
(25)

for a suitable polynomial $P$. The integral (25) can be expressed as a sum of
fundamental integrals, presented in Section A.2.

######  Remark 2.

In (22), instead of a polynomial basis we can proceed similarly if we consider
$B_{j}$ to play a role of a polynomial spline basis function. In that case we
can decompose the basis function into a sum of polynomial contributions (for
example using a Bézier extraction technique) and from thereon follow a similar
procedure as before.

Alternatively, we can approximate the smooth part of the integrand of (22)
with a spline function and write the product of the two splines as a new
spline function, for example by using formulae in [41] for the tensor product
basis functions. For a successful implementation of this idea to solve Laplace
and Helmholtz boundary integral equations in BEM see [16, 18, 20]. The model
combines a spline quasi-interpolation operator [38, 39], together with a
spline recurrence relation formula [20] to compute the integrals (25), where
$P$ plays a role of a polynomial spline function.

###  2.4 Implementation speedup

The presented singularity extraction technique can only be considered useful
in practice if it can be efficiently implemented in software. For example, to
perform a simulation using BEM, the matrix formation can induce a computation
of thousands or millions of integrals of the type (25) and it can represent
the main bottleneck in the system matrix formation. However, since the
evaluation of these integrals for different source points is embarrassingly
parallelizable, the time complexity of this step can be partially alleviated.
In this subsection we propose an additional reduction of the computational
complexity by accessing specific precomputed integrals from lookup tables.

By exchanging the order of summation and integration in (25), and writing each
polynomial $P$ as a weighted sum of monomials (the weights depend on the local
geometry around the point ${\boldsymbol{s}}$ and can be efficiently computed
in a preprocess step), we are left to analyze how to efficiently integrate
functions $R^{p}x^{q}y^{r}$. Unfortunately, analytical formulae for the
indefinite integrals (33) can have very long symbolic expressions that can
result in undesired lengthy evaluation, especially for higher values of
$-p+q+r$. To overcome this problem, we apply several simplifications of the
expressions to reduce the number of parameters, influencing the integral
values, and precompute values of particular simpler integrals. We can later on
access these values in the stored lookup tables to compute the definite
integrals.

####  2.4.1 Rectangles

The goal is to derive to an efficient procedure to compute the building blocks
of the definite integral of (21),

$\displaystyle
I:=\int_{y_{0}}^{y_{1}}\int_{x_{0}}^{x_{1}}R^{p}x^{q}y^{r}\mathop{dx}\mathop{dy},$

on aligned rectangular domains $D=[x_{0},x_{1}]\times[y_{0},y_{1}]$, where
$R(a,b,c,(x,y))=\sqrt{ax^{2}+bxy+cy^{2}}$, as in (33). Unfortunately, the
integrals depend on too many “continuous” variables (e.g.,
$a,b,c,x_{0},x_{1},y_{0},y_{1}$) to be storable in lookup tables for a dense
enough grid of values of variables. Therefore, we need to simplify the
integrals by exploiting all the correlations between the variables.

First, let us write $I$ as a linear combination of at most 4 integrals, called
$I_{0}$, which contain the same integrand as $I$, but defined on rectangles
with $(0,0)$ as one of their corners. In Figure 1 we can see 5 different
possible cases of the initial rectangle with respect to the source point
$(0,0)$. In (a) and (b), the initial rectangle $D$ is split into 4 and 2
smaller rectangles, respectively. Case (c) is a special case, when $I=I_{0}$
and no treatment is needed. In (d), the source point is outside the initial
rectangle but one of rectangle’s edge coordinate is 0 – rectangle $D$ is
extended to $D_{1}$ and then $D_{2}$ is subtracted from the augmented one. In
(e), $D_{2}$ and $D_{3}$ are subtracted from $D_{1}$, and their intersection
$D_{4}$ is added back.

$D_{3}$ $D_{4}$ $D_{1}$ $D_{2}$ (a) Interior source point:
“$D=D_{1}+D_{2}+D_{3}+D_{4}$”

$D_{1}$ $D_{2}$ (b) Edge source point:
“$D=D_{1}+D_{2}$” (c) Corner source point

$D_{2}$ $D_{1}$ (d) Aligned exterior source point:
“$D=D_{1}-D_{2}$”

$D_{3}$ $D_{2}$ $D_{4}$ $D_{1}$ (e) General exterior source point:
“$D=D_{1}-D_{2}-D_{3}+D_{4}$”

Figure 1: All possible positions of the rectangle $D$ (in gray) with respect
to the source point $(0,0)$, depicted as a gray dot. Sub-rectangles are shown
in different colours and for better visualization the overlapping ones are
depicted at different vertical layers.

For each $I_{0}$, defined on a rectangle $[0,x_{2}]\times[0,y_{2}]$, we apply
a linear transformation to map the rectangle to the unit square $[0,1]^{2}$.
Then we divide $R$ by ${a}^{p/2}$. Therefore

$\displaystyle
I_{0}=x_{2}^{p+q+1}y_{2}^{r+1}a^{p/2}\int_{0}^{1}\int_{0}^{1}R\left(1,\frac{b}{a}\frac{y_{2}}{x_{2}},\frac{c}{a}\frac{y_{2}^{2}}{x_{2}^{2}},(x,y)\right)^{p}x^{q}y^{r}\mathop{dx}\mathop{dy}.$

Observe that the integral in $I_{0}$ depends only on two continuous variables
$\bar{b}:=({b}/{a})({y_{2}}/{x_{2}})$,
$\bar{c}:=({c}/{a})({y_{2}}/{x_{2}})^{2}$ and on integers $p,q,r$. For every
needed combination of integers $p,q,r$ we can beforehand compute integrals
$\int_{0}^{1}\int_{0}^{1}R\left(1,\bar{b},\bar{c},(x,y)\right)^{p}x^{q}y^{r}\mathop{dx}\mathop{dy}$,
using the formulae in Section A.2, for a dense enough grid of points
$(\bar{b},\bar{c})$ and store them in a lookup table. Intermediate values, not
stored in the table, can be computed via interpolation.

####  2.4.2 Triangles

A simplification of the integration on a general triangle $T$ with vertices
$(x_{0},y_{0}),(x_{1},y_{1}),(x_{2},y_{2})$ is more involved compared to the
previous aligned rectangular domains. Here we demonstrate one way how to
proceed for a general case as a proof of concept – clearly it is not the
optimal nor the most numerical stable approach for every $T$. Since integral

$\displaystyle I:=\int_{T}R^{p}x^{q}y^{r}\,dT$

depends on 9 “continuous” variables
($a,b,c,x_{0},x_{1},x_{2},y_{0},y_{1},y_{2}$), we can write $I$ again as
linear combination of at most 3 integrals, denoted by $I_{0}$. These integrals
contain the same integrand as $I$, but are defined on triangles with $(0,0)$
as one of their corners. The cases are shown in Figure 2. In (a) and (b), the
initial rectangle $T$ is split into 3 and 2 smaller triangles, respectively.
Case (c) is again a special case, when $I=I_{0}$ and no treatment is needed.
In (d), the source point is outside of $T$ and close to triangle’s edge –
triangle $T$ is a union of $T_{1}$ and $T_{2}$, where we subtract triangle
$T_{3}$. In (e), $T$ is extended to $T_{1}$ and then $T_{2}$ and $T_{3}$ are
subtracted.

$T_{2}$ $T_{3}$ $T_{1}$ (a) Interior source point:
“$T=T_{1}+T_{2}+T_{3}$”

$T_{2}$ $T_{1}$ (b) Edge source point:
“$T=T_{1}+T_{2}$” (c) Corner source point

$T_{2}$ $T_{3}$ $T_{1}$ (d) Exterior source point, close to an edge:
“$T=T_{1}+T_{2}-T_{3}$”

$T_{3}$ $T_{2}$ $T_{1}$ (e) Exterior source point, close to a corner:
“$T=T_{1}-T_{2}-T_{3}$”

Figure 2: All possible positions of a triangle $T$ (in gray) with respect to
the source point $(0,0)$, depicted as a black dot. Sub-triangles are shown in
different colours and for better visualization the overlapping ones are
depicted at different vertical layers.

In the next step we focus on each of those newly constructed triangles. Let
$T_{0}$ be the integration domain of $I_{0}$, hence a triangle with vertices
$(0,0),(x_{1},y_{1}),(x_{2},y_{2})$. Linear transformation that maps a
triangle with vertices $(0,0),(1,0),(0,1)$ to $T_{0}$ is represented by matrix

$\displaystyle\begin{bmatrix}x_{1}&x_{2}\\\ y_{1}&y_{2}\end{bmatrix}.$

Therefore, the integral $I_{0}$ can be written as

$\displaystyle
I_{0}:=(x_{1}y_{2}-x_{2}y_{1})\int_{0}^{1}\int_{0}^{1-y}R(\bar{a},\bar{b},\bar{c},(x,y))^{p}(x_{1}x+x_{2}y)^{q}(y_{1}x+y_{2}y)^{r}\mathop{dx}\mathop{dy},$
(26)

where

$\displaystyle\bar{a}$ $\displaystyle:=ax_{1}^{2}+bx_{1}y_{1}+cy_{1}^{2},$
$\displaystyle\bar{b}$
$\displaystyle:=2ax_{1}x_{2}+bx_{2}y_{1}+bx_{1}y_{2}+2cy_{1}y_{2},$
$\displaystyle\bar{c}$ $\displaystyle:=ax_{2}^{2}+bx_{2}y_{2}+cy_{2}^{2}.$

The polynomial part in (26) can be written as

$\displaystyle(x_{1}x+x_{2}y)^{q}(y_{1}x+y_{2}y)^{r}=\sum_{i=0}^{q+r}\alpha_{i}x^{q+r-i}y^{i}$
(27)

for suitable coefficients $\alpha_{i}$. By exchanging the order of integration
and summation and dividing $R$ by ${a}^{p/2}$, $I_{0}$ can be written as
linear combination of functions

$\displaystyle\int_{0}^{1}\int_{y}^{1}R\left(1,\frac{\bar{b}}{\bar{a}},\frac{\bar{c}}{\bar{a}},(x,y)\right)^{p}x^{q+r-i}y^{i}\mathop{dx}\mathop{dy}.$

The latter integral depends only on two continuous variables
$\bar{\bar{b}}:={\bar{b}}/{\bar{a}},\bar{\bar{c}}:={\bar{c}}/{\bar{a}}$ and
integer powers $p,q+r-i,i$, and we can evaluate it for a dense enough grid of
points $(\bar{\bar{b}},\bar{\bar{c}})$ and store the values in a lookup table.
At first glance it might seem that every $I_{0}$ in (26) needs to be computed
as a sum of $q+r+1$ integrals due to formula (27). However, observe that the
polynomial $P$ in (25) already contains different powers of monomials and thus
we can group the integrals with the same powers of $x$ and $y$ to reduce the
amount of computations.

###  2.5 Smoothness of the regularized kernels

The smoothness of the regularized kernel $\rho$ at ${\boldsymbol{s}}$
increases by one if $n$ is increased by one. First we consider the more common
singularity subtraction technique, where the analysis is very straightforward.

######  Theorem 3.

For a singular kernel ${\mathcal{K}}_{\boldsymbol{s}}$, defined in (20), let
$\rho={\mathcal{K}}_{\boldsymbol{s}}-{\mathcal{K}}_{{\boldsymbol{s}},n}$ be
the regularized kernel with $m+n\geq 1$. Then $\rho$ is a $C^{m+n-1}$ smooth
function at ${\boldsymbol{z}}=0$ and $D^{\alpha}\rho=0$ for $0\leq|\alpha|\leq
m+n-1$.

Proof. Let us write
${\mathcal{K}}_{\boldsymbol{s}}={\mathcal{K}}_{{\boldsymbol{s}},n}+\tau_{n}$,
where $\tau_{n}$ is the corresponding tail in the series expansion of
${\mathcal{K}}_{\boldsymbol{s}}$. Then $\rho=\tau_{n}$ and from (20) it is
clear that $\zeta(\rho)\geq m+n$. Apply Proposition 12 and 13 and the proof is
complete.

∎

For the regularization by the division step we need to additionally assume
$P_{m_{2}+3}^{[1]}({\boldsymbol{z}})\neq 0$ for ${\boldsymbol{z}}\in
D_{\varepsilon}$. This condition is needed to not introduce new singularities
in the regularized kernel $\rho$. In the next proposition we locally assume
$P_{m_{2}+3}^{[1]}>0$ but an analogous statement can be done for
$P_{m_{2}+3}^{[1]}<0$.

######  Corollary 4.

Let $P_{m_{2}+3}^{[1]}>0$ for ${\boldsymbol{z}}\in D_{\varepsilon}$ for some
$\varepsilon>0$. Then there exists $\varepsilon^{\prime}>0$ such that
${\mathcal{K}}_{{\boldsymbol{s}},n}({\boldsymbol{z}})>0$ for every
${\boldsymbol{z}}\in D_{\varepsilon^{\prime}}$.

Proof. Write

$\displaystyle{\mathcal{K}}_{{\boldsymbol{s}},n}=\sum_{\ell=1}^{n}R^{-2(\ell+m_{1})+1}P_{3\ell+m_{2}}^{[\ell]}=R^{-2(n+m_{1})+1}\left(R^{2n-2}P_{m_{2}+3}^{[1]}+\sum_{\ell=2}^{n}R^{2(n-\ell)}P_{3\ell+m_{2}}^{[\ell]}\right).$
(28)

Function $R^{-2(n+m_{1})+1}$ is positive, thus it is enough to analyze the
remaining factor. Homogeneous polynomial $R^{2n-2}P_{m_{2}+3}^{[1]}$ is of
degree $2n+m_{2}+1$ and positive for ${\boldsymbol{z}}\in D_{\varepsilon}$.
Each polynomial in the remaining sum is of degree greater than $2n+m_{2}+1$,
thus it decays to zero faster than $R^{2n-2}P_{m_{2}+3}^{[1]}$, when
${\boldsymbol{z}}$ goes to $\boldsymbol{0}$. We can deduce that the kernel
${\mathcal{K}}_{{\boldsymbol{s}},n}$ is positive in $D_{\varepsilon^{\prime}}$
for sufficiently small $\varepsilon^{\prime}>0$.

∎

######  Remark 5.

The additional condition $P_{m_{2}+3}^{[1]}\neq 0$ is clearly an notable
drawback of the regularization via division, since not all kernels are
necessary of the same sign near the source point ${\boldsymbol{s}}$, e.g.,
kernel ${\mathcal{H}}$. For these types of kernel the division approach cannot
be directly applied in this simple form.

The case when ${\mathcal{K}}_{\boldsymbol{s}}({\boldsymbol{z}})\neq 0$ for all
${\boldsymbol{z}}\in D_{\varepsilon}$ but
${\mathcal{K}}_{{\boldsymbol{s}},n}({\boldsymbol{z}})=0$ for some
${\boldsymbol{z}}\in D_{\varepsilon}$ is far less problematic. For example, by
adding a simple correction term $\eta\|{\boldsymbol{z}}\|_{2}^{2}$ to
${\mathcal{K}}_{{\boldsymbol{s}},n}({\boldsymbol{z}})$ for a suitable $\eta$,
we can enforce function
${\mathcal{K}}_{{\boldsymbol{s}},n}({\boldsymbol{z}})+\eta\|{\boldsymbol{z}}\|_{2}^{2}$
to not vanish inside $D_{\varepsilon}$. Namely, $\eta>0$ if
${\mathcal{K}}_{\boldsymbol{s}}>0$, and $\eta<0$ if
${\mathcal{K}}_{\boldsymbol{s}}<0$.

Despite the mentioned drawbacks, probably the main advantage of the division
technique is the improved smoothness of $\rho$ when compared to the
subtraction splitting for the same amount of regularization terms.
Furthermore, the gap between the approaches is even bigger for kernels
${\mathcal{K}}_{\boldsymbol{s}}$ with stronger singularity. Indeed, for the
same number of terms $n$ the smoothness of $\rho$ is the same, regardless of
the type of singularity, when the regularization via division step is used.

The following lemma is needed before proving the theorem for the smoothness of
$\rho$.

######  Lemma 6.

Let $\varepsilon>0$ and let

$\displaystyle
p_{m}({\boldsymbol{z}}):=\sum_{|\alpha|=m}c_{\alpha}{\boldsymbol{z}}^{\alpha}$

be a bivariate homogeneous polynomial of degree $m$, such that
$p_{m}({\boldsymbol{z}})>0$ for $\|{\boldsymbol{z}}\|_{2}<\varepsilon$. Then
$1/({R^{\ell}p_{m}})$ is a continuous function at
${\boldsymbol{z}}=\boldsymbol{0}$ if $\zeta({R^{\ell}p_{m}})=\ell+m\leq-1$.

Proof. The proof is similar to the proof of Proposition 12 but considering a
negative parameter $\zeta$ instead of a positive one.

Fix a sufficiently small $\delta>0$. Let us show that there exists a small
enough $\varepsilon$, $0<\varepsilon<1$, such that for
$\|{\boldsymbol{z}}\|_{2}\leq\varepsilon$ it follows
$1/({R({\boldsymbol{z}})^{\ell}p_{m}({\boldsymbol{z}})})\leq\delta$.

Let us again consider only the case ${\boldsymbol{z}}\in T_{1}$ (see (38)),
since the proof for the other subdomains follows similarly. We can write
$z_{2}=\gamma z_{1}$ for $0\leq\gamma\leq 1$. From Lemma 10 we can infer that
there exists a positive constant $C_{abc\ell,1}$ such that
$C_{abc\ell,1}\leq({a+b\gamma+c\gamma^{2}})^{\ell/2}$ for all $\gamma$. From
the assumptions we know that
$p_{m}({\boldsymbol{z}})=z_{1}^{m}\sum_{|\alpha|=m}c_{\alpha}\gamma^{\alpha_{2}}>0$
for $z_{1}>0$, thus there exists a positive constant $C$ such that
$C\leq\sum_{|\alpha|=m}c_{\alpha}\gamma^{\alpha_{2}}$. By combining all the
derived inequalities we can show that

$\displaystyle
R({\boldsymbol{z}})^{\ell}p_{m}({\boldsymbol{z}})={(a+b\gamma+c\gamma^{2})}^{\ell/2}z_{1}^{\ell+m}\sum_{|\alpha|=m}c_{\alpha}\gamma^{\alpha_{2}}\geq
C_{abc\ell,1}\,C\,z_{1}^{\ell+m}$

and

$\displaystyle\frac{1}{R({\boldsymbol{z}})^{\ell}p_{m}({\boldsymbol{z}})}\leq
C_{abc\ell,1}^{-1}\,C^{-1}\,z_{1}\leq
C_{abc\ell,1}^{-1}\,C^{-1}\,\varepsilon.$

By setting $\varepsilon\leq C_{abc\ell,1}\,C\,\delta$ the proof is complete.

∎

######  Theorem 7.

Let ${\mathcal{K}}_{{\boldsymbol{s}},n}$ be an approximation of
${\mathcal{K}}_{\boldsymbol{s}}$, defined in (20), with
$P_{m_{2}+3}^{[1]}({\boldsymbol{z}})\neq 0$ for ${\boldsymbol{z}}\in
D_{\varepsilon}$. Let
$\rho={\mathcal{K}}_{\boldsymbol{s}}/{\mathcal{K}}_{{\boldsymbol{s}},n}$ be
the regularized kernel. Then $\rho$ is a $C^{n-1}$ smooth function at
${\boldsymbol{z}}=\boldsymbol{0}$ and $D^{\alpha}\rho(\boldsymbol{0})=0$ for
$1\leq|\alpha|\leq n-1$.

Proof. Let us write again
${\mathcal{K}}_{\boldsymbol{s}}={\mathcal{K}}_{{\boldsymbol{s}},n}+\tau_{n}$,
where $\tau_{n}$ is the corresponding tail in the series expansion of
${\mathcal{K}}_{\boldsymbol{s}}$.

First, let us prove it for $n=1$, i.e., function
$\rho={\mathcal{K}}_{\boldsymbol{s}}/{\mathcal{K}}_{{\boldsymbol{s}},1}$ is
$C^{0}$ continuous. The function $\rho$ can be written as

$\displaystyle\rho=\frac{{\mathcal{K}}_{\boldsymbol{s}}}{{\mathcal{K}}_{{\boldsymbol{s}},1}}=1+\frac{\tau_{1}}{{\mathcal{K}}_{{\boldsymbol{s}},1}}=1+\frac{R^{-m-1}\tau_{1}}{R^{-m-1}{\mathcal{K}}_{{\boldsymbol{s}},1}}=1+\frac{R^{-m-1}\tau_{1}}{R^{-m-2m_{1}-2}P_{m_{2}+3}^{[1]}}.$

Since $\zeta(R^{-m-1}\tau_{1})\geq(-m-1)+m+1=0$ function
$R^{-m-1}\tau_{1}({\boldsymbol{z}})$ is bounded in $D_{\varepsilon}$ (see
Proposition 9) and we can find a positive constant $C_{1}$, independent of
${\boldsymbol{z}}$, such that
$|R({\boldsymbol{z}})^{-m-1}\tau_{1}({\boldsymbol{z}})|\leq C_{1}$. For the
denominator it holds $\zeta(R^{-m-1}{\mathcal{K}}_{{\boldsymbol{s}},1})=-1$,
thus it is not bounded for $\|{\boldsymbol{z}}\|_{2}\leq\varepsilon$. By
defining $\rho(\boldsymbol{0})=1$ we can see that

$\displaystyle|\rho({\boldsymbol{z}})-\rho(\boldsymbol{0})|=|\rho({\boldsymbol{z}})-1|=\left|\frac{R({\boldsymbol{z}})^{-m-1}\tau_{1}({\boldsymbol{z}})}{R({\boldsymbol{z}})^{-m-1}{\mathcal{K}}_{{\boldsymbol{s}},1}({\boldsymbol{z}})}\right|\leq\frac{C_{1}}{|R({\boldsymbol{z}})^{-m-1}{\mathcal{K}}_{{\boldsymbol{s}},1}({\boldsymbol{z}})|}.$

From Lemma 6 it follows that the function
$1/|R^{-m-1}{\mathcal{K}}_{{\boldsymbol{s}},1}|$ is continuous at
${\boldsymbol{z}}=\boldsymbol{0}$, hence so is $\rho$.

Let us prove the remaining part by induction. Assume $D^{\alpha}\rho$ is
continuous and $D^{\alpha}\rho(\boldsymbol{0})=0$ for $0\leq|\alpha|\leq n-1$.
Then we can define the following derivatives via recursion,

$\displaystyle\rho=1+\frac{\tau_{n}}{{\mathcal{K}}_{{\boldsymbol{s}},{n}}}=:1+\frac{p_{0,0}}{q_{0,0}},$
$\displaystyle
D^{i+1}_{1}D^{j}_{2}\rho=D^{1}_{1}\frac{p_{i,j}}{q_{i,j}}=\frac{D^{1}_{1}p_{i,j}\,q_{i,j}-p_{i,j}\,D^{1}_{1}q_{i,j}}{(q_{i,j})^{2}}=:\frac{p_{i+1,j}}{q_{i+1,j}},$
(29) $\displaystyle
D^{i}_{1}D^{j+1}_{2}\rho=D^{1}_{2}\frac{p_{i,j}}{q_{i,j}}=\frac{D^{1}_{2}p_{i,j}\,q_{i,j}-p_{i,j}\,D^{2}_{1}q_{i,j}}{(q_{i,j})^{2}}=:\frac{p_{i,j+1}}{q_{i,j+1}},$

for $i+j\leq n-1$. Observe that $\zeta(p_{0,0})=\zeta(\tau_{n})\geq m+n$,
while $\zeta(q_{0,0})=\zeta({\mathcal{K}}_{{\boldsymbol{s}},{n}})=m$. From
expressions (2.5) and Proposition 9 we can infer that
$\zeta(q_{i,j})=2^{i+j}m$ and $\zeta(p_{i,j})\geq 2^{i+j}m+n-i-j$.

For $|\alpha|=i+j=n-1$ we can write

$\displaystyle\left|\frac{p_{i,j}({\boldsymbol{z}})}{q_{i,j}({\boldsymbol{z}})}-\boldsymbol{0}\right|=\frac{|R({\boldsymbol{z}})^{-2^{n-1}m-1}\,p_{i,j}({\boldsymbol{z}})|}{R({\boldsymbol{z}})^{-2^{n-1}m-1}\,|q_{i,j}({\boldsymbol{z}})|}\leq\frac{C_{2}}{R({\boldsymbol{z}})^{-2^{n-1}m-1}\,|q_{i,j}({\boldsymbol{z}})|},$
(30)

where $C_{2}$ is a positive constant, since
$\zeta(R^{-2^{n-1}m-1}\,p_{i,j})\geq 0$. From the proof of Corollary 4 it is
easy to see that there exists a small enough $\varepsilon^{\prime}>0$ such
that

$\displaystyle|q_{i,j}({\boldsymbol{z}})|=|{\mathcal{K}}_{{\boldsymbol{s}},{n}}^{2^{n-1}}({\boldsymbol{z}})|\geq\left|\left(\frac{1}{2}R({\boldsymbol{z}})^{-2m_{1}-1}P_{m_{2}+3}^{[1]}({\boldsymbol{z}})\right)^{2^{n-1}}\right|$

for all ${\boldsymbol{z}}\in D_{\varepsilon^{\prime}}$. Then for
${\boldsymbol{z}}\in D_{\varepsilon^{\prime}}$ we can further estimate the
right-hand side of (30),

$\displaystyle\frac{C_{2}}{R({\boldsymbol{z}})^{-2^{n-1}m-1}\,|q_{i,j}({\boldsymbol{z}})|}\leq\frac{2^{2^{n-1}}C_{2}}{R({\boldsymbol{z}})^{-2^{n-1}(m_{2}+3)-1}\,|P_{m_{2}+3}^{[1]}({\boldsymbol{z}})^{2^{n-1}}|}.$
(31)

The function on the right-hand side of (31) is continuous at
${\boldsymbol{z}}=\boldsymbol{0}$ by Lemma 6. Hence,
$D^{\alpha}\rho={p_{i,j}}/{q_{i,j}}$ is continuous for $|\alpha|=i+j=n-1$ and
$D^{\alpha}\rho(\boldsymbol{0})=0$.

∎

##  3 Numerical tests

In all experiments we test the accuracy of numerical integration for specific
singular integrals with respect to the number of terms $n$ in the singular
kernel series expansion; we vary $n=0,1,2,3$, whereas in the case $n=0$ we do
not apply any regularization of the kernel. For regularization we considered
the singularity subtraction and division (see Section 2.3).

For geometry let us consider a section of a spheroid, parameterized by quartic
tensor product NURBS ${\boldsymbol{F}}:[0,1]^{2}\to{\mathop{{\rm
I}\kern-1.99997pt{\rm R}}\nolimits}^{3}$ (obtained by stretching a spherical
section [10]); see Fig. 3.

Figure 3: Control points and geometry of the sphereoid.

Its control points are
$\boldsymbol{c}_{ij}=(C^{x}_{ij},C^{y}_{ij},C^{z}_{ij})/C^{w}_{ij}$ by
defining the following matrices

$\displaystyle C^{x}$
$\displaystyle=1.5\begin{bmatrix}-c_{1}&-c_{3}&-c_{4}&-c_{3}&-c_{1}\\\
-c_{2}&-c_{5}&-c_{6}&-c_{5}&-c_{2}\\\ 0&0&0&0&0\\\
c_{2}&c_{5}&c_{6}&c_{5}&c_{2}\\\
c_{1}&c_{3}&c_{4}&c_{3}&c_{1}\end{bmatrix},\;\;C^{y}=-2\begin{bmatrix}c_{1}&c_{3}&c_{4}&c_{3}&c_{1}\\\
c_{3}&c_{7}&c_{8}&c_{7}&c_{3}\\\ c_{4}&c_{8}&c_{9}&c_{8}&c_{4}\\\
c_{3}&c_{7}&c_{8}&c_{7}&c_{3}\\\
c_{1}&c_{3}&c_{4}&c_{3}&c_{1}\end{bmatrix},\;\;C^{z}=\begin{bmatrix}-c_{1}&-c_{2}&0&c_{2}&c_{1}\\\
-c_{3}&-c_{5}&0&c_{5}&c_{3}\\\ -c_{4}&-c_{6}&0&c_{6}&c_{4}\\\
-c_{3}&-c_{5}&0&c_{5}&c_{3}\\\ -c_{1}&-c_{2}&0&c_{2}&c_{1}\end{bmatrix},$
$\displaystyle C_{w}$
$\displaystyle=\begin{bmatrix}w_{1}&w_{2}&c_{9}&w_{2}&w_{1}\\\
w_{2}&c_{7}&w_{3}&c_{7}&w_{2}\\\ c_{9}&w_{3}&w_{4}&w_{3}&c_{9}\\\
w_{2}&c_{7}&w_{3}&c_{7}&w_{2}\\\ w_{1}&w_{2}&c_{9}&w_{2}&w_{1}\end{bmatrix},$

where

$\displaystyle c_{1}$ $\displaystyle=4(\sqrt{3}-1),$ $\displaystyle c_{2}$
$\displaystyle=\sqrt{2},$ $\displaystyle c_{3}$
$\displaystyle=\sqrt{2}(4-\sqrt{3}),$ $\displaystyle c_{4}$
$\displaystyle=4(2\sqrt{3}-1)/3,$ $\displaystyle c_{5}$
$\displaystyle=(3\sqrt{3}-2)/2,$ $\displaystyle c_{6}$
$\displaystyle=\sqrt{2}(7-2\sqrt{3})/3,$ $\displaystyle c_{7}$
$\displaystyle=(\sqrt{3}+6)/2,$ $\displaystyle c_{8}$
$\displaystyle=5\sqrt{6}/3,$ $\displaystyle c_{9}$
$\displaystyle=4(5-\sqrt{3})/3,$ $\displaystyle w_{1}$
$\displaystyle=4(3-\sqrt{3}),$ $\displaystyle w_{2}$
$\displaystyle=\sqrt{2}(3\sqrt{3}-2),$ $\displaystyle w_{3}$
$\displaystyle=\sqrt{2}(\sqrt{3}+6)/3,$ $\displaystyle w_{4}$
$\displaystyle=4(5\sqrt{3}-1)/9.$

###  3.1 N-refinement

In the first example we consider a simplified version of a governing integral
in BIE (4) to test the effect of the regularization technique in the most
direct manner. We test the accuracy of the numerical schemes in Section 2.3 to
compute

$\displaystyle\frac{1}{4\pi}\int_{[0,1]^{2}}{\mathcal{K}}({\boldsymbol{s}},{\boldsymbol{t}})\mathop{d\boldsymbol{t}},$
(32)

where ${\mathcal{K}}\in\\{{\mathcal{G}},\bar{\mathcal{H}}\\}$ for 16 different
source points ${\boldsymbol{s}}=(s_{1},s_{2})$ with
$s_{1},s_{2}\in\\{0.6,0.7,0.8,0.9\\}$ and for different number of quadrature
points $N$. The integral can be considered as a simplified version of (22),
where $B_{j}\equiv 1$ and $v\equiv 1$.

In Fig. 4 and Fig. 5 we show the plots of the regularized kernels $\rho$ for
the subtraction and the division regularization step via kernel
${\mathcal{K}}_{{\boldsymbol{s}},n}$ for ${\boldsymbol{s}}=(0.6,0.6)$. As
analyzed in Section 2.5, for the singularity subtraction the function $\rho$
is bounded, $C^{0}$ and $C^{1}$ for $n=1,2,3$, respectively. Instead, for the
division step the regularity of $\rho$ increases by one. For the latter
regularization for $n=3$ the function $\rho$ exhibits singularities near the
edges of the domain $[0,1]^{2}$; the issue is resolved by adding the
correction term with $\eta=10$ (see Remark 5).

After the subtraction of the singularity in (32), the singular part of the
integral is computed analytically, using the recursive formulae in Section
A.2. The regular part of the integral is computed numerically, using the
tensor product Gauss-Legendre quadrature. For each iteration step we double
the number of quadrature nodes $N$ in each direction. Since the initial $N=10$
is already sufficiently high, the accuracy of the integration is expected to
be governed by the lower smoothness of the integrand at ${\boldsymbol{s}}$.
The numerical values are compared against the values obtained by highly
accurate numerical scheme that combines the Duffy transformation [15] and
Matlab’s adaptive quadrature routines. The regularization via division,
however, cannot be used in this test in such a simple manner without a further
approximation of the integrand and is thus excluded from the test. In Fig. 6
by examining the error plots with respect to $N$ we can clearly see the
benefits of using higher $n$ in the regularization. With increasing $n$, the
accuracy and the convergence order increase for both type of kernels; for
$n=1,2,3$ the error decays as $CN^{-(n+1)}$ for some constant $C>0$.

(a) Subtraction, $n=1$ (b) Subtraction, $n=2$ (c) Subtraction, $n=3$ (d)
Division, $n=1$

(e) Division, $n=2$ (f) Division, $n=3$

Figure 4: Regularized kernel $\rho$ for ${\mathcal{G}}_{\boldsymbol{s}}$ with
${\boldsymbol{s}}=(0.6,0.6)$ marked as a black dot.

(a) Subtraction, $n=1$ (b) Subtraction, $n=2$ (c) Subtraction, $n=3$ (d)
Division, $n=1$

(e) Division, $n=2$ (f) Division, $n=3$

Figure 5: Regularized kernel $\rho$ for $\bar{\mathcal{H}}_{\boldsymbol{s}}$
with ${\boldsymbol{s}}=(0.6,0.6)$ marked as a black dot.

###  3.2 h-refinement

In this experiment we perform an h-refinement for the discretization space; we
vary size $h$ of the support of basis functions $B$ in the integrals

$\displaystyle\frac{1}{4\pi}\int_{[0,1]^{2}}{\mathcal{K}}({\boldsymbol{s}},{\boldsymbol{t}})B({\boldsymbol{t}})v({\boldsymbol{t}})\mathop{d\boldsymbol{t}},$

where again ${\mathcal{K}}\in\\{{\mathcal{G}},\bar{\mathcal{H}}\\}$. To
analyze the integral of BIE with a constant boundary datum, let us fix $v=J$,
when we test kernel ${\mathcal{G}}$, and let $v\equiv 1$ for kernel
$\bar{\mathcal{H}}$ since $J$ is already included in the definition of the
kernel. Let us fix 16 source points ${\boldsymbol{s}}=(s_{1},s_{2})$ with
$s_{1},s_{2}\in\\{0.6,0.7,0.8,0.9\\}$. Let $B$ be a piece-wise constant basis
function. To analyze only truly singular integrals, let us for each
${\boldsymbol{s}}$ consider “active” only the following 5 basis functions $B$,
denoted as

$\displaystyle
B_{\boldsymbol{s}}^{(i)}({\boldsymbol{t}}):=\left\\{\begin{array}[]{ll}1,&{\boldsymbol{t}}\in
D_{\boldsymbol{s}}^{(i)}\\\ 0,&{\boldsymbol{t}}\notin
D_{\boldsymbol{s}}^{(i)}\end{array}\right.,$

with supports on squares of sizes $h\times h$,

$\displaystyle D_{\boldsymbol{s}}^{(1)}$
$\displaystyle:=[s_{1}-h/2,s_{1}+h/2]\times[s_{2}-h/2,s_{2}+h/2],$
$\displaystyle D_{\boldsymbol{s}}^{(2)}$
$\displaystyle:=[s_{1}-2h/3,s_{1}+h/3]\times[s_{2}-2h/3,s_{2}+h/3],$
$\displaystyle D_{\boldsymbol{s}}^{(3)}$
$\displaystyle:=[s_{1}-2h/3,s_{1}+h/3]\times[s_{2}-h/3,s_{2}+2h/3],$
$\displaystyle D_{\boldsymbol{s}}^{(4)}$
$\displaystyle:=[s_{1}-h/3,s_{1}+2h/3]\times[s_{2}-2h/3,s_{2}+h/3],$
$\displaystyle D_{\boldsymbol{s}}^{(5)}$
$\displaystyle:=[s_{1}-h/3,s_{1}+2h/3]\times[s_{2}-h/3,s_{2}+2h/3],$

that are inside the domain $[0,1]^{2}$.

Following the itinerary in Section 2.3, by applying the subtraction of
singularity, we derive to formula (23). The first integral is computed with
Gauss-Legendre quadrature on domain $D_{\boldsymbol{s}}^{(i)}$ with $10\times
10$ quadrature nodes. In the second one, the regular part of the integrand,
i.e., $J$ is approximated with a polynomial of bi-degree $(4,4)$ using a
least-squares method for the same $10\times 10$ evaluation sites. The second
integral of (23) is thus approximated with an integral (25) and we integrate
it analytically.

Similarly, for the singularity division, we transform the integral into (24).
The regular part $\rho$ of the integrand is approximated with a polynomial
using the least-squares approach for the same $10\times 10$ evaluation sites.
The obtained integral is of form (25) and we integrate it analytically.

As in the first experiment, we expect the accuracy of integration schemes to
be subjected to reduced smoothness of $\rho$ at ${\boldsymbol{s}}$. In Fig. 7
we can see the maximum error plots for all ${\boldsymbol{s}}$ and
$D_{\boldsymbol{s}}^{(i)}$ for the numerical evaluation of integrals with
respect to the support size $h$. Again, there is a clear benefit in using
higher $n$ for both singularity subtraction and division. Both regularization
techniques produce similar accuracy for both kernels and for all tested $h$.

(a) kernel ${\mathcal{G}}_{\boldsymbol{s}}$ (b) kernel
$\bar{\mathcal{H}}_{\boldsymbol{s}}$

Figure 6: Maximum error plots for numerical integration involving kernels
${\mathcal{G}}_{\boldsymbol{s}}$ and $\bar{\mathcal{H}}_{\boldsymbol{s}}$
using singularity subtraction and $N$-refinement.

(a) kernel ${\mathcal{G}}_{\boldsymbol{s}}$ (b) kernel
$\bar{\mathcal{H}}_{\boldsymbol{s}}$

Figure 7: Maximum error plots for numerical integration involving kernels
${\mathcal{G}}_{\boldsymbol{s}}$ and $\bar{\mathcal{H}}_{\boldsymbol{s}}$
using singularity subtraction and division, and $h$-refinement. Plots with
solid and dashed lines correspond to results using the subtraction and
division procedure, respectively.

## Conclusion

The singularity extraction is studied for singular integrals in BIE for 3D
potential problems on smooth geometries, focusing on the single and and double
layer operator. Integrals of truncated series expansions of singular kernels
are computed via recurrence formulae – to reduce the complexity of this
process, a speedup via lookup tables that exploits correlated variables is
provided.

Numerical tests demonstrate that the presented extraction technique can be a
useful prerequisite for a numerical quadrature to recover the optimal order of
convergence of the approximate solution in BEM with a small number of
quadrature nodes. Since the extraction acts directly on the starting
parametric domain, no modification of integration domain is needed – this is
an advantageous property since the same quadrature nodes can be reused for
several integrals involving neighbouring trial basis functions. The truncation
error of the numerical integration rule can be attributed solely to the
truncation error of approximating the regular part of the integrand; thus the
accuracy of the rule can be effectively controlled by the number of terms in
the series extraction and by the choice of the approximation operator for
regular functions. Singularity subtraction and division give comparable
results; the technique via division is arguably slightly more efficient since
the governing integral does not need to be split into two integrals, however
it cannot directly handle kernels that change signs.

Interesting topics for future work include a study of other kernels with
stronger types of singularities and of non-stationary nature. Efficient
implementation for computing singular integrals via lookup tables is beyond
the scope of this paper, but it should be investigated in future work. Finding
better ways to approximate kernels across non-smooth interfaces of the
geometry, for examples on multipatch domains, could improve accuracy of the
integration routines on more complex geometries.

## References

## References

  * Aimi et al. [2018] Aimi, A., Calabrò, F., Diligenti, M., Sampoli, M. L., Sangalli, G., Sestini, A., 2018. Efficient assembly based on B-spline tailored quadrature rules for the IgA-SGBEM. Comput. Methods Appl. Mech. Engrg. 331, 327–342.
  * Aimi et al. [2020] Aimi, A., Calabrò, F., Falini, A., Sampoli, M. L., Sestini, A., 2020. Quadrature formulas based on spline quasi-interpolation for hypersingular integrals arising in IgA-SGBEM. Comput. Methods Appl. Mech. Engrg. 372, 113441\.
  * Aimi and Diligenti [2002] Aimi, A., Diligenti, M., 2002. Numerical integration in 3D Galerkin BEM solution of HBIEs. Comput. Mech. 28, 233–249.
  * Atkinson [2009] Atkinson, K. E., 2009. The Numerical Solution of Integral Equations of the Second Kind. Cambridge University Press.
  * Beer et al. [2015] Beer, G., Marussig, B., Zechner, J., 2015. A simple approach to the numerical simulation with trimmed CAD surfaces. Comput. Methods Appl. Mech. Engrg. 285, 776–790.
  * Botha [2013] Botha, M. M., 2013. A family of augmented Duffy transformations for near-singularity cancellation quadrature. IEEE Antennas Wirel. Propag. Lett. 61 (6), 3123–3134.
  * Calabrò et al. [2018] Calabrò, F., Falini, A., Sampoli, M. L., Sestini, A., 2018. Efficient quadrature rules based on spline quasi-interpolation for application to IgA-BEMs. J. Comput. Appl. Math. 338, 153–167.
  * Calabrò et al. [2019] Calabrò, F., Loli, G., Sangalli, G., Tani, M., 2019. Quadrature rules in the isogeometric Galerkin method: State of the art and an introduction to weighted quadrature. In: Giannelli, C., Speleers, H. (Eds.), Advanced Methods for Geometric Modeling and Numerical Simulation. Vol. 35 of Springer INdAM Series. Springer Cham, pp. 43–55.
  * Chen and Hong [1999] Chen, J. T., Hong, H.-K., 1999. Review of dual boundary element methods with emphasis on hypersingular integrals and divergent series. Appl. Mech. Rev. 52 (1), 17–33.
  * Cobb [1988] Cobb, J. E., 1988. Tiling the sphere with rational Bézier patches. Technical report UUCS-88-009, Computer Science, University of Utah.
  * Costabel [1986] Costabel, M., 1986. Principles of boundary element methods. Techn. Hochsch., Fachbereich Mathematik.
  * Cottrell et al. [2009] Cottrell, J. A., Hughes, T. J. R., Bazilevs, Y., 2009. Isogeometric analysis: toward integration of CAD and FEA. John Wiley & Sons.
  * de Boor [2001] de Boor, C., 2001. A practical guide to splines, revised Edition. Vol. 27 of Applied Mathematical Sciences. Springer-Verlag, New York.
  * Dölz et al. [2018] Dölz, J., Harbrecht, H., Kurz, S., Schöps, S., Wolf, F., 2018. A fast isogeometric BEM for the three dimensional Laplace and Helmholtz problems. Comput. Methods Appl. Mech. Engrg. 330, 83–101.
  * Duffy [1982] Duffy, M. G., 1982. Quadrature over a pyramid or cube of integrands with a singularity at a vertex. SIAM J. Numer. Anal. 19 (6), 1260–1262.
  * Falini et al. [2022] Falini, A., Giannelli, C., Kanduč, T., Sampoli, M. L., Sestini, A., 2022. A collocation IGA-BEM for 3D potential problems on unbounded domains. In: Manni, C., Speleers, H. (Eds.), Springer INdAM volume ‘Geometric Challenges in Isogeometric Analysis’. Vol. 49 of Springer INdAM Series. Springer Cham., p. In press.
  * Falini et al. [2019] Falini, A., Giannelli, C., Kanduč, T., Sampoli, M. L., Sestini, A., 2019. An adaptive IgA-BEM with hierarchical B-splines based on quasi-interpolation quadrature schemes. Int. J. Numer. Methods Engrg. 117 (10), 1038–1058.
  * Falini et al. [2020] Falini, A., Kanduč, T., Sampoli, M. L., Sestini, A., 2020. Cubature rules based on bivariate spline quasi-interpolation for weakly singular integrals. In: Fasshauer, G. E., Neamtu, M., Schumaker, L. L. (Eds.), Approximation Theory XVI. AT 2019. Vol. 336 of Springer in Mathematics & Statistics. Springer Cham, pp. 73–86.
  * Falini and Kanduč [2019] Falini, A., Kanduč, T., 2019. A study on spline quasi-interpolation based quadrature rules for the isogeometric Galerkin BEM. In: Giannelli, C., Speleers, H. (Eds.), Advanced Methods for Geometric Modeling and Numerical Simulation. Vol. 5 of Springer INdAM Series. Springer, pp. 193–227.
  * Falini et al. [2021] Falini, A., Kanduč, T., Sampoli, M. L., Sestini, A., 2021. Isogeometric BEM collocation for 3D Laplace and Helmholtz problems. in preparation.
  * Farin [2002] Farin, G., 2002. Curves and surfaces for computer-aided geometric design, 5th Edition. Computer Graphics and Geometric Modeling. Academic Press Inc., San Diego, CA.
  * Fata [2009] Fata, S. N., 2009. Explicit expressions for 3D boundary integrals in potential theory. Int. J. Numer. Meth. Eng. 78, 32–47.
  * Feischl et al. [2017] Feischl, M., Gantner, G., Haberl, A., Praetorius, D., 2017. Optimal convergence for adaptive IGA boundary element methods for weakly-singular integral equations. Numer. Math. 136, 147–182.
  * Gantner and Praetorius [2020] Gantner, G., Praetorius, D., 2020. Adaptive BEM for elliptic PDE systems, part I: abstract framework, for weakly-singular integral equations. Appl. Anal., 1–34.
  * Gao [2002] Gao, X.-W., 2002. The radial integration method for evaluation of domain integrals with boundary-only discretization. Eng. Anal. Bound. Elem. 26 (10), 905–916.
  * Gao [2010] Gao, X.-W., 2010. An effective method for numerical evaluation of general 2D and 3D high order singular boundary integrals. Comput. Methods Appl. Mech. Engrg. 199 (45), 2856–2864.
  * Ginnis et al. [2014] Ginnis, A., Kostas, K., Politis, C., Kaklis, P., Belibassakis, K., Gerostathis, T., Scott, M., Hughes, T., 2014. Isogeometric boundary-element analysis for the wave-resistance problem using T-splines. Comput. Methods Appl. Mech. Engrg. 279, 425–439.
  * Gradshteyn and Ryzhik [2007] Gradshteyn, I. S., Ryzhik, I. M., 2007. Table of integrals, series, and products, 7th Edition. Academic Press, Elsevier.
  * Guiggiani et al. [1992] Guiggiani, M., Krishnasamy, G., Rudolphi, T. J., Rizzo, F. J., 09 1992. A general algorithm for the numerical solution of hypersingular boundary integral equations. J. Appl. Mech. 59 (3), 604–614.
  * Heltai et al. [2014] Heltai, L., Arroyo, M., DeSimone, A., 2014. Nonsingular isogeometric boundary element method for Stokes flows in 3D. Comput. Methods Appl. Mech. Engrg. 268, 514–539.
  * Heltai et al. [2017] Heltai, L., Kiendl, J., DeSimone, A., Reali, A., 2017. A natural framework for isogeometric fluid–structure interaction based on BEM–shell coupling. Comput. Methods Appl. Mech. Engrg. 316, 522–546.
  * Hughes et al. [2005] Hughes, T. J. R., Cottrell, J. A., Bazilevs, Y., 2005. Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement. Comput. Methods Appl. Mech. Engrg. 194 (39-41), 4135–4195.
  * Järvenpää et al. [2003] Järvenpää, S., Taskinen, M., Ylä-Oijala, P., 2003. Singularity extraction technique for integral equation methods with higher order basis functions on plane triangles and tetrahedra. Int. J. Numer. Meth. Eng. 58 (8), 1149–1165.
  * Johnston et al. [2013] Johnston, B. M., Johnston, P. R., Elliott, D., 2013. A new method for the numerical evaluation of nearly singular integrals on triangular elements in the 3D boundary element method. J. Comput. Appl. Math. 245, 148–161.
  * Khayat et al. [2008] Khayat, M. A., Wilton, D. R., Fink, P. W., 2008. An improved transformation and optimized sampling scheme for the numerical evaluation of singular and near-singular potentials. IEEE Antennas Wirel. Propag. Lett. 7, 377–380.
  * Klaseboer et al. [2009] Klaseboer, E., Fernandez, C. R., Khoo, B. C., 2009. A note on true desingularisation of boundary integral methods for three-dimensional potential problems. Eng. Anal. Bound. Elem. 33 (6), 796–801.
  * Kostas et al. [2018] Kostas, K., Fyrillas, M., Politis, C., Ginnis, A., Kaklis, P., 2018. Shape optimization of conductive-media interfaces using an IGA-BEM solver. Comput. Methods Appl. Mech. Engrg. 340, 600–614.
  * Mazzia and Sestini [2009] Mazzia, F., Sestini, A., 2009. The BS class of Hermite spline quasi-interpolants on nonuniform knot distributions. BIT 49 (3), 611–628.
  * Mazzia and Sestini [2012] Mazzia, F., Sestini, A., 2012. Quadrature formulas descending from BS Hermite spline quasi-interpolation. J. Comput. Appl. Math. 236, 4105–4118.
  * Monegato and Sloan [1997] Monegato, G., Sloan, I. H., 1997. Numerical solution of the generalized airfoil equation for an airfoil with a flap. SIAM J. Numer. Anal. 34 (6), 2288–2305.
  * Mørken [1991] Mørken, K., 1991. Some identities for products and degree raising of splines. Constr. Approx. 7, 195–208.
  * Mousavi and Sukumar [2010] Mousavi, S. E., Sukumar, N., 2010. Generalized Duffy transformation for integrating vertex singularities. Comput. Mech. 45, 127–140.
  * Peng et al. [2017] Peng, X., Atroshchenko, E., Kerfriden, P., Bordas, S., 2017. Isogeometric boundary element methods for three dimensional static fracture and fatigue crack growth. Comput. Methods Appl. Mech. Engrg. 316, 151–185, special Issue on Isogeometric Analysis: Progress and Challenges.
  * Reid et al. [2015] Reid, M. T. H., White, J. K., Johnson, S. G., 2015. Generalized Taylor-Duffy method for efficient evaluation of Galerkin integrals in boundary-element method computations. IEEE Antennas Wirel. Propag. Lett. 63 (1), 195–209.
  * Rong et al. [2014] Rong, J., Wen, L., Xiao, J., 2014. Efficiency improvement of the polar coordinate transformation for evaluating BEM singular integrals on curved elements. Eng. Anal. Bound. Elem. 38, 83–93.
  * Sauter and Schwab [2011] Sauter, S. A., Schwab, C., 2011. Boundary element methods. Vol. 39 of Springer Series in Computational Mathematics. Springer-Verlag, Berlin, Heidelberg.
  * Schumaker [2007] Schumaker, L. L., 2007. Spline functions: basic theory, 3rd Edition. Cambridge University Press.
  * Scuderi [2009] Scuderi, L., 2009. A new smoothing strategy for computing nearly singular integrals in 3D Galerkin BEM. J. Comput. Appl. Math. 225 (2), 406–427.
  * Simpson et al. [2014] Simpson, R., Scott, M., Taus, M., Thomas, D., Lian, H., 2014. Acoustic isogeometric boundary element analysis. Comput. Methods Appl. Mech. Engrg. 269, 265–290.
  * Taus et al. [2019] Taus, M., Rodin, G. J., Hughes, T. J., Scott, M. A., 2019. Isogeometric boundary element methods and patch tests for linear elastic problems: Formulation, numerical integration, and applications. Comput. Methods Appl. Mech. Engrg. 357, 112591.
  * Taus et al. [2016] Taus, M., Rodin, G. J., Hughes, T. J. R., 2016. 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.
  * Telles [1987] Telles, J. C. F., 1987. A self-adaptive co-ordinate transformation for efficient numerical evaluation of general boundary element integrals. Int. J. Numer. Meth. Eng. 24 (5), 959–973.
  * Venås and Kvamsdal [2020] Venås, J. V., Kvamsdal, T., 2020. Isogeometric boundary element method for acoustic scattering by a submarine. Comput. Methods Appl. Mech. Engrg. 359, 112670\.
  * Zhang et al. [2015] Zhang, Y., Li, X., Sladek, V., Sladek, J., Gao, X., 2015. A new method for numerical evaluation of nearly singular integrals over high-order geometry elements in 3D BEM. J. Comput. Appl. Math. 277, 57–72.

##  Appendix A The underlying integrals

In this section we analyze some basic properties of the integrand functions
for the underlying integrals in intrinsic coordinates. Then we provide
recursive formulae to evaluate these kind of integrals analytically. Following
the isoparametric paradigm, the most common surface integrals are defined
either on rectangular or triangular domains. The definite integrals are
computed as a linear combination of indefinite ones, evaluated at appropriate
boundary points.

Suppose $f$ is a general integrable bivariate function we would like to
integrate analytically. For both type of domains we need a primitive function
of $f$, $F(x,y):=\int\int f(x,y)\mathop{dx}\mathop{dy}$. Since indefinite
integrals are not uniquely defined, we can neglect all byproducts that do not
contribute to the final values of the definite integrals, e.g., functions $g$
that satisfy $\partial^{2}g/(\mathop{\partial x}\mathop{\partial y})\equiv 0$.
For the triangular domain let us also define ${F^{[\rm t]}(y)}:=\int(\int
f(x,y)\mathop{dx})|_{x=1-y}\mathop{dy}$. The surface integral on rectangle
$[x_{0},x_{1}]\times[y_{0},y_{1}]$ can be simply computed by evaluating the
indefinite integral at the 4 corners,

$\displaystyle\int_{y_{0}}^{y_{1}}\int_{x_{0}}^{x_{1}}f(x,y)\mathop{dx}\mathop{dy}=F(x_{1},y_{1})-F(x_{0},y_{1})-F(x_{1},y_{0})+F(x_{0},y_{0}).$

To compute the surface integral of $\tilde{f}$ on a general triangle $T$, a
preliminary step is to apply an affine mapping $\varphi$ to map, for example,
a triangle with vertices $(0,0),(1,0),(0,1)$ to $T$. Then we need to again
insert appropriate boundary values into the two indefinite integrals,

$\displaystyle\int_{T}\tilde{f}(\tilde{x},\tilde{y})\,dT=\int_{0}^{1}\int_{0}^{1-y}f(x,y)\mathop{dx}\mathop{dy}=F(0,0)-F(0,1)+F^{[\rm
t]}(1)-F^{[\rm t]}(0),$

where $f:=(\tilde{f}\circ\varphi)\cdot J$ and constant $J$ is the area of the
infinitesimal surface element.

The two derived formulae allow us to analytically compute surface integrals on
a rectangle and triangle, if the primitive function of $f$ is known. In the
remaining part of the section we focus only on a special family of these
functions that are the building blocks for the singularity extraction
technique, and provide formulae for their indefinite integrals.

###  A.1 Main properties

Our goal in this section is to derive recursive formulae to compute integrals
of a type

$\displaystyle I_{p,q,r}=I_{p,q,r}(x,y):=\int\int
R^{p}x^{q}y^{r}\mathop{dx}\mathop{dy},$ (33)

where $R=R(a,b,c,(x,y)):=\sqrt{ax^{2}+bxy+cy^{2}}$ is a function,
corresponding to the first fundamental form of a smooth surface. Here we
assume $a,c>0$ and $b\in{\mathop{{\rm I}\kern-1.99997pt{\rm R}}\nolimits}$,
such that $4ac-b^{2}>0$, which ensures the quadratic function to be positive
and hence $R>0$. Then, $p$ is a negative odd integer, while $q,r$ are non-
negative integers. Double integrals (33) can be computed via analytical
formulae for single inner and outer integrals,

$\displaystyle I_{p,q}^{(1)}=I_{p,q}^{(1)}(x):=\int
R^{p}x^{q}\mathop{dx},\qquad I_{p,r}^{(2)}=I_{p,r}^{(2)}(y):=\int
R^{p}y^{r}\mathop{dy},$ (34)

namely $I_{p,q,r}=\int y^{r}I_{p,q}^{(1)}\mathop{dy}.$ To compute integrals on
a triangle, we also need formulae to compute

$\displaystyle I_{p,q,r}^{[\rm t]}=I_{p,q,r}^{[\rm t]}(y):=\int
y^{r}(I_{p,q}^{(1)}|_{x=1-y})\mathop{dy}.$ (35)

######  Remark 8.

Since $x$ and $y$ appear symmetrically in (33) it is sufficient to study
analytical formuale for $I_{p,q}^{(1)}$ only. Formulae for $I_{p,r}^{(2)}$ are
therefore simply transformed from the derived expressions for $I_{p,q}^{(1)}$
by exchanging variables $x$ and $y$, and parameters $a$ and $c$, and $q$ and
$r$.

Smoothness (regularity) of the integrand $R^{p}x^{q}y^{r}$ in (33) at
$(x,y)=(0,0)$ is easily characterized by the parameter $\zeta$,

$\displaystyle\zeta(R^{p}x^{q}y^{r}):=p+q+r.$ (36)

For the study of the smoothness we relax the condition on $p$ to be an
arbitrary integer. When we consider a sum of integrands of this type, it is
also convenient to use the following definition

$\displaystyle\zeta\left(\sum_{\ell}c_{\ell}R^{p_{\ell}}x^{q_{\ell}}y^{r_{\ell}}\right):=\min_{\ell}\zeta(R^{p_{\ell}}x^{q_{\ell}}y^{r_{\ell}}),$

where we assume that the sum is written in irreducible form for some
coefficients $c_{\ell}$; it is also meaningful to set $\zeta=\infty$ if we sum
over an empty set of indices.

Without a proof we mention two simple properties of the parameter $\zeta$.

######  Proposition 9.

Let $f,g$ be functions of a type $R^{p}x^{q}y^{r}$. Then

  * 1.

$\zeta(fg)=\zeta(f)+\zeta(g)$,

  * 2.

if $\zeta(f)=m$, then $n$-th derivative of $f$ can be written as
${\partial^{n}f}/{(\mathop{\partial x}^{i}\mathop{\partial
y}^{n-i})}=\sum_{\ell=0}^{n}f_{\ell}$, $i=0,1,\dots,n$, where $f_{\ell}$ are
also of a type $R^{p}x^{q}y^{r}$ and $\zeta(f_{\ell})=m-n$ for every $\ell$.

Before analyzing the smoothness of the integrands, we need to prove the
following lemma.

######  Lemma 10.

Function $\gamma\mapsto(a+b\gamma+c\gamma^{2})^{p/2}$ is positive and bounded
for $|\gamma|\leq 1$. Namely, there exist two positive constants
$C_{abcp,1},C_{abcp,2}$, independent of $\gamma$, such that

$\displaystyle C_{abcp,1}\leq(a+b\gamma+c\gamma^{2})^{p/2}\leq C_{abcp,2}.$

Proof. There exist two positive constants
$\tilde{C}_{abc,1},\tilde{C}_{abc,2}$, independent of $\gamma$, such that

$\displaystyle\tilde{C}_{abc,1}\leq{a+b\gamma+c\gamma^{2}}\leq\tilde{C}_{abc,2}.$
(37)

The upper bound $\tilde{C}_{abc,2}$ exists since $\gamma\mapsto
a+b\gamma+c\gamma^{2}$ is continuous and thus bounded on $|\gamma|\leq 1$. The
existence of the lower bound $\tilde{C}_{abc,1}$ follows from assumption
$0<4ac-b^{2}$. By raising (37) to the power $p/2$ we get the sought
expression. ∎

######  Proposition 11.

Function $R^{p}x^{q}y^{r}$ is bounded if $\zeta(R^{p}x^{q}y^{r})=0$.

Proof. Let us split ${\mathop{{\rm I}\kern-1.99997pt{\rm R}}\nolimits}^{2}$
into eight subdomains, formed by cutting the plane with lines $x=0$, $y=0$,
$y=x$ and $y=-x$. Let us prove the proposition only for the subdomain

$\displaystyle T_{1}:=\\{(x,y)\in{\mathop{{\rm I}\kern-1.99997pt{\rm
R}}\nolimits}^{2}:0\leq y\leq x\\},$ (38)

since for the other 7 cases the proof is similar.

Let us choose a point $(x,y)\in T_{1}$. By writing $y=\gamma x$ for
$0\leq\gamma\leq 1$ we get

$\displaystyle
R^{p}x^{q}y^{r}={(a+b\gamma+c\gamma^{2})}^{p/2}\gamma^{r}x^{p+q+r}.$

By considering $p+q+r=0$, from Lemma 10 there exist positive constants
$C_{abcp,1},C_{abcp,2}$ such that

$\displaystyle 0\leq
C_{abcp,1}\gamma^{r}\leq{(a+b\gamma+c\gamma^{2})}^{p/2}\,\gamma^{r}\leq
C_{abcp,2}\gamma^{r}\leq C_{abcp,2}$

for all $0\leq\gamma\leq 1$. Since this estimate holds for every $(x,y)\in
T_{1}$, the function $R^{p}x^{q}y^{r}$ is bounded on $T_{1}$.

∎

######  Proposition 12.

Let $f(x,y):=R^{p}x^{q}y^{r}$ with $\zeta(f)\geq 1$. Then $f$ is $C^{0}$
continuous and $f(0,0)=0$.

Proof. The function is analytical on ${\mathop{{\rm I}\kern-1.99997pt{\rm
R}}\nolimits}^{2}\backslash\\{(0,0)\\}$, hence it is enough to check it is
continuous at $(0,0)$. Fix sufficiently small $0<\delta$. Let us show that
there exists sufficiently small $\varepsilon$, $0\leq\varepsilon\leq 1$, such
that from $\|(x,y)\|_{2}\leq\varepsilon$ it follows
$|R^{p}x^{q}y^{r}|\leq\delta$.

By splitting ${\mathop{{\rm I}\kern-1.99997pt{\rm R}}\nolimits}^{2}$ again
into eight subdomains and assuming $(x,y)\in T_{1}$ (see (38)), we can write
$y=\gamma x$ for $0\leq\gamma\leq 1$. Then we can again write the function as

$\displaystyle
R^{p}x^{q}y^{r}={(a+b\gamma+c\gamma^{2})}^{p/2}\gamma^{r}x^{p+q+r}.$

From Lemma 10 there exists a positive constant $C_{abcp,2}$, independent of
$\gamma$, and we can take $\varepsilon$ small enough, so that

$\displaystyle R^{p}x^{q}y^{r}\leq C_{abcp,2}\,\varepsilon^{p+q+r}\leq
C_{abcp,2}\,\varepsilon\leq\delta.$

A similar proof can be asserted for the other 7 subdomains.

∎

######  Proposition 13.

Let $f(x,y):=R^{p}x^{q}y^{r}$ with $\zeta(f)\geq n+1$. Then $f$ is $C^{n}$
continuous and ${\partial^{n}f}/{(\mathop{\partial x}^{i}{\mathop{\partial
y}}^{n-i})}(0,0)=0$ for $i=0,1,\dots,n$.

Proof.

Let $\zeta(f)=p+q+r\geq 2$. It is straightforward to check that

$\displaystyle\frac{\partial}{\partial
x}(R^{p}x^{q}y^{r})=paR^{p-2}x^{q+1}y^{r}+\frac{pb}{2}R^{p-2}x^{q}y^{r+1}.$
(39)

From Proposition 12 it follows that both functions on the right-hand side of
(39) are continuous. A similar argument can be stated for
$\partial(R^{p}x^{q}y^{r})/{\partial y}$. Thus $R^{p}x^{q}y^{r}$ is a $C^{1}$
smooth function.

If $p+q+r\geq n+1$, we can apply formula (39), or its variation for the
derivative in variable $y$, $n$ times and get

$\displaystyle\frac{\partial^{n}}{\mathop{\partial x}^{i}\mathop{\partial
y}^{n-i}}(R^{p}x^{q}y^{r})=\sum_{j=0}^{n}c_{j}R^{p-2n}x^{q+j}y^{r+n-j}$ (40)

for suitable coefficients $c_{j}$ and $i=0,1,\dots,n$. Since
$\zeta(R^{p-2n}x^{q+j}y^{r+n-j})\geq 1$ and by applying Proposition 12,
functions on the right-hand side of (40) are continuous for all mixed $n$-th
derivatives, hence $R^{p}x^{q}y^{r}$ is a $C^{n}$ smooth function. From
Proposition 12 also follows directly the condition for the all $n$-th
derivatives to be zero at $(0,0)$. ∎

###  A.2 Recursive formulae for $I_{p,q,r}$ and $I_{p,q,r}^{[\rm t]}$

Formulae for double integrals are derived in two steps. First, we focus on the
inner integral of (33), $I_{p,q}^{(1)}$. The outer integral is computed
similarly, by using formulae for the i̧nner integral, and by suitably
replacing the parameters $a,c$ and $q,r$ and variables $x,y$.

First, we examine a border case of $I_{p,q}^{(1)}$ in (34) when $p=-1$.

######  Lemma 14.

Let $q\geq 0$. Then

$\displaystyle I_{-1,q}^{(1)}=RP_{q-1}+C_{q-1}y^{q}I_{-1,0}^{(1)},$ (41)

where

$\displaystyle I_{-1,0}^{(1)}=\frac{1}{\sqrt{a}}\log(2\sqrt{a}R+2ax+by).$ (42)

The bivariate homogeneous polynomials $P_{q}=P_{q}(x,y)$ of degree $q$ and
coefficients $C_{q}$ are defined as

$\displaystyle P_{q}=\left\\{\begin{array}[]{ll}0,&q=-1\\\ \frac{1}{a},&q=0\\\
\frac{1}{a(q+1)}\left(x^{q}-\frac{b(2q+1)}{2}yP_{q-1}-cqy^{2}P_{q-2}\right),&q\geq
1\end{array}\right.,\qquad C_{q}=\left\\{\begin{array}[]{ll}1,&q=-1\\\
-\frac{b}{2a},&q=0\\\
\frac{1}{a(q+1)}\left(-\frac{b(2q+1)}{2}C_{q-1}-cqC_{q-2}\right),&q\geq
1\end{array}\right..$

Proof. We refer to [28], Section 2.26 for the special case $q=0$,
$I_{-1,0}^{(1)}=\int R^{-1}\mathop{dx}$. From the same reference we also know
that

$\displaystyle I_{-1,1}^{(1)}=\frac{1}{a}R-\frac{b}{2a}yI_{1,0}^{(1)}.$

When $q\geq 2$ the considered integral can be rewritten in the following way

$\displaystyle a\int R^{-1}{x^{q}}\mathop{dx}$ $\displaystyle=\int
R^{-1}x^{q-2}(ax^{2}+(bxy+cy^{2})-(bxy+cy^{2}))\mathop{dx}$
$\displaystyle=\int Rx^{q-2}\mathop{dx}-\int
R^{-1}x^{q-2}(bxy+cy^{2})\mathop{dx}.$ (43)

For the first integral in (A.2) we apply integration by parts,

$\displaystyle I_{1,q-2}^{(1)}=\int
Rx^{q-2}\mathop{dx}=\frac{1}{q-1}Rx^{q-1}-\frac{1}{2(q-1)}\int
R^{-1}x^{q-1}(2ax+by)\mathop{dx}.$ (44)

By combining recursive formulae (A.2) and (44) we obtain the formula for
$q\geq 2$,

$\displaystyle
I_{-1,q}^{(1)}=\frac{1}{aq}\left(Rx^{q-1}-\frac{b(2q-1)}{2}yI_{-1,q-1}^{(1)}-c(q-1)y^{2}I_{-1,q-2}^{(1)}\right).$

By considering two stopping cases of the recursion, $q=0,1$, after some
simplification we can express $I_{-1,q}^{(1)}$ with formula (41).

∎

Let us analyze another border case of $I_{p,q}^{(1)}$, when $p\leq-3$ and
$q=0,1$.

######  Lemma 15.

Let $p\leq-3$. Then

$\displaystyle
I_{p,0}^{(1)}=\sum_{\ell=p+2,p+4,\dots,-1}\bigg{(}\prod_{\ell_{2}=p,p+2,\dots,\ell-4}D_{\ell_{2}}\bigg{)}\bigg{(}\prod_{\ell_{2}=p,p+2,\dots,\ell-2}E_{\ell_{2}}\bigg{)}R^{\ell}(2ax+by)y^{p-\ell},$
(45)

where

$\displaystyle D_{\ell}=-a(2\ell+6),\qquad
E_{\ell}=-\frac{2}{(4ac-b^{2})(\ell+2)},$

and

$\displaystyle
I_{p,1}^{(1)}=\frac{1}{a(p+2)}R^{p+2}-\frac{b}{2a}yI_{p,0}^{(1)}.$ (46)

Proof. For $p\leq-3$ we integrate $I_{p+2,0}^{(1)}$ by parts,

$\displaystyle 2a\int R^{p+2}\mathop{dx}=R^{p+2}(2ax+by)-\frac{p+2}{2}\int
R^{p}(2ax+by)^{2}\mathop{dx}.$

The expression can be rewritten as

$\displaystyle 0$ $\displaystyle=R^{p+2}(2ax+by)-\frac{p+2}{2}\int
R^{p}(2ax+by)^{2}\mathop{dx}-2a\int R^{p+2}\mathop{dx}$
$\displaystyle=R^{p+2}(2ax+by)-\frac{p+2}{2}\int
R^{p}\left((2ax+by)^{2}-4aR^{2}\right)\mathop{dx}-a(2(p+2)+2)\int
R^{p+2}\mathop{dx}$
$\displaystyle=R^{p+2}(2ax+by)+\frac{(4ac-b^{2})y^{2}(p+2)}{2}\int
R^{p}\mathop{dx}-a(2p+6)\int R^{p+2}\mathop{dx},$

hence

$\displaystyle I_{p,0}^{(1)}$ $\displaystyle=\int
R^{p}\mathop{dx}=\frac{2}{(4ac-b^{2})(p+2)}\left(a(2p+6)y^{-2}I_{p+2,0}^{(1)}-R^{p+2}(2ax+by)y^{-2}\right).$

By continuing with the recursion for higher $p$ we obtain after some
simplifications the expression (45).

From

$\displaystyle R^{p+2}=\int\frac{\partial}{\partial
x}R^{p+2}\mathop{dx}=\int\frac{p+2}{2}R^{p}(2ax+by)\mathop{dx}=a(p+2)I_{p,1}^{(1)}+\frac{(p+2)b}{2}yI_{p,0}^{(1)}$

we get formula (46).

∎

When $p\leq-3$ and $2\leq q$ the integral $I_{p,q}^{(1)}$ is computed
recursively from integrals of higher and lower powers $p$ and $q$,
respectively, by applying formula
$aR^{p}x^{q}=R^{p}x^{q-2}(ax^{2}+(bxy+cy^{2})-(bxy+cy^{2}))$,

$\displaystyle
I_{p,q}^{(1)}=\frac{1}{a}I_{p+2,q-2}^{(1)}-\frac{b}{a}yI_{p,q-1}^{(1)}-\frac{c}{a}y^{2}I_{p,q-2}^{(1)},$
(47)

see Fig. 8. The recursion terminates when either the case $p=-1$ is reached
and we use the formula (41), or the case $p\leq-3,q=0,1$ is encountered and we
use the expressions (45) and (46), respectively.

$I_{p,q}^{(1)}$ $I_{p,q-1}^{(1)}$ $I_{p,q-2}^{(1)}$ $I_{p+2,q-2}^{(1)}$
formula (41) for $I_{-1,q}^{(1)}$ formula (45) for $I_{p,1}^{(1)}$ formula
(46) for $I_{p,0}^{(1)}$ Figure 8: Illustration of the recursion for
$I_{p,q}^{(1)}$ in (47).

We are left to show how to compute the outer integration of $I_{p,q,r}$, using
similar formulae to the inner integration.

######  Proposition 16.

Integral $I_{p,q,r}$ is a linear combination of integrals
$y^{\zeta_{0}+2}I_{-1,0}^{(1)}$ and
$x^{\zeta_{0}-\hat{p}-\hat{r}+1}I_{\hat{p},\hat{r}}^{(2)}$ for
$\hat{p}=p+2,p+4,\dots,-1$, $\hat{r}=0,1,\dots,\zeta_{0}-\hat{p}+1$, where
$\zeta_{0}:=p+q+r$ .

Proof. Using recursion (47) on $y^{r}I_{p,q}^{(1)}$ we obtain

$\displaystyle I_{p,q,r}=\int
y^{r}I_{p,q}^{(1)}\mathop{dy}=\sum_{\hat{q}=0,1,\dots,q}\alpha_{\hat{q}}\int
y^{\zeta_{0}-\hat{q}+1}I_{-1,\hat{q}}^{(1)}\mathop{dy}+\sum_{\hat{p}=p,p+2,\dots,-1}\int\beta_{\hat{p}}y^{\zeta_{0}-\hat{p}}I_{\hat{p},0}^{(1)}+\gamma_{\hat{p}}y^{\zeta_{0}-\hat{p}-1}I_{\hat{p},1}^{(1)}\mathop{dy},$
(48)

for suitable coefficients $\alpha_{\hat{p}},\beta_{\hat{p}},\gamma_{\hat{q}}$.
Observe that $\zeta$ parameter is preserved (equal to $\zeta_{0}$) for all
integrands inside double integrals on the right-hand side of (48).

By applying formulae (45) and (46) on the second term of (48), the term writes
as

$\displaystyle\sum_{\hat{p}=p+2,p+4,\dots,-1}\int
R^{\hat{p}}P_{\zeta_{0}-\hat{p}+1}^{[1]}\mathop{dy}$ (49)

for some bivariate homogeneous polynomials $P_{\zeta_{0}-\hat{p}+1}^{[1]}$ of
degree $\zeta_{0}-\hat{p}+1$, hence the sum is a linear combination of
$x^{\zeta_{0}-\hat{p}-\hat{r}+1}I_{\hat{p},\hat{r}}^{(2)}$ for
$\hat{p}=p+2,p+4,\dots,-1$ and $\hat{r}=0,1,\dots,\zeta_{0}-\hat{p}+1$.

On each integral of the first term of (48) we apply formula (41). Each non-
logarithmic term can be rewritten as

$\displaystyle\int RP_{\zeta_{0}}^{[2]}\mathop{dy}=\int
R^{-1}(R^{2}P_{\zeta_{0}}^{[2]})\mathop{dy}=:\int
R^{-1}P_{\zeta_{0}+2}^{[3]}\mathop{dy},$

for suitable polynomials $P_{\zeta_{0}}^{[2]},P_{\zeta_{0}+2}^{[3]}$, hence it
can be treated in the same manner as the addend in (49) for $\hat{p}=-1$. The
remaining integrals of (48) are integrated by parts using

$\displaystyle\int
y^{\zeta_{0}+1}I_{-1,0}^{(1)}\mathop{dy}=\frac{y^{\zeta_{0}+2}}{\zeta_{0}+2}I_{-1,0}^{(1)}-\frac{1}{\zeta_{0}+2}\int\left(\frac{1}{\sqrt{a}}-xR^{-1}\right)y^{\zeta_{0}+1}\mathop{dy}=\frac{1}{\zeta_{0}+2}\left(y^{\zeta_{0}+2}I_{-1,0}^{(1)}-\frac{y^{\zeta_{0}+2}}{\sqrt{a}(\zeta_{0}+2)}+xI_{-1,\zeta_{0}+1}^{(2)}\right).$

We can ignore the second addend since it does not depend on variable $x$.

∎

The integral $I_{p,q,r}^{[\rm t]}$ from (35) is computed similarly.

######  Corollary 17.

Let $\hat{a}=a$, $\hat{b}:=(-2a+b)$, $\hat{c}:=a+b+c$. Then integral
$I_{p,q,r}^{(t)}$ is a linear combination of integrals
$y^{\zeta_{0}+2}I_{-1,0}^{(1)}|_{x=1-y}$ and $I_{\hat{p},\hat{r}}^{(2)}$,
involving integrand function $R(\hat{a},\hat{b},\hat{c},(1,y))$,
$\hat{p}=p+2,p+4,\dots,-1$, $\hat{r}=0,1,\dots,\zeta_{0}-\hat{p}+1$ and
$\zeta_{0}:=p+q+r$.

Proof. By applying similar arguments as in the proof of Proposition 16,
$I_{p,q,r}^{(t)}=\int y^{r}I_{p,q}^{(1)}|_{x=1-y}\mathop{dy}$ is computed from
relation (48) and by using $x=1-y$. Instead of (49) we get

$\displaystyle\sum_{\hat{p}=p+2,p+4,\dots,-1}\int(R^{\hat{p}}P_{\zeta_{0}-\hat{p}+1}^{[1]})|_{x=1-y}\mathop{dy}=\sum_{\hat{p}=p+2,p+4,\dots,-1}\int{R(\hat{a},\hat{b},\hat{c},(1,y))}^{\hat{p}}P_{\zeta_{0}-\hat{p}+1}^{[1]}|_{x=1-y}\mathop{dy},$
(50)

Thus we can write each term in the sum (50) as a linear combination of
$I_{\hat{p},\hat{r}}^{(2)}$, for $\hat{p}=p+2,p+4,\dots,-1$,
$\hat{r}=0,1,\dots,\zeta_{0}-\hat{p}+1$ and for the modified parameters
$\hat{a},\hat{b},\hat{c}$.

For the remaining integral it holds

$\displaystyle\int
y^{\zeta_{0}+1}I_{-1,0}^{(1)}|_{x=1-y}\mathop{dy}=\frac{y^{\zeta_{0}+2}}{\zeta_{0}+2}I_{-1,0}^{(1)}|_{x=1-y}-\frac{1}{\zeta_{0}+2}\int\left(\frac{1}{\sqrt{a}}-R(\hat{a},\hat{b},\hat{c},(1,y))^{-1}\right)y^{\zeta_{0}+1}\mathop{dy}$

and

$\displaystyle
I_{-1,0}^{(1)}|_{x=1-y}=\frac{1}{\sqrt{a}}\log(2\sqrt{a}R(\hat{a},\hat{b},\hat{c},(1,y))+(-2a+b)y+2a).$

∎