An efficient method for evaluating BEM singular integrals on curved elements
with application in acoustic analysis
Junjie Rong, Lihua Wen∗, Jinyou Xiao
College of Astronautics, Northwestern Polytechnical University, Xi’an 710072, P. R. China

arXiv:1306.0282v1 [] 3 Jun 2013

Abstract
The polar coordinate transformation (PCT) method has been extensively used to treat various singular integrals in
the boundary element method (BEM). However, the resultant integrands of the PCT tend to become nearly singular
when (1) the aspect ratio of the element is large or (2) the field point is closed to the element boundary; thus a large
number of quadrature points are needed to achieve a relatively high accuracy. In this paper, the first problem is
circumvented by using a conformal transformation so that the geometry of the curved physical element is preserved
in the transformed domain. The second problem is alleviated by using a sigmoidal transformation, which makes the
quadrature points more concentrated around the near singularity.
By combining the proposed two transformations with the Guiggiani’s method in [M. Guiggiani, et al. A general algorithm for the numerical solution of hypersingular boundary integral equations. ASME Journal of Applied
Mechanics, 59(1992), 604-614], one obtains an efficient and robust numerical method for computing the weakly-,
strongly- and hyper-singular integrals in high-order BEM with curved elements. Numerical integration results show
that, compared with the original PCT, the present method can reduce the number of quadrature points considerably,
for given accuracy. For further verification, the method is incorporated into a 2-order Nyström BEM code for solving
acoustic Burton-Miller boundary integral equation. It is shown that the method can retain the convergence rate of the
BEM with much less quadrature points than the existing PCT. The method is implemented in C language and freely
available.
Keywords: singular integrals, boundary element method, Nyström method, acoustics

1. Introduction
The boundary element method (BEM) has been a most important numerical method in science and engineering.
Its unique advantages includes the highly accurate solution on the boundary, the reduction of dimensionality, and
the incomparable superior in solving infinite or semi-infinite field problems, etc. Historically, relatively low-order
discretizations have been used with geometries modelled using first order elements and surface variables modelled to
zero or first order on those elements. Recently, however, there has been increasing interest in the use of high order
methods, in order to obtain extra digits of precision with comparatively small additional effort. Successful usage
are seen in acoustics [1], electromagnetics [2, 3], elasticity [4], aerodynamics [5], to name a few. One of the main
difficulties in using high order BEM is the lack of efficient methods for accurately evaluating various singular integrals
over curved elements. For the well-developed low order methods, there is no great difficulty since robust analytical
and numerical integration schemes are generally available for planar elements [6, 7, 9]. When concerned with high
order elements, however, a fully numerical method is required; various techniques have been proposed, for example,
the singularity subtraction [8], special purpose quadrature [10, 3, 11, 12] and the variable transformation methods
[13, 14, 15, 16].
∗ Corresponding

author
Email addresses: rxrjj@126.com (Junjie Rong), lhwen@nwpu.edu.cn (Lihua Wen), xiaojy@nwpu.edu.cn (Jinyou Xiao)

Preprint submitted to Elsevier

June 4, 2013

In [8], Guiggiani et al proposed a unified formula (Guiggiani’s method for short) for treating various order of
singular integrals on curved elements. It is a singularity subtraction method, and has found extensive use in BEM. In
this method the singular parts are extracted from the integrand and treated analytically. The remaining parts are regular
so that conventional Gaussian quadrature can be employed, but the number of quadrature points needed would be large,
depending on the regularity of the associated integrands. The special purpose quadrature can be used to substantially
reduce the number of quadrature points, besides it is more robust and highly accurate. Most recently, James Bremer
[3] proposed such a method which can achieve machine precision. The problem with the special purpose quadrature
is that its construction somewhat complicated and time-consuming. Variable transformation method, also known as
singularity cancelation, eliminate the singularity of the integrand by the null Jacobian at the field point through a
proper change of variables. Although simple to implement, it is generally hard to be used in handling hypersingular
integrals.
In most existing integration methods, including those mentioned above, polar coordinates transformation (PCT)
always serves as a common base [8, 3, 12, 13, 16, 21]. It converses the surface integral into a double integral in radial
and angular directions. Many works have been done on dealing with the singularity in the radial direction; numerical
integration on the angular direction, however, still deserves more attention. In fact, after singularity cancelation or
subtraction, although the integrand may behaves very well in the radial direction, its behavior in the angular direction
would be much worst, so that too many quadrature points are needed. Especially when the field point lies close
to the boundary of the element, one can clearly observe near singularity of the integrand in the angular direction.
Unfortunately, this case is frequently encountered in using high order elements, especially for non-conform elements
as will be demonstrated by a Nyström BEM in this paper. Similar problems have been considered in work about
the nearly singular BEM integrals. Effective methods along this line are the subdivision method [17], the Hayami
transformation [18], the sigmoidal transformation [19], etc. For the singular BEM integrals, the angular transformation
has been used for computing weakly singular integrals over planar element [13].
Another problem of the PCT is how to find a proper planar domain to establish the polar coordinates. In usual
practice, the integral is carried out over standard reference domain (i.e., triangle in this paper) in intrinsic coordinates.
However, the reference triangle is independent of the shape of the curved element, and thus the distortion of the
element is brought into the integrand, which can cause near singularity in the angular direction. Consequently, the
performance of quadratures is highly sensitive to the shape of element. The above two problems in the angular
direction were considered and resolved by special purpose quadrature rule in [3]; whereas, as mentioned before, the
algorithm has its own overhead.
In this paper, two strategies are proposed to overcome two problems in the angular direction, respectively. First, a
conformal transformation is carried out to map the curved physical element onto a planar integration domain. Since it
is conformal at the field point, the resultant integration domain perseveres the shape of the curved element. Second, a
new sigmoidal transformation is introduced to alleviate the near singularity caused by the closeness of the field point
to the element boundary. The two proposed techniques, when combined with Guiggiani’s method, lead to a unified,
efficient and robust numerical integration methods for various singular integrals in high order BEM. As a byproduct
of the conformal transformation, the line integral in Guiggiani’s method can be evaluated in close form.
The paper is organized as follows. The Nyström BEM for acoustics and singular integrals encountered are reviewed in section 2. Section 3 describes the Guiggiani’s unified framework for treating various order of singular
integrals, with emphasis on the two reasons which render the poor performance of polar coordinates methods. Two
efficient transformations which represent the main contribution of this paper are proposed in section 4. Numerical
examples are given in section 5 to validate the efficiency and accuracy of the present methods. Section 6 concludes
the paper with some discussions.
2. Model problem: acoustic BIEs and Nyström discretization
The method presented in this paper will be used in solving the acoustic Burton-Miller BIE, which is briefly
recalled here. The BIE is solved by using the Nyström method with the domain boundary being partitioned into
curved quadratic elements. Various singular integrals that will be treated in the following sections are summarized.

2

2.1. Acoustic Burton-Miller formulation
The time harmonic acoustic waves in a homogenous and isotropic acoustic medium Ω is described by the following
Helmholtz equation
∇2 u(x) + k2 u(x) = 0, ∀x ∈ Ω,
(1)

where, ∇2 is the Laplace operator, u(x) is the sound pressure at the point x = (x1 , x2 , x3 ) in the physical coordinate
system, k = ω/c is the wavenumber, with ω being the angular frequency and c being the sound speed. For static case
k = 0, (1) becomes the Laplace equation. By using Green’s second theorem, the solution of Eq. (1) can be expressed
by integral representation
Z
Z
∂G(x, y)
u(x) +
G(x, y)q(y)dS (y) + uI (x), ∀x ∈ Ω,
(2)
u(y)dS (y) =
S ∂n(y)
S

where x is a field point and y is a source point on boundary S ; q(y) = ∂u(y)/∂n(y) is the normal gradient of sound
pressure; n(y) denotes the unit normal vector at the source point y. The incident wave uI (x) will not be presented for
radiation problems. The three dimensional fundamental solution G is given as
G(x, y) =

eikr
,
4πr

(3)

with r = |x − y| denoting the distance between the source and the field points.
Before presenting the BIEs it is convenient to introduce the associated single, double, adjoint and hypersingular
layer operators which are denoted by S, D, M and H, respectively; that is,
Z
G(x, y)q(y)dS (y),
(4a)
S q(y) =
ZS
∂G(x, y)
D u(y) =
u(y)dS (y),
(4b)
S ∂n(y)
Z
∂G(x, y)
q(y)dS (y),
(4c)
M q(y) =
S ∂n(x)
Z
∂2G(x, y)
H u(y) =
u(y)dS (y).
(4d)
S ∂n(x)∂n(y)
The operator S is weakly singular and the integral is well-defined, while the operators D and M are defined in Cauchy
principal value sense (CPV). The operator H, on the other hand, is hypersingular and unbounded as a map from the
space of smooth functions on S to itself. It should be interpreted in the Hadamard finite part sense (HFP). Denoting a
vanishing neighbourhood surrounding x by sε , the CPV and HFP integrals are those after extracting free terms from a
limiting process to make sε tends to zero in deriving BIEs [8].
Letting the field point x approach the boundary S in Eq. (2) leads to the CBIE
C(x)u(x) + Du(x) = Sq(x) + uI (x),

x ∈ S,

(5)

where, C(x) is the free term coefficient which equals to 1/2 on smooth boundary. By taking the normal derivative of
Eq. (2) and letting the field point x go to boundary S , one obtains the hypersingular BIE (HBIE)
C(x)q(x) + Hu(x) = Mq(x) + qI (x),

x ∈ S.

(6)

Both CBIE and HBIE can be applied to calculate the unknown boundary values of interior acoustic problems. For an
exterior problem, they have a different set of fictitious frequencies at which a unique solution can’t be obtained. However, Eqs. (5) and (6) will always have only one solution in common. Given this fact, the Burton-Miller formulation
which is a linear combination of Eqs. (5) and (6) (CHBIE) should yield a unique solution for all frequencies [20]
h
i
C(x)u(x) + (D +α H) u(y) − uI (x) = (S +α M) q(y) − α c(x)q(x) − qI (x) , x ∈ S ,
(7)
where, α is a coupling constant that can be chosen as i/k.

3

2.2. Nyström method and the singular integrals
In this paper, the Nyström method is used to discretize the BIEs. Let K be one of the integral operators in (4)
and K be the associated kernel function. Divide the problem into two regions, a region near and far from the field
point x. If the element locates in far field, Nyström method replaces the integral operator K with a summation under
a quadrature rule
Z
X
K(xi , y)u(y)dS (y) 
ω j K(xi , y j )u(y j), ∆S ∈ S \Di ,
(8)
∆S

j

where, xi and y j are quadrature points, Di is the near field of xi , ω j is the jth weight over element ∆S . Such quadrature
rules can be obtained by mapping the Gaussian quadrature rules onto the parameterization of ∆S .
If the element locates in near field, however, the kernels exhibit singularities or even hypersingularities. As a
result, conventional quadratures fail to give correct results. In order to maintain high-order properties, the quadrature
weights are adjusted by a local correction procedure. Thus (8) becomes
Z
X
K(xi , y)u(y)dS (y) 
ω̄ j (xi )u(y j), ∆S ∈ Di ,
(9)
∆S

j

where ω̄ j (xi ) represents the modified quadrature weights for specialized rule at the singularity. The local corrected
procedure is performed by approximating the unknown quantities using linear combination of polynomial basis functions which are defined on intrinsic coordinates (Fig.1). Modified weights are obtained by solving the linear system
Z
X
(n)
K(xi , y)φ(n) (y)dS (y), ∆S ∈ Di ,
(10)
ω̄ j φ (y j ) =
∆S

j

where φ(n) are polynomial basis functions. For 2th Nyström method used in this paper, φ(n) employed are of the form
φ(n) (ξ1 , ξ2 ) = ξ1p ξ2q ,

p + q ≤ 2,

(11)

where, p and q are integers, ξ1 and ξ2 denote local intrinsic coordinates.
The right hand side integrals in Eq. (10) are of crucial importance to the accuracy of the Nyström method. They
are referred to as singular integrals when xi lies on the element, and nearly singular integrals when xi is closed to but
not on the element. This paper deals with the integrals in the first case, the nearly singular integrals are treated via an
recursive subdivision quadrature. When concerned with the first three operators in equation (4), the integrals have a
weak singularity of r−1 , while the other operator H is hypersingular of order r−3 , as r → 0.
3. An unified framework for singular integrals
Since various order of singularities appear in Eq. (7) or many other BIEs, it is advantage to find a unified formula
to treat these integrals in the same framework. In such a way, these integrals can be implement in just one program
so that the computational cost will be reduced. By expansion of the singular integrands in polar coordinates, the three
types of singular integrals considered in this paper can be handled in a unified manner by using the formula proposed
by Guiggiani [8]. Despite of this advantage the method can be of low efficiency in practical usage; the very reasons
are then explained.
3.1. Polar coordinates transformation
Following a common practice in the BEM, the curved element is first mapped onto a region ∆ of standard shape
in the parameter plane (Fig. 1). In this case, the integral must be evaluated in Eq. (10) are of the form
Z
K(x, y(ξ))φ(ξ) |J(ξ)| dξ1 dξ2 ,
(12)
I=
∆

where |J| is the transformation Jacobian from global coordinates to the intrinsic coordinates,
#
"
∂y
∂y
∂y
∂y
,
|J(ξ)| =
.
×
J=
∂ξ1 ∂ξ2
∂ξ1 ∂ξ2
4

ξ2
1

2
u2
4x

0.5

5

(0.5, 0.5)

u1
6

1

3

0

0.5

1 ξ1

Figure 1: Description of curved triangle. Left hand image: curved triangle and a field point; right hand image:the reference triangle and typical
field points distribution of 2th order Nyström method in intrinsic coordinate system

Then, polar coordinates (ρ, θ) centered at ξ s (the image of x on intrinsic plane) are defined in the parameter space
(Fig. 2)

(s)


 ξ1 = ξ1 + ρ cos θ


 ξ2 = ξ(s) + ρ sin θ
2

so that dξ1 dξ2 = ρdρdθ. Due to piecewise smooth property of the boundary of ∆, the triangle is split into three
sub-triangles. The associated integral of Eq. (12) now becomes
I = lim

ρ(ε)→0

3 Z
X
j=1

θj
θ j−1

Z

ρ̂(θ)

ρ(ε)

K(ρ, θ)φ(ρ, θ) |J(ρ, θ)| ρdρdθ,

(13)

where ρ̂(θ) gives a parametrization of the boundary of ∆ in polar coordinates, (θ j−1 , θ j ) are three intervals on which
ρ̂(θ) is smooth. The limiting process means CPV or HFP integral mentioned before, although for well-defined integral
operator S, the limiting is not necessary. According to the geometrical relationship in Fig.2
hj
.
cos θ̄

ρ̂(θ) =

(14)

where, h j is the perpendicular distance from ξ s to jth side of the planar triangle and θ̄ is the angle from the perpendicular to a point ξ (Fig.2). In each sub-triangle, θ̄ equals to θ minus a constant.
ξ2

θ̄

θ̄ ≈ 90◦

θ

ρ̂(θ)
h

ξs

ξs
ξ1
Figure 2: Integration under polar coordinates. Left: Polar coordinates transformation; Right: near singularity as the field point approaches the
boundary

For singularity of order no more than 3, the integrand of (13), denoted by F(ρ, θ), can be expressed as series
expansion under polar coordinate ([23, 8])
F(ρ, θ) = K(ρ, θ)φ(ρ, θ) |J| ρdρ =

∞

X
f−2 (θ) f−1 (θ)
2
+
f
(θ)
+
ρ
f
(θ)
+
ρ
f
(θ)
+
·
·
·
=
ρi fi (θ),
+
0
1
2
ρ
ρ2
i=p
5

(15)

where, fi are just functions of θ, integer p is determined by the order of singularity,


0,
weakly singular;



−1,
strongly singular;
p=


 −2, hyper-singular.

First, let’s consider the hypersingular integrals. Due to the appearance of the two terms with (i = −2, −1), the
integration of F(ρ, θ) must be performed in the HFP sense; for more details see Guiggiani’s work [8]. The resultant
formula for hypersingular integrals is given by
3 Z
X

I = I1 + I2 =

θ j−1

j=1

|

+

θj

Z

ρ̂(θ)

0

"

f−2 (θ) f−1 (θ)
F(ρ, θ) −
+
ρ
ρ2
{z

!#

dρdθ

I1

3 Z
X

θj

θ j−1

j=1

|

!
1
dθ .
f−1 (θ) ln ρ̂(θ) − f−2 (θ)
ρ̂(θ)
{z
}

}

(16)

I2

The strongly and weak singular integrals can be treated in a similar manner based on expansion (15). The resultant
computing formulas can also be written in the above form except that for strongly singular integrals f−2 = 0, and for
weak singular integrals f−2 = f−1 = 0 (thus I2 = 0). Therefore, formula (16) actually provide an unified approach to
evaluate singular integrals in BEM.
In (16) the double integral I1 and single integral I2 are all regular, thus one may conclude that it is sufficient to guarantee numerical accuracy in evaluating these two integrals, it would be nevertheless too expensive in practical usage,
especially when used in high order Nyström method considered in this paper. The difficulties and the corresponding
solutions will be presented in the following sections.
3.2. Difficulties in evaluating I1
It is the computation of I1 that accounts for the main cost for evaluating the various singular integrals. The
integrand of I1 can be approximated by polynomials in ρ; that is,
I1 =

3 Z
X
j=1

θj

θ j−1

Z

ρ̂(θ)

0




f0 (θ) + ρ f1 (θ) + ρ2 f2 (θ) + . . . dρdθ.

(17)

The number of terms in this approximation is determined by the kernel function, the order of basis function and the
flatness of the associated element. The relative flatness of element is a basic requirement in BEM in order to guarantee
the accuracy. Consequently, it appears that the integrand of I1 can be well approximated by low order polynomials
in ρ in solving many problems including Laplace, Helmholtz, elasticity and so forth, using quadratic elements. This
implies that low order Gaussian quadratures are sufficient for numerical integration in ρ.
In angular θ direction, however, two difficulties are frequently encountered which severely retard the convergence
rate of Gaussian quadratures. To show this, preforming integration in ρ in (17) one obtains
I1 =

3 Z
X
j=1

θj

θ j−1

!
1 2
ρ̂(θ) f0 (θ) + ρ̂ (θ) f1 (θ) + · · · dθ.
2

(18)

The difficulties in computation of I1 are caused by the near singularities of fi (θ) and ρ̂(θ).
3.2.1. Near singularity in fi (θ)
Functions fi (θ) can be expressed as (Appendix A)
fi (θ) =
6

f˜i (θ)
,
Aα (θ)

(19)

where f˜i are regular trigonometric functions and α is integer determined by the subscript i. Function A(θ) is depend
on the shape of the element and parametric coordinate system (see [8] for a definition of A(θ)). Specifically, let u1 , u2
be the two column vectors of the Jacobian matrix J which spans the space tangent to the element at the point x, i.e.
∂y
∂ξ1

u1 =

and u2 =
y=x

∂y
∂ξ2

.

(20)

y=x

Then A(θ) is given by
q
A(θ) = |u1 |2 cos2 θ + u1 · u2 sin 2θ + |u2 |2 sin2 θ
r


1
1
= u1 · u2 sin 2θ + |u1 |2 − |u2 |2 cos 2θ + |u1 |2 + |u2 |2
2
2
r




1
=
|u1 |2 + |u2 |2 µ sin(2θ + ϕ) + 1 .
2

If we let λ = |u1 | / |u2 | and cos γ =

u1 ·u2
|u1 ||u2 | ,

(21)

then

ϕ = arctan
and
µ=

s

1−

λ2 − 1
,
2λ cos γ

4 sin2 γ
< 1.
(λ + λ−1 )2

(22)

It can be seen from (21) that, if µ → 1, there exist two points θ ∈ [0, 2π] such that A(θ) → 0. Thus function fi
tends to be nearly singular. The circumstance µ → 1 occurs in two cases according to Eq. (22): (1) the aspect ratio
of the element is large, i.e., λ approaches 0 or ∞; (2) peak or big obtuse corner appear in the element which lead to
sin γ → 0. Both of these two cases indicate a distorted shape of the element.
To illustrate the influence of the element aspect ratio on the behavior of function 1/A(θ) which represents the
smoothness of fi , consider the curved element in figure 7; see section 5.1 for more detailed descriptions. The aspect
ratio of the element is controlled by s; larger s implies more distorted shape. Let b be the singular point. The plots of
1/A(θ) in the interval corresponding to the sub-triangle (2-b-3) in figure 7 with various s are exhibited in Fig. 3. It is
clear that 1/A(θ) varies more acutely as the increase of the aspect ratio.
The above analysis shows how element shape affects fi and thus the integrand in (18). In an intuitionistic manner,
the near singularity in fi is induced by the fact that the integration plane ∆ in intrinsic coordinate system is independent
on the shape of element, as a result the distortion is brought into the integrand. One can suppose that if the integration
is performed over another planar triangle which reflects the distortion of the element, the near singularity should be
eliminated. This is the main idea of the our new method in section 4.1.
3.2.2. Near singularity in ρ̂(θ)
In addition to the near singularity in fi , another obstacle retards the convergence of the numerical quadrature for
(18) is the near singularity in ρ̂(θ), which can be clearly seen from (14). When the field point x lies near to the
boundary of the element, the restrict of θ̄ approaches ±π/2, thus the denominator cos θ̄ is close to 0 at the two ends of
the interval (θ j−1 , θ j ) (see Fig.2). The effect of the near singularity in ρ̂(θ) on the total behavior of the integrand F(ρ, θ)
is demonstrated by the left plots of Fig. 6 (a) and (b), where the integrand have large peaks near the two ends of the
interval.
Unfortunately, the situation that causes the near singularity in ρ̂(θ) is ubiquitous in using high order Nyström
method. See Fig. 3.1 for a typical distribution of field points in 2th order Nyström method.
4. Efficient transformation methods
How to effectively resolve the above mentioned two difficulties is crucial to the accurate evaluation of the BEM
singular integrals. More recently, special quadratures are constructed for this purpose in [3]. Although accurate and
7

1.5
s=0.5
s=2.0
s=4.0
s=10

1/A(θ)

1

0.5

0

−0.5

0

0.5

θ (rad)

1

1.5

2

Figure 3: The plot of 1/A(θ) for various value of s. The aspect ratio of the associated element increased with s. The curved element is described in
section 5.1, with field point (b)

robust numerical results are reported, it is noticed that the abscissas and weights of the special quadratures are depend
on the singular (field) point as well as the element on which the integral is defined. The construction of the quadrature
for integral can be rather complicated and time-consuming.
In this section, however, we propose a more simple yet efficient method. First, the intrinsic coordinates (ξ1 , ξ2 )
is transformed onto a new system which results in constant A(θ). An additional benefit of constant A(θ) is that the
line integral I2 can be evaluated in closed form; see section 4.3. Then the sigmoidal transformation is introduced to
alleviate the near singularity caused by ρ̂(θ).
4.1. Conformal transformation
First, the near singularity caused by A(θ) in fi is considered and resolved. The idea is to introduce a new transformation under which A(θ) becomes constant. It can be seen from (21) that if
u1 · u2 = 0

and

|u1 | = |u2 | ,

(23)

then A(θ) will be constant, i.e. A(θ) = |u1 | = |u2 |.
Relations (23) indicate a conformal mapping from the element to the plane at field point x, i.e., both angle and
the shape of the infinitesimal neighborhood at x are preserved. However, this condition generally can not be satisfied
under intrinsic coordinates. In [18] Hayami proposed to connect the three corners of the curved element to establish a
planar triangle which preserves the shape of the element. This operation, nevertheless, can not satisfy condition (23)
exactly.
Here, we propose a transformation in which the curved physical element is mapped to a triangle ∆¯ in plane
(2)
(2)
(η1 , η2 ) as shown in Fig. 4. The coordinates of the three corners of ∆¯ are (0, 0), (η(2)
1 , η2 ) and (1, 0), with η2 > 0.
The transformation from system (ξ1 , ξ2 ) to (η1 , η2 ) can be realized by linear interpolation
η=

3
X

φ̄i (ξ1 , ξ2 )η(i) ,

(24)

i=1

¯ φ̄i are linear interpolating functions
where, η(i) are coordinates of the three corners of ∆,
φ̄1 = 1 − ξ1 − ξ2 ,
φ̄2 = ξ2 ,
φ̄3 = ξ1 .

8

(25)

ξ2
1
0.5

2

x3

5

T

1 ξ1
T−1

u1
6

1

(0.5, 0.5)

0.5

0

J

u2
4x

ξ (s)

(2)

(2)

(η1 , η2 )

η2

3

θ

x2
x1
o

η

(s)

η1

Figure 4: Transformations between coordinate systems. Left: curved element in global coordinates; Top right: intrinsic coordinates and the
reference element; Bottom right: the parametric plane to perform integral

Plugging (25) into (24) yields
η = Tξ,

(26)

where, T is the transformation matrix
"
#
1 η(2)
1
T=
0 η(2)
2

−1

and T



1 − η(2)
1 

(2) 

 .
η2 
= 

1
0
(2)
η

(27)

2

The Jacobian matrix at x from physical coordinates to (η1 , η2 ) coordinates, denoted by [ū1 ū2 ], can be written as


i
i h
h
−η(2)
u2
1 u1
ū1 ū2 = u1 u2 T−1 = u1
(28)
+
(2)
(2) .
η
η
2

2

In order to obtain constant A(θ), ū1 and ū2 have to satisfy condition (23), thus one has
η(2)
1 =

cos γ
,
λ

η(2)
2 =

sin γ
,
λ

(29)

where, λ and γ are defined in section 3.2.
By using relation (26) the integral (12) can be transformed onto (η1 , η2 ) plane. One can then employ the polar
coordinate transform in Section 3.1. The origin of the polar system is set as the image of x on (η1 , η2 ) plane. Then
integral I1 becomes
!#
3 Z θ j Z ρ̂(θ) "
X
f−2 (θ) f−1 (θ)
T−1 ρdρdθ,
(30)
F(ρ, θ) −
I1 =
+
2
ρ
ρ
j=1 θ j−1 0
in which the near singularity caused by A(θ) has been successfully removed.
Numerical results show that the above transformation can always improve the numerical integration regardless the
shapes of the elements being regular or irregular.
4.2. Sigmoidal transformation
When the field point x approaches the element boundary, function ρ̂(θ) and thus the integral I1 becomes nearly
singular. Here a sigmoidal transformation is introduced to alleviate this problem. The sigmoidal transformation
was first proposed to calculate the singular and nearly singular integrals in two dimensional BEM [22]. Recently,
9

1
0.5

1
Gaussian points
clustered points

Gaussian points
clustered points

0.5

θ̄

0

θ̄

0
−0.5

−0.5

−1

−1

0

0.2

0.4

0.6

0.8

1

0

0.2

0.4

w

0.6

0.8

1

w

Figure 5: Sigmoidal transformation with m = 3 to cluster 10 Gaussian quadrature points (x-axis) to the endpoints of the interval. The sub-triangle
(2-b-3) in Fig. 7 is considered, whose θ̄ ∈ [−1.466, 1.005]. Left: A naive use of sigmoidal transformation; Right: Sigmoidal transformation in this
paper.

this approach was adopted as an angular transformation in dealing with nearly singular integrals in 3D BEM [19].
One should notice that the angular transformation proposed by Khayat and Wilton [13] can be alternatively used,
but our numerical experience indicates a better overall performance by using sigmoidal transformation, especially in
hypersingular case.
A sigmoidal transformation can be thought of as a mapping of the interval [0, 1] onto itself whose graph is Sshaped. It has the effect of translating a grid of evenly spaced points on [0, 1] onto a non-uniform grid with the node
points clustered at the endpoints. A typical sigmoidal transformation is given by [22]
σ(w) =

wm

wm
,
+ (1 − w)m

σ, w ∈ [0, 1], m ≥ 1.

(31)

Consider one sub-triangle in figure 2 in which θ̄ ∈ [θ̄ j−1 , θ̄ j ], j = 1, 2, 3. A naive use of the above transformation can
be (see Fig. 5)
θ̄ − θ̄ j−1
= σ(w).
(32)
θ̄ j − θ̄ j−1

However, this would be of low efficiency. Since the integral I1 tends to be nearly singular (due to ρ̂(θ)) only when
θ̄ j−1 or θ̄ j is close to right angle or both, it is thus more reasonable to cluster the quadrature nodes according to the
discrepancy of θ̄ j−1 and θ̄ j to right angles, respectively. In addition there are cases where both θ̄ j−1 and θ̄ j are not very
close to right angles and thus the near singularity is not severe. For these cases, one can use the Gauss quadratures
directly without any transformation. Nevertheless, numerical examples in this paper show that the modification put
forward below can always achieve more accurate results.
The transformation used in this paper is based on the fact that the singularity occurs at θ̄ = ±π/2 and is given by
1
π
θ̄ +
= σ(z), z = (z j − z j−1 )w + z j−1 , w ∈ [0, 1],
(33)
π
2

where, z(θ̄ j−1 ) and z(θ̄ j ) are the values of z corresponding to θ̄ j−1 and θ̄ j which can be easily obtained by (31). Obviously, both transformation (33) and (32) project [θ̄ j−1 , θ̄ j ] onto [0, 1], but (33) can cluster the quadrature points
adaptively according to closeness of θ̄ to right angle, which is shown in Fig. 5. The left hand image demonstrate the
part of sigmoidal transformation used in dealing with the sub-triangle (2-b-3) in Fig. 7. Compared with full sigmoidal
transformation shown left, the distribution of quadrature points are more reasonable since they are more clustered on
θ̄ j−1 which is more closely to right angle. Besides, the points are not too away from the center. A comparison between
the original sigmoidal transformation and the modification is illustrated in table 1. Without any doubt, the modified
transformation is more powerful when dealing with the integrals in this paper.
The degree of cluster is impacted by parameter m, see [22]. Although numerical examples indicates not too much
difference in number of quadrature points when m takes the value between 2 and 3, it should be point out that the
optimum value of m is 3 for weakly singular case and 2 for hypesingular case.
10

Table 1: Comparison of relative errors using different methods in k = 0 case. The curved element is show in Fig. 7, with field point (b). “n” denotes
the number of quadrature points used in angular direction.

n
5
8
10

single
Eq. (32)
3.19e-3
1.16e-5
3.06e-6

Guiggiani
6.24e-4
1.02e-3
3.20e-4

0.8

Present
9.83e-5
6.44e-7
1.26e-8

Present
4.09e-05
1.49e-07
3.71e-09

0.3

integand
8th degree

integrand
8th degree

0.25

0.6

0.2
integrand

integrand

hyper
Eq. (32)
1.57e-04
1.87e-06
5.22e-07

Guiggiani
1.80e-03
1.69e-04
2.39e-05

Norm of residuals = 0.21819

0.4

0.15
Norm of residuals = 0.00738

0.1
0.05

0.2

0
0

-0.5

0

0.5

1

T (rad)

1.5

-0.05

2

-0.5

0

0.5

T (rad)

1

1.5

2

(a) Laplace single kernel, m = 3
4

integrand
8th degree

integrand
8th degree

0.6
0.5

2

integrand

integrand

3
Norm of residuals = 0.83882

1

Norm of residuals = 0.00946

0.4
0.3
0.2
0.1

0

0
-1

-0.5

0

0.5

1

T (rad)

1.5

-0.1

2

-0.5

0

0.5

T (rad)

1

1.5

2

(b) Laplace hypersingular kernel, m = 2
Figure 6: Integrand of I1 in angular direction corresponding to sub-triangle (2-b-3) in Fig. 7;the basis function is selected as ξ22 . 8th order
polynomial is plotted to fit the image. Left: Image of the original Guiggiani’s method; Right: Image after transformations proposed in this paper.

Finally, the integral has the formula
I1 =

3 Z
X
j=1

0

1

Z

ρ̂(θ)

0

"

F(ρ, w) −

f−2 (w) f−1 (w)
+
ρ
ρ2

!#

θ̄′j (w) T−1 ρdρdw,

(34)

where, θ̄ j (w) denotes transformation (33) on interval [θ̄ j−1 , θ̄ j ].
The effects of the two transformations in sub-sections 4.1 and 4.2 are demonstrated in Fig. 6. It can be seen that
after the transformations, the integrands become more regular and can be well approximated by polynomials of order
8, which means that low order Gaussian quadratures can achieve accurate results.
4.3. Analytical integration of the one-dimensional integrals
P
Consider the evaluation of the regular one dimensional integral I2 in Eq. (16), i.e., I2 = 3j=1 I2( j) and
( j)

I2 =

Z

θj

θ j−1

f−1 (θ) ln ρ̂(θ) − f−2 (θ)

!
1
dθ.
ρ̂(θ)

(35)

One notes that the integrand of I2 suffers from the same difficulties as that of the double integral I1 as stated
in section 3.2. The above proposed transformations can also be used to promote the efficiency of the numerical
quadratures.
11

In this section close form expression for I2 is derived. It is attributed to the fact that after the conformal transformation in section 4.1 the denominators of fi , i.e. A(θ), become constant. For Burton-Miller equation, it is readily
to see from (A.7) that f−2 (θ) is constant on each sub-triangle thanks to the constant A(θ), and that f−1 (θ) is typically
made up of combinations of elementary trigonometric functions of θ [8]
f−1 (θ) =c1 cos3 θ + c2 cos2 θ sin θ + c3 cos θ sin2 θ + c4 sin3 θ
+ d1 cos θ + d2 sin θ.

(36)

The coefficients ci and di are given in (A.8). Here we notice that the expressions in (A.8) are also applicable for Laplace
equations. For other problems, such as elasticity, Stokes flow, the expressions of f−2 and f−1 are more complicated
and can be derived analogously as in Appendix A.
To further simplify the derivation, the coordinate system (η1 , η2 ) is rotated so that the axis oη1 is parallel to the
¯ Thus θ = θ̄ over each sub-triangle. By substituting ρ̂(θ), f−2 and f−1 into (35), the
perpendicular line of the side of ∆.
integral with f−2 becomes
Z θj
Z
f ( j) θ j
1
f−2 (θ)
dθ = −2
cos θdθ.
ρ̂(θ)
h j θ j−1
θ j−1
The integral with f−1 consists of terms
ln h j

Z

θj

cos p θ sinq θdθ,

p+q=3

or

p + q = 1,

(37)

θ j−1

and

Z

I¯i =

θj

cos p θ sinq θ ln cos θdt,

p+q=3

or

p + q = 1.

(38)

θ j−1

The explicit expressions of integral (37) can be easily obtained and are omitted. Here the expression of I¯i is derived.
By changing the integral variable
t = sin θ
√
cos θ = 1 − t2
(39)
√
dt = 1 − t2 dθ,
one obtains
I¯1 =

Z

θj

3

cos θ ln cos θdθ =

θ j−1

Z

tj

t j−1

√
√
t 1 − t2 ln 1 − t2 dt = I˜1 − I˜3 ,

√
t − t3
ln 1 − t2 dt = I˜2 − I˜4 ,
√
θ j−1
t j−1
1 − t2
Z tj
Z θj
√
t2 ln 1 − t2 dt = I˜3 ,
cos θ sin2 θ ln cos θdθ =
I¯3 =

I¯2 =

Z

I¯4 =

Z

θj

cos2 θ sin θ ln cos θdθ =

θ j−1
θj

Z

sin3 θ ln cos θdθ =

I¯6 =

θj

cos θ ln cos θdθ =
θj

θ j−1

Z

tj

t j−1

θ j−1

Z

tj

t j−1

θ j−1

I¯5 =

Z

sin θ ln cos θdθ =

Z

Z

tj

t j−1
tj

t j−1

t3

ln

√

1 − t2 dt = I˜4 ,

√
1 − t2
√
ln 1 − t2 dt = I˜1 ,

√
t
ln 1 − t2 dt = I˜2 ,
√
1 − t2

12

(40)

where,

t
√
 √

t+1 j
,
1 − t2 dt = t ln 1 − t2 − 1 + ln √
t j−1
1 − t2 t j−1
√
Z tj
√

√
 tj
t ln 1 − t2
dt = 1 − t2 1 − ln 1 − t2
=
,
√
t j−1
t j−1
1 − t2
Z tj
 1
√
√
1 
1 1 + t tj
,
t2 ln 1 − t2 dt = t3 3 ln 1 − t2 − 1 − t + ln
=
9
3
6 1 − t t j−1
t j−1
"
# tj
Z tj 3 √
√
√
1
1
t ln 1 − t2
2
2
2
2
(8 + t ) − (2 + t ) ln 1 − t
dt = 1 − t
.
=
√
9
3
t j−1
1 − t2
t j−1

I˜1 =
I˜2
I˜3
I˜4

Z

tj

ln

(41)

5. Numerical examples

Figure 7: Curved element and the field points in example 1

Two numerical examples are provided to demonstrate the efficiency and robustness of the method in this paper.
The first one is used to illustrate the superior of the method in this paper to the method in 3.1, i.e., the classic polar
coordinates transformation and Guiggiani’s method [8], using a sample curved element. The second one is used to
verify the accuracy and convergence of the method when employed in solving Burton-Miller equation. The method is
implemented in C language. The numerical integration code are freely available.
The focus of this paper is to improve the performance of the quadrature in angular direction, while in radial
direction, as explained in section 3.2, the integrand behaves well so that only a small number of quadrature points
are sufficient. In section 5.1 6 points are used in every test in radial direction. However, when implemented in BEM
programs, fewer points are needed because as the subdivision of the model the diameter of each element will decrease
so that fewer terms are needed in expansion (15) to approximate the integrand. As a result, the rate of convergence of
BEM will not be affected. For the 2th order Nyström BEM in section 5.2, we use 3 points in polar direction.
5.1. Accuracy of the method on a curved element
The accuracy, effectiveness and robustness of the method in this paper are investigated. The integrals is performed
over a curved element extracted from a cylinder surface. The size of the cylinder and the element are given in Fig.7.
Distance s in Fig. 7 is designed to change the aspect ratio of the element. s = 0.5 represents a quite regular element;
the quality of the element becomes bad as the increase of s. Four positions for the field point x are considered, which
correspond to intrinsic coordinates (a) ξ = (0.3, 0.3), (b) ξ = (0.1, 0.8), (c) ξ = (0.45, 0.45), (d) ξ = (0.64, 0.31),
respectively. One should note that points (b) and (c) match approximately with the nodes of the 2th order Nyström
13

Table 2: methods to compare

methods
Guiggiani
Gui+sig
Present
Present+a

description
Guiggiani’s method in section 3.1
Guiggiani’s method with a sigmoidal transformation applied in angular direciton
Method in this paper, i.e. Guiggiani’s method with the two transformations in section 4
Similar to “Present”, except that the line integral is calculated analytically

discretization, while (d) match with that of the 3th order. In addition, for the rather centered point (a) the original
Guiggiani’s method in section 3.1 are believed to work well.
Since no analytical solutions are available, the (relative) error, given by
relative error =

|Icalc − Iref |
|Iref |

is evaluated by comparison with the results Iref of the Guiggiani’s methods using Gauss quadratures with 256 quadrature points in angular direction and 32 points in radial direction. Icalc denotes the result calculated by the various
methods. The single layer operator S and the hypersingular operator H are tested. The parameter m in sigmoidal
transformation takes the value of 3 in weakly singular case and 2 in hypersingular case. The basis function φ is
selected as the second order monomial ξ22 .
A description of the methods tested are listed in table 2. The method in this paper are indicated by “Present” in
which the line integral I2 in (16) is computed by Gauss quadratures; while in “Present+a” the line integral is computed
by using the closed formulations in section 4.3. Note that “Present+a” only validates for hypersingular integrals.
5.1.1. Results for regular element
First the performance of the present method for regular element are verified; for this purpose, we let s = 0.5 and
the wave number k = 0 (Laplace kernels).
Convergence behavior of the original methods and the present method for the four field points are illustrated in
Fig.8. The results of the present method are marked by triangles or circles. The results for single layer kernel are
plotted with dashed line, while real lines are for hypersingular integrals. A comparison between case (a) and the other
three cases shows that the original Guiggiani’s method tends to converge slowly when the field point is closed to the
element boundary; while the present method can achieve almost the same fast convergence. Even in regular case (a)
in which the original method performs better than other three cases, the benefit of the present method is also obvious.
The effect of each transformation scheme in section 4 are studied and exhibited in table 3. Both field points (b)
and (d) are tested. For the weak singular (single layer) integral, the accuracy can be improved by only using the
sigmoidal transformation; whereas the present method, say the combination of the two transformations in section 4,
can further promote the accuracy significantly. One notice that for the hypersingular case the sole use of the sigmoidal
transformation can not improve the accuracy. However, this problem can be completely overcome by using the
parametric coordinates transformation. The analytical formulas for the line integrals can be better than the numerical
quadratures.
Now consider the Helmholtz kernels. The results are similar to that of the Laplace kernels, so only the results for
field point (b) are given in table 4. The wavenumber k is selected as 2.0 so that the wavelength is about 3 times size of
the triangular element. In practical BEM implementation such a mesh is rather coarse. It is obvious that the present
method can greatly improve the accuracy.
5.1.2. Results for distorted element
As have been mentioned in section 4.1, the present method should be robust for elements with high aspect ratio;
here this is validated. The integrals are evaluated over elements with different s and k = 0. Point (a) is taken to be
the field point. Both the original and present methods are used to compute the integrals with relative error 10−8 . The
needed number of Gauss quadrature points in angular direction are listed in table 5. It is seen that this number of the
original method increases drastically as the increase of s, which clearly indicates that the performance of the original
method depend heavily on the shape of the element. On the contrary, the present method behaves more robust with the
14

0

0

10

10
Guiggiani (single)
Present (single)
Guiggiani (hyper)
Present+a (hyper)

−2

10

10

−4

10

relative error

relative error

−4

−6

10

−8

10

−10

−6

10

−8

10
10

−12

2

10

−10

10
10

Guiggiani (single)
Present (single)
Guiggiani (hyper)
Present+a (hyper)

−2

−12

4
6
8
10
12
14
16
18
number of integration points in angular direction

10

20

2

4
6
8
10
12
14
16
18
number of integration points in angular direction

(a) field point (a)

(b) field point (b)

0

0

10

10
Guiggiani (single)
Present (single)
Guiggiani (hyper)
Present+a (hyper)

−2

10

−2

10

−4

relative error

relative error

−4

10

−6

10

−8

10

−10

−6

10

−8

10
10

−12

2

10

−10

10
10

20

Guiggiani (single)
Present (single)
Guiggiani (hyper)
Present+a (hyper)

−12

4
6
8
10
12
14
16
18
number of integration points in angular direction

10

20

2

4
6
8
10
12
14
16
18
number of integration points in angular direction

(c) field point (c)

(d) field point (d)

Figure 8: Convergence comparison for the k = 0 case.
Table 3: Relative errors for different methods. “n” denotes the number of Gauss quadrature points in angular direction.

n

Guiggiani

single
Gui+sig

Present

5
8
10
12

6.24e-04
1.02e-03
3.20e-04
5.65e-05

1.22e-04
2.83e-05
3.08e-07
8.91e-08

9.83e-05
6.44e-07
1.26e-08
2.95e-11

5
8
10
12

4.05e-02
1.89e-02
1.07e-02
5.83e-03

1.18e-03
3.20e-05
5.85e-07
1.50e-08

6.06e-04
6.35e-06
2.23e-07
5.16e-09

Guiggiani
field point (b)
1.81e-03
1.70e-04
2.38e-05
3.01e-06
field point (d)
5.33e-03
2.35e-03
8.11e-04
3.45e-04

hyper
Gui+sig
Present

Present+a

5.51e-03
2.62e-04
2.72e-05
1.12e-05

1.28e-03
5.41e-06
5.64e-08
4.20e-10

4.1e-05
1.49e-07
3.71e-09
6.70e-10

6.28e-03
1.48e-04
2.65e-04
3.15e-05

1.11e-02
7.32e-05
1.93e-06
5.34e-08

1.68e-04
4.08e-06
2.92e-07
1.57e-08

Table 4: Relative errors in dynamic case with k = 2.0, the field point is selected as (b).

n
5
8
10
12

single
Guiggiani Present
2.52e-03 2.00e-04
1.16e-03 2.44e-06
4.09e-04 1.09e-07
1.26e-04 4.66e-09
15

hyper
Guiggiani Present+a
1.92e-03
3.64e-04
2.08e-04
3.22e-07
2.04e-05
6.93e-09
4.36e-06
4.71e-10

20

Table 5: Number of integration points needed in angular direction to make the relative error under 10−8

s
0.5
1.5
2.0
4.0
10.0

single
Guiggiani Present
15
8
19
9
22
9
42
11
82
13

Guiggiani
18
27
31
62
142

hyper
Present
11
12
13
14
15

Present+a
11
11
12
13
16

Figure 9: Geometry model

change of s. For element with aspect ratio larger than 10, high accurate results can still be achieved with a moderate
increase of quadrature points.
5.2. An exterior sound radiation problems
The overall performance of the present integration method is tested by incorporating it into a Nyström BEM code
for solving the Burton-Miller equation. The geometry of the problem consists of three sections of cylinders, as shown
in Fig. 9. The origin of coordinate system lies on the middle point of the symmetry axis of the cylinder. All the
surfaces of the cylinder vibrate with given velocity q(x) = ∂u/∂x (Neumann problem), where u is chosen to be the
fundamental solution u(x) = G(x, y) in (3) with y = (0, 0, 0).
The Nyström BEM with 2th order basis functions is used to discretize Burton-Miller equation. The wave number
k = 2.5. The surface is meshed by quadratic triangular element and refined four times; the finest mesh consists of
1638 elements. In evaluating the singular integrals, both the original Guiggiani’s method and the present method are
tested. A common m = 2.5 is used in the sigmoidal transformation for all the singular integrals, and 3-point Gauss
quadrature are used in radial direction.
Three different quadratures in the angular direction are used and compared, that is, the present method using 4
points and the Guiggiani’s method with 10 and 15 points. The L2 relative errors of the boundary values of u(x) are
demonstrated in Fig. 10. It is shown that the accuracy of BEM is retarded when using the Guiggiani’s method with
10 points. For the Guiggiani’s method with 15 points, the convergence of the BEM tends to slow down with further
refinement of the mesh. The present method with 4 points, however, can achieve the theoretical converge rate.
6. Conclusions
Highly efficient methods with high accuracy and low computational cost are crucial for high-order boundary
element analysis. In [8], Guiggiani proposed an unified framework to treat the singular integrals of various orders
in BEM. It is based on the polar coordinate transformation which has been extensively used in dealing with BEM
singular integrals. However, the performance of the polar coordinate transformation deteriorates when the field point
is close to the element boundary or the aspect ratio of the element becomes large.
16

−2

O(N )
Present(n=4)
Guiggiani(n=10)
Guiggiani(n=15)

−3

2

L error

10

−4

10

−5

10

3

4

10

10
Number of unkowns
Figure 10: L2 -errors of u in example 5.2.

In this paper, first, a conformal transformation is introduced to circumvent the near singularity caused by large
aspect ratio of element. This transformation maps a curved physical element onto a planar triangle. Since it is
conformal at the field point, the resultant integration domain (planar triangle) perseveres the shape of the curved
element. Then, a sigmoidal transformation is applied to alleviate the near singularity due to the closeness of the
field point to the element boundary. The rationale behind this is that the sigmoidal transformation can cluster the
quadrature points to area of near singularity. The combination of the two transformations can effectively alleviate
the two problems of the existing polar coordinate transformation method, and thus leads to considerable reduction of
quadrature points in angular direction.
The efficiency and robustness of the present method are illustrated by various singular integrals on a curved
quadratic element which is typical in high order Nyström BEM. It is shown that highly accurate results with relative
error 10−8 can be achieved with 10 quadrature points in angular direction. Moreover, the method is more stable with
the change of element aspect ratio. For further verification, the method is applied to a 2-order Nyström BEM for
solving acoustic Burton-Miller equation. Theoretical convergence rate of the BEM can be retained with much less
quadrature points than the existing quadrature methods.
Acknowledgements
This work was supported by National Science Foundations of China under Grants 11074201 and 11102154 and
Funds for Doctor Station from the Chinese Ministry of Education under Grants 20106102120009 and 20116102110006.
Appendix A. Coefficients in Eq. (36)
Consider a general integrand in BEM in polar coordinates
F(ρ, θ) =

ρF̄(ρ, θ)
,
rβ

(A.1)

where, β is the order of singularity, F̄(ρ, θ) is a regular function. There holds
∞

X
1
= ρ−β
S ν−β (θ)ρν ,
β
r
ν=0
17

(A.2)

F̄(ρ, θ) =

∞
X

av (θ)ρν ,

(A.3)

ν=0

where, S ν (θ) is given by ([23], Theorem 4)

S ν−β (θ) = A−β−2ν (θ)g3ν (θ),

(A.4)

g3ν (θ) are homogeneous trigonometric polynomials of order 3v.
For hypersingular integrand, β = 3, thus,
f−2 (θ) = S −3 (θ)a0 (θ),
f−1 (θ) = S −2 (θ)a0 (θ) + S −3 (θ)a1 (θ).

(A.5)

See [8] for the expressions of S −2 (θ) and S −3 (θ).
For hypersingular kernel of Burton-Miller equation in this paper, the expressions ai can be obtained similarly with
[8]. Let Jk be the kth component of the vector
∂y
∂y
×
.
∂η1 ∂η2
Then Jk can be expanded as
#
"
∂Jk
∂Jk
cos θ +
sin θ = Jk0 + ρJk1 (θ) + O(ρ2 ).
Jk = Jk0 + ρ
∂η1 η=ηs
∂η2 η=ηs
The basis function φ can be expanded analogously as
"
∂φ
∂φ
φ = φ0 + ρ
cos θ +
∂η1 η=ηs
∂η2

#

sin θ = φ0 + ρφ1 (θ) + O(ρ2 ).
η=ηs

Then, a0 and a1 (θ) in Eq. (A.5) can be expressed as (repeated indicies imply summation)
ni (x)
Ji0 φ0 ,
4π

ni (x) 
Ji1 (θ)φ0 + Ji0 φ1 (θ) .
a1 (θ) =
4π
Substituting Eq. (A.6) into Eq. (A.5), yields
a0 =

(A.6)

f−1 (θ) =c1 cos3 θ + c2 cos2 θ sin θ + c3 cos θ sin2 θ + c4 sin3 θ
+ d1 cos θ + d2 sin θ,
ni (x)
f−2 (θ) =
Ji0 φ0 .
4πA3
with

3ni (x)Ji0 φ0 ∂yk ∂2 yk
,
8πA5 ∂η1 ∂η21


3ni (x)Ji0 φ0  ∂yk ∂2 yk
1 ∂yk ∂2 yk 
 ,

c2 = −
+

∂η1 ∂η1 ∂η2 2 ∂η2 ∂η21
4πA5


1 ∂yk ∂2 yk 
3ni (x)Ji0 φ0  ∂yk ∂2 yk
 ,
+
c3 = −

∂η2 ∂η1 ∂η2 2 ∂η1 ∂η22
4πA5

(A.7)

c1 = −

3ni (x)Ji0 φ0 ∂yk ∂2 yk
c4 = −
,
8πA5 ∂η2 ∂η22

!
∂φ
ni (x) ∂Ji
,
φ
+
J
0
i0
∂η1
4πA3 ∂η1
!
∂φ
ni (x) ∂Ji
.
φ
+
J
d2 =
0
i0
∂η2
4πA3 ∂η2
Note that all the derivatives are evaluated at the field point x.
d1 =

18

(A.8)

References
[1] L. F. Canino, J. J. Ottusch, M. A. Stalzer, et al. Numerical solution of the Helmhotz equation in 2D and 3D using a high-order Nyström
discretization. Journal of computational physics, 146 (1998), 627-633.
[2] S. Rjasanow, L. Weggler. ACA accelerated high order BEM for maxwell probelms. Computational mechanics, 51 (2013), 431-441.
[3] J. Bremer, Z. Gimbutas. A Nyström method for weakly singular integral operators on surfaces. Journal of computational physics, 231 (2012),
4885-4903.
[4] X. W. Gao, T. G. Davies. Boundary Element Programming in Mechanics. Cambridge University Press (ISBN: 052177359-8), 2002.
[5] D. J. Willis, J. Peraire, J. K. White. A quadratic basis function, quadratic geometry, high order panel method. In 44th AIAA Aerospace sciences
meeting, number AIAA-2006-1253, 2006.
[6] H. J. Wu, Y. J. Liu, W. K. Jiang. A low frequency fast multiploe boundary element method based on analytical integration of the hypersingular
integral for 3D acoustic problems. Engineering analysis with boundary elements, 37 (2013), 309-318.
[7] S. N. Fata. Explicit expressions for 3D boundary integrals in potential theory. International journal for numerical methods in engineering, 78
(2009), 32-47.
[8] M. Guiggiani, G. Krishnasamy, T. J. Rudolphi, F. J. Rizzo. A general algorithm for the numerical solution of hypersingular boundary integral
equations. ASME Journal of applied mechanics, 59 (1992), 604-614.
[9] S. Järvenpää, M. Taskinen, P. Ylä-Oijala. Singularity extraction technique for integral equation methods with higher order basis functions on
plane triangles and tetrahedra. International journal for numerical methods in engineering, 58 (2003), 1149-1165.
[10] P. Kolm, V. Rokhlin. Numerical quadratures for singular and hypersingular integrals. Computers and Mathematics with Applications, 41
(2001), 327-352.
[11] R. D. Graglia. G. Lombardi. Machine precision evaluation of singular and nearly singular potential integrals by use of Gauss quadrature
formulas for rational functions. IEEE transactions on antennas and propagation, 56 (2008), 981-998.
[12] M. Carley. Numeriacl quadratures for singular and hypersingular integrals in boundary element methods. SIAM journal on scientific computing, 29 (2007), 1207-1216.
[13] M. A. Khayat, D. R. Wilton, P. W. Fink. An improved transformation and optimized sampling scheme for the numberical evaluation of
singular and near-singular potentials. IEEE transactions on antennas and propagation, 7 (2008), 377-380.
[14] M. A. Khayat, D. R. Wilton. Numerical evaluation of singular and near-singular potential integrals. IEEE transactions on antennas and
propagation, 53 (2005), 3180-3190.
[15] M. G. Duffy. Quadrature over a pyramid or cube of integrands with a singularity at a vertex. SIAM Journal on Numerical Analysis, 19(6):
1260-1262, 1982.
[16] B. M. Johnston, P. R. Johnston. A comparsion of transformation methods for evaluating two-dimensional weakly singular integrals. International journal for numerical methods in engineering, 56 (2003), 589-607.
[17] L. Scuderi. On the compuation of nearly singular integrals in 3D BEM collocation. International Journal for numerical methods in engineering, 74 (2008), 1733-1770.
[18] K. Hayami. Variable transformation for nearly singular integrals in the boundary element method. Research institute for mathematical sciences, Kyoto university, 41 (2005), 821-842.
[19] B. M. Johnston, P. R. Johnston, D. Elliott. A new method for the numerical evaluation of nearly singular integrals on triangular elements in
the 3D boundary element method. Journal of Computational and Applied Mathematics, 245 (2013), 148-161.
[20] A. J. Burton, G. F. Miller. The application of integral equation methods to the numerical solution of some exterior boundary-value problems.
Proceedings of the Royal Society of London, Series A, 323 (1971), 201-210.
[21] X. W. Gao. An effective method for numerical evaluation of gneral 2D and 3D high order singular boundary integrals. Computer methods in
applied mechanics and engineering, 199 (2010), 2856-2864.
[22] P. R. Johnston. Application of sigmoidal transformations to weakly singular and near-singular boundary element integrals. International
journal for numerical methods in engineering, 45 (1999), 1333-1348.
[23] C. Schwab, W. L. Wendland. Kernel properties and representations of boundary integral operators. Mathematische Nathrichten, 156 (1992),
187-218.

19

