10 August 2022
POLITECNICO DI TORINO
Repository ISTITUZIONALE
Machine Precision Evaluation of Singular and Nearly Singular Potential Integrals by Use of Gauss Quadrature Formulas
for Rational Functions / Graglia, Roberto; Lombardi, Guido. - In: IEEE TRANSACTIONS ON ANTENNAS AND
PROPAGATION. - ISSN 0018-926X. - STAMPA. - 56:4(2008), pp. 981-998. [10.1109/TAP.2008.919181]
Original
Machine Precision Evaluation of Singular and Nearly Singular Potential Integrals by Use of Gauss
Quadrature Formulas for Rational Functions
Publisher:
Published
DOI:10.1109/TAP.2008.919181
Terms of use:
openAccess
Publisher copyright
(Article begins on next page)
This article is made available under terms and conditions as specified in the corresponding bibliographic description in
the repository
Availability:
This version is available at: 11583/1513279 since:
IEEE
IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. 56, NO. 4, APRIL 2008
981
Machine Precision Evaluation of Singular and
Nearly Singular Potential Integrals by Use of Gauss
Quadrature Formulas for Rational Functions
Roberto D. Graglia, Fellow, IEEE, and Guido Lombardi, Member, IEEE
Abstract—A new technique for machine precision evaluation of
singular and nearly singular potential integrals with
1
singular-
ities is presented. The numerical quadrature scheme is based on a
new rational expression for the integrands, obtained by a cancella-
tion procedure. In particular, by using library routines for Gauss
quadrature of rational functions readily available in the literature,
this new expression permits the exact numerical integration of sin-
gular static potentials associated with polynomial source distribu-
tions. The rules to achieve the desired numerical accuracy for sin-
gular and nearly singular static and dynamic potential integrals
are presented and discussed, and several numerical examples are
provided.
Index Terms—Boundary element methods, finite element
methods, finite volume methods, integral equations, method of
moment (MoM).
I. INTRODUCTION
P
OTENTIAL integrals with unbounded singular kernels
arise in the moment-method solution of the integral equa-
tions of electromagnetism whenever the source and the testing
subdomain coincide, that is, in the so-called self-term case. The
integrals become nearly singular if the source and testing subdo-
mains are very close to each other, but do not overlap. All these
integrals are usually computed by the singularity cancellation
[1]–[9] or the singularity subtraction method [10]–[19], with
the near-singular case particularly considered in [1], [6]–[9],
[14], [20]–[22]. The cancellation method is based on variable
transformations whose Jacobian cancels out the singularity of
the integral kernel. The superiority of the cancellation method
with respect to the subtraction method is established in [1].
Despite its convenience, however, the technique in [1] does not
permit one to anticipate the precision of the numerical results,
even in the simpler case of static potential integrals, and does
not discuss reasons for the higher degree of difficulty inherent
in the near-singular case.
Due to the reduced cost and increased computation speed of
modern computers, the time is now ripe for machine precision
Manuscript received December 5, 2006; revised September 28, 2007. This
work was supported by NATO in the framework of the Science for Peace Pro-
gramme under the grant CBP.MD.SFPP 982376—Electromagnetic Signature
of Edge-Structures for Unexploded Ordnance Detection.
R. D. Graglia is with the Dipartimento di Elettronica, Politecnico di Torino,
10129 Torino, Italy and also with the ISMB-Instituto Superiore Mario Boella,
10138 Torino, Italy (e-mail: roberto.graglia@polito.it).
G. Lombardi is with the Dipartimento di Elettronica, Politecnico di Torino,
10129 Torino, Italy (e-mail: guido.lombardi@polito.it).
Digital Object Identifier 10.1109/TAP.2008.919181
evaluation of potential integrals. This paper describes a new nu-
merical technique to compute singular and nearly singular po-
tential integrals with machine precision; above all, it explains
how to obtain it and why it is not possible to predict the preci-
sion of the other existing integration techniques.
The most recent advances in the area of computational elec-
tromagnetics are not adequately represented in many of the pre-
viously quoted papers because of their age. In particular, the
majority of those papers prefer to redefine the integration do-
main of the potential integral by rotating and translating the
original given domain. Conversely, in this paper, the singularity
cancellation is carried out directly by working in the
parent ref-
erence-frame from the outset, thereby using the language and
the techniques well known to finite element analysts [23]. With
respect to the previous approaches, this has some important ad-
vantages that are described in Sections II and III. Section IV then
considers potential integrals on surface elements, explains the
new integration technique and elucidates the difficulties in eval-
uating near-singular integrals as compared to the singular ones.
Our quadrature scheme is based on a new rational representation
of the integrands resulting from the presented cancellation pro-
cedure. This expression permits the exact numerical integration
of the static potentials associated with polynomial source distri-
butions using library routines for Gauss quadrature of rational
functions that are readily available in the literature [24]. Inte-
grals on three-dimensional source domains are briefly discussed
in Section V, whereas Section VI provides several numerical re-
sults for potential integrals on surface elements. Section VI also
explains, in considerable detail and with several examples, how
to achieve machine precision accuracy for static and dynamic
potential integrals. Preliminary results of this work have been
reported in [25], [26].
II. O
BJECT AND
PARENT ELEMENTS FOR COMPUTATIONAL
ELECTROMAGNETICS
Modern electromagnetic (EM) codes model the geometry of
a given problem as the union of subdomains of different but
simple geometrical shapes. EM problems are then numerically
solved by expanding the unknowns in terms of vector or scalar
functions locally defined on these subdomains. The expansion
functions are conveniently defined on rectilinear domains of a
parent space, with all subdomains of the global geometry ob-
tained by properly mapping one or a few parent domains onto
the global object-space. In its parent space a domain is described
in terms of a set
of normalized (dimensionless) parametric
0018-926X/$25.00 © 2008 IEEE
982 IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. 56, NO. 4, APRIL 2008
coordinates ; the parent domain is recti-
linear even though the corresponding object-space subdomains
might be curvilinear. The word element indicates a subdomain
together with a set of expansion functions and associated de-
grees of freedom defined on it [27]. However, the term element
is often misused to indicate just the subdomain itself, and we
are not ourselves faultless on this score.
Several elements were defined and used in previous works.
Two-dimensional (2D) triangular and quadrilateral elements,
as well as three-dimensional (3D) tetrahedral and brick ele-
ments are, for example, discussed in [28], whereas prism and
pyramidal elements are given in [29] and [30], respectively. In
the following we assume the reader to be comfortable with the
definitions given in those papers, and adopt the same notation
used there to discuss application of the singularity cancellation
method to evaluate potential integrals.
All parametric coordinates are positive inside the element,
and every point outside of the element has one or more nega-
tive parent coordinates. For nested integration over the element
region
(1)
the
upper bounds must be obtained explicitly through the de-
pendency relationships (see [28]–[30]), thereby ensuring non-
negative coordinates within the integration region.
The number
of parametric coordinates used to describe a
given element is the size of the element, or equivalently, the size
of the set
. The size of a two-dimensional element is the number
of its edges whereas, for three-dimensional elements, the size
is the number of the element faces. The sizes of triangular and
quadrilateral elements are three and four, respectively; the sizes
of tetrahedrons, triangular-prisms, pyramids, and brick elements
are four, five, five, and six, respectively.
In the object space the element geometry is defined by
in-
terpolation (or control) points
, where
is a multi-index array of size , with integer entries
for . The element dependency relations further
constrain the integer entries
of the multi-index array, these
constraints depending on the shape, size and dimension of the
element [28]–[30]. The position vector in the object space is
then expressed in terms of
shape functions , usually
of polynomial form, attached to each interpolation or control
point
(2)
In the global object space, the element region is the entire set of
points
obtained by the mapping (2) of all the points of the
parent region
defined in (1).
III. S
INGULARITY CANCELLATION IN PARENT COORDINATES
Potential integrals on a given element are normally evaluated
by subdividing the object-space element region into subdomains
obtained by joining by a line each vertex of the entire domain
to the given observation point
or, for 2D elements, to its pro-
jection onto the element or its extension [5], [11]. For 2D el-
ements the subdomains are triangles whereas, in 3D, they are
tetrahedrons or pyramids (see also [1]). This subdivision nor-
mally requires the re-parameterization of each subdomain with
a sub-mapping of the kind given in (2).
Since the expansion functions are often given in terms of the
parent coordinates of the entire element, while the potential in-
tegral kernels are expressed in terms of the global coordinates
of the integration point, it is convenient to subdivide, whenever
possible, the entire element directly in its parent space so to
permit one to evaluate the kernel and the expansion functions
using only one mapping. Direct subdivision is possible when
the lines used to subdivide the object element are straight and
are mapped, in the parent element, again by straight lines joining
the parent vertices to a common point
. This happens, for ex-
ample, whenever the object element is defined by an affine trans-
formation of the parent element. In the following, the object el-
ements of this kind are called hyper-straight. For example, this
is the case for rectilinear object-elements defined as in (2) by
a mapping whose Jacobian is a non-zero constant (e.g., recti-
linear triangles and tetrahedra, parallelograms, hexahedra with
three pairs of parallel faces, etc.). In several applications, the
majority of the elements are of the kind just described, although
it is true that all the other cases have to deal with the more gen-
eral, but more computationally intensive subdivision technique
in the object space. Three fundamental properties are valid for
hyper-straight elements: 1) the integration coordinates defined
by our cancellation procedure define the parametric coordinates
of all sub-elements via the same formula (see (3), (4), (38),
(40) below); 2) the point about which the element is subdivided
is the center of similarity transformations (see Section V and
Appendix III); 3) the potential integrals on 3D hyper-straight
elements are simplified by the application of these transforma-
tions.
A general procedure to subdivide an element directly in its
parent space is presented in Appendix I. In the case of 2D el-
ements, this procedure yields to the results of Section III-A.
A discussion on application of this procedure to the case of
3D elements is deferred to Section V, with details reported in
Appendix III.
In general, for a given element of size
, a potential integral
on
over the element region is subdivided into sub-in-
tegrals, with integral subdomains obtained by joining a point
to each vertex of the parent domain (see
Fig. 1).
is the arbitrarily located common origin of different
local pseudo-radial frames introduced to locally perform each
sub-integral by properly changing the integration variables; the
Jacobian of each variable transformation vanishes at
.
In applications involving 3D elements,
is the parent point
that maps the observation point
of the global object-space; no-
tice that in the object space we use no superscript for the obser-
vation point
, whereas the source (i.e., integration) point is
primed. In applications involving planar 2D elements,
is the
parent point that maps, in the global space, the normal projec-
tion
of the observation point onto the element surface, or
onto its planar extension. If the 2D element is non-planar, the 2D
potential integrals are performed by working with a rectilinear
planar patch of the object space that is tangent to the original
curved one; in this case
is the point that maps the normal
projection
of the observation point onto this tangent patch.
GRAGLIA AND LOMBARDI: MACHINE PRECISION EVALUATION OF POTENTIAL INTEGRALS 983
Fig. 1. A domain of size
is broken into
subdomains by joining the point
=(
; ;
...
;
)
to each domain vertex. The figure illustrates the case
for elements of size four: (a) the two-dimensional quadrilateral element is sub-
divided about the point
into four triangular subdomains; (b) the three-di-
mensional tetrahedral element is broken about the point
into four tetrahedral
subdomains.
A. Transformation Formulas for 2D Elements
When applied to 2D elements, the transformation formulas
of Appendix I introduce a new pseudo-radial variable
, and
a pseudo-angular variable
. For 2D-elements, the pseudo-ra-
dial variable plays the same role as a radial variable of a cylin-
drical reference frame (centered at
); for this reason, we call
the pseudo-radial variable of 2D-elements, whereas we prefer to
call
the pseudo-radial variable of 3D elements. We now con-
sider the parent triangular domain
, as well as the quadrilateral
domain
shown in Fig. 1(a). As stated, these two domains are
split into triangular subdomains with a common vertex at
.
Each subdomain retains only one edge of the original element,
and
indicates the sub-triangle that retains the edge.
Our procedure maps each subdomain
into the square domain
.
The triangular domain
is split into three subdomains
( , 2, 3) with subscripts counted modulo three. For the
th subdomain, the variable transformation formulas (29) of
Appendix I read as follows
(3)
with a Jacobian
. The above result is certainly not new.
In fact, the current literature usually refers to the above transfor-
mation as the Duffy [4] or, sometimes, the Graglia transforma-
tion [5], [14], although the first user of this transformation is lost
in the mists of time. In fact, for example, the same transforma-
tion was used in 1971 by Tracey [2], and by Stern and Becker
in 1978 [3].
With reference to Appendix I, (3) is obtained by setting
as per (31); the dependency relation (see [28])
then yields . In this case, (3)
is obtained from (29) by dropping the subscript in
, that is
by setting
.
For the quadrilateral domain
, split into four subdomains
( , 2, 3, 4), one counts the subscripts modulo four to
write the two dependency relationships
and
(see [28]). For the th triangular subdomain
of
, the variable transformation formulas (29) of Appendix I
read as follows:
(4)
Fig. 2. A four-sided planar patch of the object space is broken into four tri-
angular subdomains about the normal projection
rr
r
of the observation point
rr
r
onto the plane of the patch: (a) the distance from
rr
r
to the plane of the patch is
d
; (b)
R
and
R
are the distances of the observation point from the two
vertices of the
k
th side, whose length is
`
.
with a Jacobian . The coordinate lines and
of the th subdomain of a quadrilateral patch are shown
in Fig. 2(b). With reference to Appendix I, (4) is obtained by
setting
as per (31); the dependency relationships then
yield
and . Equations (4) are then
obtained from (29) by simply dropping the subscript in
.
Notice that (3) and (4) fully comply with the dependency re-
lationships of the triangular and quadrilateral element, respec-
tively. Both (3) and (4) can deal with a point
located outside
the parent domain, or on its border. In fact, the integral on
and
is the algebraic sum of the sub-integrals on the subdomains
(for ) where, for located outside of the ele-
ment, at least one of the parametric coordinates of
(say )is
negative, so that the corresponding Jacobian
is also
negative whereas, for
lying on the border , the Jaco-
bian
vanishes together with the sub-integral on .
IV. P
OTENTIAL INTEGRALS ON 2D ELEMENTS
With reference to Fig. 2, let us consider a planar patch and a
nearby observation point
, where is the unit vector
normal to
, is the normal projection of onto the plane of
the patch, and
vanishes whenever lies on the patch-surface or
its extension. The singularity cancellation procedure is applied
to potential integrals of the form
(5)
where
is a vector or scalar basis function, is
the vector distance from integration to observation point, and
.
For each sub-triangle
, and in terms of the variables
defined in Section III-A, the distance from integration to obser-
vation point normalized with respect to the magnitude
of the
edge vector
reads as follows:
(6)
where
indicates the distance, normalized w.r.t. ,
of the observation point from the plane of the patch, whereas
and depend on the (unnormalized) distances , of
984 IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. 56, NO. 4, APRIL 2008
the observation point from the two vertices of the sub-triangle
opposite [see Fig. 2(b)]
(7)
with
(8)
and where
and
. Notice that in Fig. 2(b) increases clock-wise on
the sub-triangle
from to , because of (3) and
(4).
Successive variable transformations into parent (2) and then
into pseudo-radial (3), (4) coordinates yield
(9)
where
is the Jacobian of the transformation (2) between
global and parametric
-coordinates, and where it is henceforth
understood that the distance
from the observation to the inte-
gration point is expressed in terms of the integration variables
in use (see for example (6)). The pseudo-radial transformation
(3) or (4) produces a
factor in the integrand that cancels
out the singularity of (9) at
in the non-displaced case
of
. The modified Euler’s substitution (see [31])
reported in Table I, first column, reduces (9) to
(10)
where
is the new integration variable replacing , and with
(11)
As reported in Table I,
is a real function of independent of
, with and for all in the integration interval
[0,1]. The integral (10) is evaluated numerically by integrating
first along
, that is for , and then on . This integral
simplifies considerably when the observation points lies on the
patch-surface (self-element integration) or on its extension, that
is for
. In this case, the mapping reported at bottom of
Table I yields a constant value for
that does not depend
on
(12)
and the distance
simplifies into the rational function (see
Table I)
(13)
It is interesting to observe that the static form of (10) in the case
of a constant basis function
immediately yields
(14)
with
(15)
at
. Thus, for , the free-space static potential of a con-
stant source distribution is exactly integrated over a triangular
and quadrilateral element by using only three (
,2,3) and
four
sampling points, respectively. The result
(15) is equivalent to corresponding closed form results that can
be obtained from several others published papers [1], [10]–[16],
[19].
The pole in the kernel of (10) is easily cancelled by the tran-
scendental transformation reported in the right-hand column of
Table I, that yields
(16)
with
(17)
in the general case of
, and in
the case
. A variable transformation formula based on
the use of the hyperbolic sine function to regularize (5) directly
into an integral of the form of (16) has already been used in [1],
thereby obtaining a result very similar to (16). Transcendental
transformations of the same kind as [1] are also discussed in [8],
[9].
Our Euler modified transformation permits one to analyti-
cally explain the reasons why numerical integration of near-
singular integrals for small
values is much more difficult
than integration of the non-displaced singular integral, thereby
showing, very clearly, the different degree of difficulty one has
to deal with singular and near-singular displaced integrals. First
of all, the pole of the kernel of (10) varies with
in the near-sin-
gular displaced case, whereas
is -invariant in the singular
case. Secondly, after integration on , the remaining in-
tegral in
always exhibits a kernel with a large dynamic range
in the case of very small, non-vanishing
values. This latter re-
sult is apparent if one considers the static logarithmic kernel of
(14) and, in view of the integration by parts rule, it is intuitively
expected to hold also for potentials associated with non-con-
stant
functions. To further elaborate, the integral along in
