Engineering Analysis with Boundary Elements 27 (2003) 57–71

www.elsevier.com/locate/enganabound

A boundary cloud method with a cloud-by-cloud polynomial basis

Gang Lia, N.R. Alurub,*
a
Department of Mechanical and Industrial Engineering, Beckman Institute for Advanced Science and Technology,
University of Illinois at Urbana-Champaign, Urbana, IL 61801, USA
b
Department of General Engineering, Beckman Institute for Advanced Science and Technology, University of Illinois at Urbana-Champaign,
Urbana, IL 61801, USA
Received 18 April 2002; revised 6 September 2002; accepted 18 September 2002

Abstract
We have recently presented a boundary cloud method (BCM) [Comput. Meth. Appl. Mech. Engng 191 (2002) 2337], which combines
boundary integral formulations with scattered point interpolation techniques. A generalized least-squares approach, which requires
information about the outward normal to the boundary, is employed to construct interpolation functions. Since an outward normal is not well
defined for geometries with corners for 2D problems (or for corners and edges for 3D problems), points could not be placed at corners when
discretizing the surface of the object. In this paper, we introduce a new implementation of the BCM, which uses a varying base interpolating
polynomial to construct interpolation functions. The key idea is to define an appropriate polynomial basis which ensures linear completeness.
The polynomial basis can change from cloud to cloud depending on the definition of the cloud at each point. The new implementation can
handle points at corners and is much simpler and at least an order of magnitude faster compared to our original implementation. The original
implementation can be more accurate and can give higher order convergence rates, but is limited because it cannot handle points at corners.
Numerical results comparing the original and the new implementation are shown for several potential and electrostatic problems.
q 2002 Published by Elsevier Science Ltd.
Keywords: Boundary cloud method; Least-squares approach; Cloud by cloud; Polynomial basis; Meshless method

1. Introduction Boundary integral formulations when combined with

scattered point or meshless techniques (see Refs. [7 – 9] for
Boundary integral formulations [2], especially when overview on meshless methods) are attractive compu-
combined with fast algorithms based on multi-pole expan- tational methods for exterior problems as they alleviate
sions [3], fast Fourier transform (FFT) [4] and singular the mesh generation requirements. Recently, several mesh-
value decomposition (SVD) [5,6], are powerful compu- less methods for boundary-only analysis have been
tational techniques for rapid analysis of exterior problems. proposed in the literature. Some of the methods include
Boundary integral formulations have the advantage of the boundary node method [10 –12], the hybrid boundary
reducing the original dimensionality of the problem. For node method [13] and the boundary knot method [14]. The
example, boundary integral formulations reduce a 2D boundary node method is a combined boundary integral/
problem to a 1D problem and a 3D problem to a 2D meshless approach for boundary only analysis of partial
problem. The computational complexity of fast boundary differential equations. In a boundary node method, given a
element methods is typically OðN log NÞ; where N is the scattered set of points, meshless interpolation functions are
number of panels or elements used to discretize the surface constructed by using a moving least-squares approach [15].
of the object. A major burden in fast boundary element A straight-forward use of Cartesian coordinates in a moving
methods is the mesh generation task, which can be very least-squares approach can lead to singular coefficient
expensive and complicated for complex objects. matrices if the points on the boundary lie along a straight
line (see Ref. [1] for a more detailed discussion). The
* Corresponding author. boundary node method circumvents this problem by
E-mail address: [email protected] (N.R. Aluru). employing a cyclic coordinate (or curvilinear coordinates
0955-7997/03/$ - see front matter q 2002 Published by Elsevier Science Ltd.
PII: S 0 9 5 5 - 7 9 9 7 ( 0 2 ) 0 0 1 1 0 - 8
58 G. Li, N.R. Aluru / Engineering Analysis with Boundary Elements 27 (2003) 57–71

in 3D) in the moving least-squares approach to construct with the boundary integral formulations are summarized
interpolation functions. Recently, we have introduced a below.
boundary cloud method (BCM), which is also a combined
boundary-integral/scattered point approach for boundary 2.1. Potential problem
only analysis of partial differential equations. The BCM
employs a generalized least-squares approach to construct Consider an arbitrary domain V; as shown in Fig. 1, on
interpolation functions. The generalized least-squares which the potential equation along with the boundary
approach enables the use of Cartesian coordinates, which conditions is to be solved, i.e.
is attractive. However, since the approach requires infor-
72 u ¼ 0 in V ð1Þ
mation about the outward normal to the boundary, points
could not be placed at the corners (for 2D) or along the u¼g on Gg ð2Þ
edges and corners (for 3D) as the normal at these points is
not well-defined. ›u
q¼ ¼h on Gh ð3Þ
In this paper, we introduce a new implementation of ›n
BCM where points can be placed along corners and edges. where u is the unknown potential; Gg ; the portion of the
The basic idea is to use a varying polynomial basis to boundary where Dirichlet boundary conditions are speci-
eliminate the singularity encountered in the moving least- fied; Gh ; the portion of the boundary where Neumann
squares approach when all the points lie along a straight boundary conditions are specified; g and h, the prescribed
line. Specifically, when a weighting function is centered at a Dirichlet and Neumann boundary conditions, respectively;
point, a cloud for that point is defined. If all the points in the q, the normal derivative of u; n is the outward normal vector.
cloud lie along a straight line (or in the same plane for 3D), A boundary integral equation for the potential problem is
the moment matrix becomes singular when a complete given by (see, e.g. Ref. [2] for details)
linear base interpolating polynomial is used, i.e. {1; x; y}: ( )
ð ›ln rðP; QÞ
However, if the basis is chosen carefully, then the ln rðP; QÞqðQÞ 2 ½uðQÞ 2 uðPÞ dSQ ¼ 0
singularity can be avoided. For example, instead of using dV ›nQ
{1; x; y}; the use of either {1; x} or {1; y} depending on the ð4Þ
definition of the cloud can ensure the linear completeness as
well as eliminate the problems encountered in the moving where lnðrÞ is the Green’s function for the 2D potential
least-squares approach. Compared to the generalized least- equation, P and Q are the source and field points,
squares approach, the use of a varying polynomial basis is respectively, rðP; QÞ is the distance between P and Q and
attractive for two reasons: first, since the varying basis nQ is the outward normal vector at a field point Q on
approach does not require information about the outward boundary dV:
normal, points can be placed at corners (and edges in 3D);
second, the varying basis approach is much simpler and 2.2. Electrostatics
faster to use to compute interpolation functions. However,
as demonstrated in the results section, the generalized least- Consider a two conductor system as shown in Fig. 2. V1
squares approach is more accurate and produces higher and V2 denote the geometry of conductors 1 and 2,
convergence rates compared to the varying basis approach. respectively, dV1 and dV2 denote the surface or boundary
The rest of the paper is outlined as follows: Section 2 of conductors 1 and 2, respectively, and V
is the domain
summarizes the governing equations and the boundary exterior to V1 and V2 : A potential of g1 is applied on
integral formulations for potential and electrostatic pro- conductor 1 and a potential of g2 is applied on conductor 2.
The objective is to compute the surface charge density on
blems; Section 3 summarizes the generalized least-squares
approach and introduces the varying basis approach to
compute interpolation functions; Section 4 describes the
numerical implementation of the method; numerical results
for potential and electrostatic problems are shown in Section
5; conclusions are given in Section 6.

2. Governing equations and boundary-integral

formulations

We will focus on two-dimensional problems in this

paper, but the approach can be extended for three-
dimensional problems. The governing equations along Fig. 1. An arbitrary domain on which the potential equation is to be solved.
 G. Li, N.R. Aluru / Engineering Analysis with Boundary Elements 27 (2003) 57–71 59

the point t are approximated by

uðx; yÞ ¼ pT ðx; yÞat ð10Þ
T
›p ðx; yÞ
qðx; yÞ ¼ at ð11Þ
›n
where p is the base interpolating polynomial; at ; the
unknown coefficient vector for point t; n is the direction of
the outward normal to the boundary (note that n can be
different at every point). The linear basis interpolating
polynomial and its normal derivatives are given by
Fig. 2. A two conductor electrostatic system. pT ðx; yÞ ¼ ½ 1 x y m¼3 ð12Þ

the two conductors. The governing equation along with the  

›pT ðx; yÞ ›x ›y
boundary conditions for the exterior electrostatic problem is ¼ 0
given by [16] ›n ›n ›n

72 f ¼ 0 in V
ð5Þ ¼ ½0 b
cosðnx
nxÞ b 
cosðny
nyÞ m¼3 ð13Þ

f ¼ g1 on dV1 ð6Þ b and ðny

where ðnx
nxÞ b are the angles between the outward
nyÞ
f ¼ g2 on dV2 ð7Þ normal direction and the positive x- and y-axis, respectively.
For a point t, the unknown coefficient vector at is
A boundary integral equation for the electrostatic problem is computed by minimizing
given by [17]
X
NP X
NP
ð Jt ¼ wi ðxt ; yt Þ½pT ðxi ; yi Þat 2 u^ i 2 þ ~ i ðxt ; yt Þ
w
fðPÞ ¼ GðP; QÞsðQÞdSQ þ C ð8Þ i¼1 i¼1
dV
ð " #2
sðQÞdSQ ¼ CT ð9Þ ›pT ðxi ; yi Þ
 at 2 q^ i ð14Þ
dV ›n
where s is the unknown surface charge density; P, the where NP is the number of nodes, wi ðxt ; yt Þ is the weighting
source point; Q, the field point; G is the Green’s function function centered at ðxt ; yt Þ and evaluated at node i whose
and dV ¼ dV1 < dV2 : In 2D electrostatics, GðP; QÞ ¼ 2 coordinates are ðxi ; yi Þ: The weighting function for the
lnlP 2 Ql=ð2pe 0 Þ; where lP 2 Ql is the distance between second term in Eq. (14) is different from the weighting
the source point P and the field point Q, e 0 is the function for the first term to account for dimensionality
permittivity of the free space, CT is the total charge of the analysis. Typically, w ~ i ¼ wi for convenience. This means
system ðCT ¼ 0Þ and C is an unknown variable which needs that the scaling factor is set to one. u^ i and q^ i are nodal
to be computed. The unknown variable C arises because the parameters. The weighting function is nonzero when the
potential at infinity is not zero for two-dimensional location of node i is within a certain distance from point
problems. ðxt ; yt Þ: Thus, for a point ðxt ; yt Þ; the weighting function is
nonzero only for a few other nodes in the vicinity of it. The
region where the weighting function is nonzero is called a
3. Scattered point interpolation cloud. The weighting function is typically a cubic spline or a
Gaussian function. In this paper, a cubic spline is used
The first step in the BCM is to construct interpolation whose definition is given by
functions over a scattered set of points. As shown in Fig. 1,

the surface of the object is represented by scattered points. 1 x 2 xi
wi ðxt Þ ¼ ^ t
w ð15Þ
The objective is to approximate unknown quantities by dx dx
interpolation functions. The Hermite-type approach that 8
> 0 z , 22
was introduced in Ref. [1] is first quickly summarized >
>
>
>
> ðzÞ =6
3
before the varying basis approach is presented to construct > 22 # z # 21
>
>
interpolation functions. >
< 2=3 2 z ð1 þ z=2Þ 21 # z # 0
2
^ ¼
wðzÞ ð16Þ
>
> 2=3 2 z2 ð1 2 z=2Þ 0 # z # 1
3.1. Hermite-type interpolation >
>
>
>
>
> 2ðzÞ3 =6 1#z#2
In a Hermite-type interpolation, given a point t, >
>
:
the unknown and its normal derivative in the vicinity of 0 z.2
60 G. Li, N.R. Aluru / Engineering Analysis with Boundary Elements 27 (2003) 57–71

where z ¼ ðxt 2 xi Þ=dx and dx denotes the support size of the

weighting function in the x-direction. A multi-dimensional
weighting function can be constructed as a product of one-
dimensional weighting functions. In two dimensions, the
weighting function is given by

 !
1 xt 2 xi 1 yt 2 yi
wi ðxt ; yt Þ ¼ ^
w w^ ð17Þ
dx dx dy dy

where dy is the support size in the y direction.

By minimizing Jt the necessary interpolation functions Fig. 3. Definition of singular and nonsingular cloud.
can be computed. Once the interpolation functions are
computed, the unknowns are evaluated by

X
NP X
NP
uðx; yÞ ¼ ~ I ðx; yÞ^uI þ
M N~ I ðx; yÞ^qI ð18Þ
shown in Fig. 3, the cloud of point 1 is singular and the
I¼1 I¼1
cloud for point 2 is nonsingular as the points do not lie along
a straight line. In this paper, we use a linear polynomial
X
NP X
NP basis. The base interpolating polynomial is given by
qðx; yÞ ¼ S~ I ðx; yÞ^uI þ T~ I ðx; yÞ^qI ð19Þ
I¼1 I¼1
pTv ðx; yÞ
where M ~ I ðx; yÞ; N~ I ðx; yÞ; S~ I ðx; yÞ and T~ I ðx; yÞ are the (
½1 x y m ¼ 3 ðbasis for nonsingular cloudÞ
interpolation functions and are given by ¼
½ 1 x  or ½ 1 y  m ¼ 2 ðbasis for singular cloudÞ
t pðxI ; yI ÞwI ðxt ; yt Þ
~ I ðx; yÞ ¼ pT ðx; yÞC21
M ð20Þ ð25Þ

0 For a point t, the unknown coefficient vector bt is computed

t p ðxI ; yI ÞwI ðxt ; yt Þ
N~ I ðx; yÞ ¼ pT ðx; yÞC21 ð21Þ
by minimizing

›pT ðx; yÞ 21 X
NP
S~ I ðx; yÞ ¼ Ct pðxI ; yI ÞwI ðxt ; yt Þ ð22Þ Jt ¼ wi ðxt ; yt Þ½pTv ðxi ; yi Þbt 2 u^ i 2 ð26Þ
›n
i¼1

›pT ðx; yÞ 21 0 where NP is the number of nodes, wi ðxt ; yt Þ is the weighting

T~ I ðx; yÞ ¼ Ct p ðxI ; yI ÞwI ðxt ; yt Þ ð23Þ function centered at ðxt ; yt Þ and evaluated at node i whose
›n
coordinates are ðxi ; yi Þ and u^ i is a nodal parameter. The
~ I ðx; yÞ; N~ I ðx; yÞ; S~ I ðx; yÞ and T~ I ðx; yÞ are multi-
Note that M weighting function is again the cubic spline defined in Eqs.
valued. We limit the interpolation functions to a single value (15) –(17).
by computing the interpolations at the point where the The stationary of Jt leads to
weighting function is centered (see Refs. [1,19,20] for
details). ðPTv WPv Þbt ¼ PTv Wu^ ð27Þ

Eq. (27) can be rewritten as

3.2. Varying basis least-squares approximation
C0t bt ¼ A0t u^ ð28Þ
In this section, we introduce a varying basis least-squares
approach to approximate the unknown quantity, i.e.
bt ¼ C0t 21 A0t u^ ð29Þ
uðx; yÞ ¼ pTv ðx; yÞbt ð24Þ
where C0t is an m £ m matrix and A0t is an m £ NP matrix,
where pv is the varying base interpolating polynomial and bt whose definitions are given by
is the unknown coefficient vector for point t. To construct
the varying basis interpolation functions, clouds are C0t ¼ PTv WPv ð30Þ
classified into two types: singular and nonsingular. When
all the points inside a cloud lie along a straight line, the
cloud is defined as singular, otherwise, it is nonsingular. As A0t ¼ PTv W ð31Þ
 G. Li, N.R. Aluru / Engineering Analysis with Boundary Elements 27 (2003) 57–71 61

where Remarks.
2 3 2 3
pTv ðx1 ; y1 Þ u^ 1 1. A fixed least-squares approach [18] is used to illustrate the
6 7 6 7
6 pT ðx ; y Þ 7 6 u^ 7 construction of interpolation functions by a varying basis
6 v 2 2 7 6 2 7
6 7 6 7
Pv ¼ 6 7; u^ ¼ 6 . 7; approach. A moving least-squares approach can also be
6 .. 7 6 . 7
6 . 7 6 . 7 used to define interpolation functions. The interpolation
4 5 4 5
functions constructed by either the fixed least-squares
pv ðxNP ; yNP Þ
T
u^ NP approach or by the moving least-squares approach are
2 3 ð32Þ
identical. However, the derivatives of the interpolation
w1 ðxt ; yt Þ 0 0 0
6 7 function can be different. Since the derivatives of the
6 w2 ðxt ; yt Þ 0 7
6 0 0 7 interpolation function are not required for the work
6 7
W¼6 7 presented in the paper, both moving least-squares and
6 . . 7
6 0 0 . 0 7 fixed least-squares approaches produce identical results.
4 5
2. Both the Hermite-type approach and the varying basis
0 0 0 wNP ðxt ; yt Þ
approach avoid singular moment matrices when all the
points in the cloud lie along a straight line. Unlike the
As shown in Ref. [1], when a cloud is singular, the moment
Hermite-type approach, the varying basis approach over-
matrix C0t becomes singular if the full polynomial basis
comes the difficulty by appropriately redefining the basis
½1 x y is used. In such cases, the appropriate base
for the singular clouds.
interpolating polynomial can be determined from the
3. Both the Hermite-type and the varying basis interpola-
information of points contained in the cloud. For example,
tions functions are meshless interpolation functions. No
if all the points in the singular cloud are defined by a line
connectivity information among the scattered points is
y ¼ constant; then pTv ¼ ½1 x can be used. If the singular
required to compute the interpolation functions.
cloud is defined by a line x ¼ constant; then pTv ¼ ½1 y can
4. Both the Hermite-type and the varying basis interpolants
be used. In the general case of a cloud defined by y ¼
are consistent, i.e. any linear combination of the basis
ax þ b; either pTv ¼ ½1 x or pTv ¼ ½1 y can be used.
polynomial can be interpolated exactly. The varying basis
Substituting the definition of bt from Eq. (29) into Eq.
interpolants are consistent even when the basis is reduced.
(24), the approximation for the unknown can be rewritten as
For example, assigning to each u^ i the values of basis
uðx; yÞ ¼ pTv ðx; yÞbt ¼ pTv ðx; yÞC0t 21 A0t u^ ð33Þ functions gives
u^ ¼ Pv ð40Þ
In short form,
Substituting Eq. (40) to Eq. (33)
uðx; yÞ ¼ Nðx; yÞu^ ð34Þ
uðx; yÞ ¼ pTv ðx; yÞC0t 21 A0t u^ ¼ pTv ðx; yÞC0t 21 PTv WPv
where Nðx; yÞ is 1 £ NP interpolation vector given by
¼ pTv ðx; yÞC0t 21 C0t ¼ pTv ðx; yÞ ð41Þ
Nðx; yÞ ¼ pTv ðx; yÞC0t 21 A0t ð35Þ
Without losing generality, assume a singular cloud is
Nðx; yÞ is again multi-valued and the interpolation function represented by a reduced basis pTv ¼ ½1 x when all the
can be limited to a single value by computing the points in the cloud lie along a straight line defined by
interpolations at the point where the weighting function is y ¼ kx: Assuming the unknown u takes the form u ¼
centered. Eq. (35) can be rewritten in a matrix form as, a þ bx þ cy (a, b and c are constants), for the points in the
u ¼ Nu^ ð36Þ singular cloud, u is given by u ¼ a þ ðb þ kcÞx: Since any
linear combination of pTv ðx; yÞ ¼ ½1 x can be interpolated
where N is an NP £ NP interpolation matrix. At any point exactly, the unknown u is interpolated exactly by the
ðx; yÞ; the unknowns are evaluated by reduced basis ½1 x:
5. The varying basis least-squares technique is a much
X
NP
uðx; yÞ ¼ NI ðx; yÞ^uI ð37Þ simpler approach compared to the Hermite-type
I¼1 approach. The varying basis approach is a much faster
technique to compute interpolation functions compared to
where
the Hermite-type approach. For a typical potential
NI ðx; yÞ ¼ pTv ðx; yÞC0t 21 pv ðxI ; yI ÞwI ðxt ; yt Þ ð38Þ problem with 512 nodes, the solution of the new
implementation is about 10 times faster than the original
The normal derivative of u can be interpolated by using the Hermite-type approach.
same interpolation function, i.e. 6. In a Hermite-type approach, the outward normal to the
X
NP boundary is required to compute interpolation func-
qðx; yÞ ¼ NI ðx; yÞ^qI ð39Þ tions. Hence, points could not be placed at corners (or
I¼1 edges in 3D) where the normal to the boundary is not
62 G. Li, N.R. Aluru / Engineering Analysis with Boundary Elements 27 (2003) 57–71

well-defined. A varying basis approach eliminates this NC ð

X
problem as it does not require the outward normal and Gði; jÞ ¼ ln rðPi ; Qk ÞNj ðQk ÞdSQ ð47Þ
k¼1 dSk
points can be placed at corners (and edges in 3D).
7. In a Hermite-type approach, the normal derivative of In a general mixed boundary value problem, either u or q is
the unknown (i.e. q ) is also interpolated exactly when specified at each point. By substituting all the known
the unknown function (i.e. u ) is interpolated exactly. quantities given by the boundary conditions into a vector x;
For example, if the base interpolating polynomial is a 2NP £ 2NP linear system is obtained
linear, and if u is linear then a constant q is " #" # " #
interpolated exactly. If q is constant, but is discontinu- F 2G u^ 0
0 00
¼ ð48Þ
ous from one boundary edge to the other boundary N N q^ x
edge, it can still be interpolated exactly in a Hermite-
type approach. A varying basis approach does not where N0 and N00 are obtained by appropriately selecting the
possess this property and typically error is introduced rows of N corresponding to the known us and qs. u^ and q^ are
in approximating a constant q. An example demon- computed by solving the above linear system. Once u^ and q^
strating this aspect is discussed in Section 5. are known, u and q can be computed from Eqs. (44) and (45).

4.2. Electrostatic problem

4. Numerical implementation
The discretized boundary integral equation for exterior
As described in Ref. [1], the boundary of the domain is electrostatics is given by
discretized into cells for integration purpose. Each cell NC ð
X
contains a certain number of nodes and the number of nodes fðPÞ ¼ GðP; Qk ÞsðQk ÞdSQ þ C ð49Þ
k¼1 dSk
can vary from cell to cell. The concept of a cell is different
from that of element or panel in boundary-element methods. NC ð
X
The cell can be of any shape or size and the only restriction sðQk ÞdSQ ¼ CT ð50Þ
k¼1 dSk
is that the union of all the cells equal to the boundary of the
domain. The varying basis interpolation of the surface charge density
is given by
4.1. Potential problem
X
NP
sðx; yÞ ¼ NI ðx; yÞs^I ð51Þ
The discretized boundary integral equation for the I¼1
potential problem is given by
Substituting Eq. (51) into Eqs. (49) and (50) for each source
NC ð
X NC ð
X ›ln rðP; Qk Þ node, the final matrix form is given by
lnrðP; Qk ÞqðQk ÞdSQ 2
k¼1 dSk k¼1 dSk ›nQ k Ax ¼ b ð52Þ
NC ð
X where
›ln rðP; Qk Þ
 uðQk ÞdSQ þ uðPÞ dSQ ¼ 0 ð42Þ 8
dSk ›nQ k > XNC ð
k¼1 >
>
>
> Aði; jÞ ¼ GðPi ; Qk ÞNj ðQk ÞdSQ i ¼ 1; …; NP
where NC is the total number of cells and Qk is the field >
> k¼1 dSk
>
>
point in cell k. >
< NC ð
X
Substituting Eqs. (37) and (39) into Eq. (42) for each
> AðNP þ 1; jÞ ¼ Nj ðQk ÞdSQ j ¼ 1; …NP ð53Þ
node, the final matrix form of the BCM with a varying basis >
> k¼1 dSk
>
>
is given by >
> Aði; NP þ 1Þ ¼ 1 i ¼ 1; …; NP
>
>
>
:
Fu^ ¼ Gq^ ð43Þ AðNP þ 1; NP þ 1Þ ¼ 0
u ¼ Nu^ ð44Þ 8 9 8 9
> s1 > > f1 >
>
> >
> >
> >
>
>
> >
> >
> >
>
q ¼ Nq^ ð45Þ >
> s 2 >
> >
> f 2 >
>
>
< >
= >
< >
=
where F and G are NP £ NP matrices given by .. ..
x¼ . > ; b ¼ . > ð54Þ
>
> > >
> >
NC ð
X >
> >
> >
> >
>
›ln rðPi ; Qk Þ >
> s >
> >
> s NP >
>
Fði; jÞ ¼ Nj ðQk ÞdSQ 2 Nj ðPi Þ >
>
NP >
> >
> >
>
k¼1 dSk ›nQ k : ; : ;
C CT
NC ð
X ›ln rðPi ; Qk Þ
dSQ ð46Þ The integrals in Eqs. (46), (47), and (53) are weakly singular
k¼1 dSk ›nQk and near singular when the source point P and the field point
 G. Li, N.R. Aluru / Engineering Analysis with Boundary Elements 27 (2003) 57–71 63

Q coincide or very close to each other, respectively. Regular The results are compared with the results of BCM with
Gauss quadrature becomes inaccurate in these cases. Special Hermite-type interpolations. Unless otherwise mentioned
treatment is employed to compute the various integrals. we use the 2D cubic spline weighting function as given in
Given a source point P, cells are classified into three Eq. (17). The convergence of the method is measured by
categories: (1) for cells which are far away from the source using a global error measure
point, a regular Gauss quadrature is used to compute the vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
u
integral, (2) for cells which are nearby to the source point, 1 u 1 X NP
e ¼ ðeÞ t ½uðeÞ 2 uIðcÞ 2 ð55Þ
the weights of the nodes are computed more accurately by lu lmax NP I¼1 I
using the Nyström [21] and the SVD [22] techniques, and
(3) for a cell which includes the source point, the integral is where e is the error in the solution and the superscripts (e)
computed analytically. Extensive details on computing and (c) denote, respectively, the exact and the computed
these integral can be found in Ref. [1]. solutions.

5. Numerical examples 5.1. 2D Potential examples

In this section we present results for several problems Consider the solution of the potential equation on various
using the BCM with varying base interpolating polynomial. domains as shown in Fig. 4. Fig. 4(a) is a 1 £ 1 square

Fig. 4. (a) 1 £ 1 square domain on which the potential equation is solved, (b) a circular domain, and (c) an elliptical domain.
64 G. Li, N.R. Aluru / Engineering Analysis with Boundary Elements 27 (2003) 57–71

domain where Dirichlet or Neumann boundary conditions the boundary of the domain is described by
are specified along the four edges of the boundary. Four
examples with different boundary conditions are presented. x ¼ a·cosðuÞ; y ¼ b·sinðuÞ ð56Þ
For a uniform distribution of nodes with a linear basis, the
where a ¼ 1:0 and b ¼ 0:5: A Neumann problem is solved
cloud sizes are chosen to be rx ¼ ry ¼ 1:17Ds; where Ds is
over the elliptical domain. The cell structure and the cloud
the spacing between the points. A uniform cell structure is size are identical to those used for the circular domain.
used. Fig. 4(b) is a circular domain with a radius of 1.0. A The first example is a Dirichlet problem on the square
Dirichlet problem is solved over the circular domain. Nodes domain with a linear exact solution for u and a constant
are uniformly distributed along the boundary of the domain. exact solution for q. Dirichlet boundary conditions are
A uniform cell structure with two Gauss points per cell is specified along the four edges of the boundary. The
considered. The cloud sizes are chosen to be rx ¼ ry ¼ 1:17 boundary conditions are
Ds where Ds is the perimeter of the circle divided by the
number of points. Fig. 4(c) is an elliptical domain and u¼xþy for x ¼ 0; y ¼ 0; x ¼ 1; y ¼ 1 ð57Þ

Fig. 5. (a) Comparison of the computed and exact solution, (b) convergence of q with varying basis BCM. Note that no points are put at the corners. The
boundary is defined between 0 and 1 and no points are put at either 0 or 1. Even though points can be put at the corners with the new implementation, no points
were put at the corners in order to compare the solution with the original implementation.
 G. Li, N.R. Aluru / Engineering Analysis with Boundary Elements 27 (2003) 57–71 65

Fig. 6. Comparison of convergence of BCM with varying basis interpolation and Hermite-type interpolation for a Dirichlet problem with quadratic exact
solution for u.

Fig. 7. Comparison of convergence of BCM with varying basis interpolation and Hermite-type interpolation for a Neumann problem with cubic exact solution
for u.
66 G. Li, N.R. Aluru / Engineering Analysis with Boundary Elements 27 (2003) 57–71

Fig. 8. Convergence of BCM with varying basis interpolation and a random point distribution for a Dirichlet problem with quadratic exact solution for u.

Fig. 9. Comparison of results of BCM with varying basis interpolation and exact solution for a Dirichlet problem with quadratic exact solution for u.
 G. Li, N.R. Aluru / Engineering Analysis with Boundary Elements 27 (2003) 57–71 67

Fig. 10. Convergence of BCM for a Dirichlet problem on a circular domain.

The exact solution for q is given by Note that the normal derivative of u is discontinuous at
( points ð0; 1Þ and ð1; 0Þ: At these two points, q jumps from
›u 21 for y ¼ 0; x ¼ 0
q¼ ¼ ð58Þ 2 1 to 1. For this example, Gaussian quadrature is used for
›n 1 for x ¼ 1; y ¼ 1 all the cells. However, the weights of the Gauss points are

Fig. 11. Convergence of BCM for a Neumann problem on an elliptical domain.

68 G. Li, N.R. Aluru / Engineering Analysis with Boundary Elements 27 (2003) 57–71

Fig. 12. A two plate conductor system. Also shown is the discretization of the geometry into points.

computed such that the integration of lnðP; QÞ; x·lnðP; QÞ interpolation functions give a better performance by
and y·lnðP; QÞ are exact. Thus, in this case, the numerical including the normal derivative of the basis into the
error comes only from the interpolation. We obtain the exact interpolation. The convergence rate of the BCM with
solution by using the Hermite-type interpolations. However, Hermite-type interpolation is 1.41.
interpolation error arises near the two corners where q is The third example is a Neumann problem on the square
discontinuous when the varying basis least-squares interp- domain with a cubic exact solution for u. The specified
olation is used. Comparison of the numerical results and the boundary conditions are
exact solution is shown in Fig. 5(a). When the point 8
distribution is refined (8, 16, 32, 64 points per edge and 2 >
> 23y2 for x ¼ 0
>
>
points per cell), the numerical error is reduced. As shown in >
< 23x2
›u for y ¼ 0
Fig. 5(b), the convergence rate is 0.5. This example clearly ¼ ð61Þ
›n >
> 2
3y þ 6y 2 3 for x ¼ 1
shows one advantage of the Hermite-type interpolation, i.e. >
>
>
: 2
the capability to capture the discontinuity in the normal 3x þ 6x 2 3 for y ¼ 1
derivative of the solution around the corner.
The second example is a Dirichlet problem on the square u ¼ 2x3 2 y3 þ 3x2 y þ 3xy2 for x ¼ 0:97; y ¼ 0 ð62Þ
domain with a quadratic solution for u. Dirichlet boundary
conditions are specified along the four edges of the The exact solution for u is given by
boundary. The boundary conditions are u ¼ 2x3 2 y3 þ 3x2 y þ 3xy2
ð63Þ
u ¼ x2 2 y2 for x ¼ 0; y ¼ 0; x ¼ 1; y ¼ 1 ð59Þ
for x ¼ 0; y ¼ 0; x ¼ 1; y ¼ 1
The exact solution for q is given by The convergence rate of the BCM with varying basis and the
8 BCM with Hermite-type interpolations is shown in Fig. 7
>
> 0 for y ¼ 0; x ¼ 0
›u < employing a uniform distribution of 8, 16, 32 and 64 points
q¼ ¼ 2 for x ¼ 1 ð60Þ per edge and 4 points per cell. The accurate integration
›n >
>
: distance is set to be 0.3. The convergence rate of the BCM
22 for y ¼ 1
with varying basis is found to be 2.09. The convergence rate
The exact solution of u is quadratic which is outside the is found to be 2.06 with the Hermite-type interpolation.
linear base interpolating polynomial. This problem is Note that for this example the potential, u, and its normal
solved by the varying basis BCM with a uniform derivative, q, are continuous along the boundary.
distribution of 8, 16, 32 and 64 points per edge and 2 The fourth example is identical to the second example
points per cell to study the convergence behavior. The except that the nodes are distributed randomly on the
accurate integration distance, defined in Ref. [1], is set to boundary. An accurate integration distance of 0.3 is used for
be 0.5. A comparison of convergence of the BCM with this example. The convergence plot with a random point
varying basis and the BCM with Hermite-type interp- distribution of 8, 16, 32 and 64 nodes per edge is shown in
olation is shown in Fig. 6. The normal derivative of the Fig. 8. The exact solution and the numerical solution
potential is not continuous at the points ð1; 0Þ; ð0; 1Þ and obtained with a random distribution of 64 nodes per edge are
ð1; 1Þ: This discontinuity is not captured by the varying compared in Fig. 9. The cloud size in this case is rx ¼ ry ¼
basis BCM. The convergence rate of the varying basis 1:67Ds; where Ds is the average spacing between the points.
BCM is found to be 0.43. In comparison, the Hermite The convergence rate is found to be 1.35.
 G. Li, N.R. Aluru / Engineering Analysis with Boundary Elements 27 (2003) 57–71 69

Fig. 13. Surface charge density distribution on the conductors obtained by the varying basis BCM.

The next example is a Dirichlet potential problem on the The accurate integration distance for the circular domain is
circular domain shown in Fig. 4(b). For the Dirichlet set to be 0.8. The convergence of the BCM is shown in Fig. 10.
problem, the boundary conditions are The convergence rate is found to be 3.61 for the Hermite-type
interpolation and 1.85 for varying basis interpolation. The
u ¼ r·cosðuÞ ¼ x on all the boundary ð64Þ
varying basis interpolation shows a lower convergence rate
The exact solution for q is compared to the regular Hermite-type interpolation.
The next example is a Neumann potential problem on the
›u
¼ cosðuÞ on all the boundary ð65Þ elliptical domain shown in Fig. 4(c). The boundary
›n

Fig. 14. Convergence of the surface charge density at point pðsp Þ:

70 G. Li, N.R. Aluru / Engineering Analysis with Boundary Elements 27 (2003) 57–71

Fig. 15. A two conductor system with a complex geometry for the ground plane. Also shown is the point distribution.

Fig. 16. Surface charge density distribution on the conductors obtained by the varying basis BCM.

conditions are The exact solution for u is

›u b·cosðuÞ u¼x on all the boundary ð68Þ
b ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
¼ cosðnx
nxÞ
›n b cos ðuÞ þ a2 sin2 ðuÞ
2 2
ð66Þ The accurate integration distance for the elliptical domain is
set to be 0.8. The convergence of the BCM is found to be
on all the boundary
2.61 (Fig. 11) for the Hermite-type BCM and 2.58 for the
u¼x at two points on the boundary ð67Þ varying basis BCM.
 G. Li, N.R. Aluru / Engineering Analysis with Boundary Elements 27 (2003) 57–71 71

5.2. 2D Exterior electrostatic analysis References

The first example is a two plate conductor system as [1] Li G, Aluru NR. Boundary cloud method: a combined scattered point/
shown in Fig. 12. Each conductor has dimensions of boundary integral approach for boundary-only analysis. Comput Meth
Appl Mech Engng 2002;191(21/22):2337–70.
100 £ 2 and the gap between the two plates is 1. The [2] Kane JH. Boundary element analysis in engineering continuum
applied voltage is set to be 1 for the upper plate and 0 for the mechanics. Englewood Cliffs, NJ: Prentice-Hall; 1994.
lower plate. The permittivity of the free space e 0 is also set [3] Greengard L, Rokhlin V. A fast algorithm for particle simulations.
to be 1 for simplicity. Since the plate is very long compared J Comput Phys 1987;73(2):325–48.
to the gap the surface charge density s at the center point [4] Phillips JR, White JK. A precorrected-FFT method for electrostatic
analysis of complicated 3-D structures. IEEE Trans Comput-Aid Des
p is given by Integr Circuits Syst 1997;16(10):1059–72.
[5] Kapur S, Long DE. IES3: a fast integral equation solver for efficient
e 0 Vab 3-dimensional extraction. IEEE Computer Aided Design, 1997,
s¼ ð69Þ
d Digest of Technical Papers, IEE/ACM International Conference;
1997. p. 448–55.
where Vab is the potential difference between the two plates [6] Shrivastava V, Aluru NR. A fast boundary cloud method for exterior
2-D electrostatics. Int J Numer Meth Engng 2003.
and d is the gap. The surface charge density distribution is
[7] Belytschko T, Krongauz Y, Organ D, Fleming M, Krysl P. Meshless
shown in Fig. 13 and the convergence of the surface charge methods: an overview and recent developments. Comput Meth Appl
density at point p is shown in Fig. 14. Mech Engng 1996;139:3–47.
The second example shows that the BCM with the [8] Li S, Liu WK. Meshfree and particle methods and their applications.
varying basis interpolation is attractive for solving compli- Appl Mech Rev 2002;55:1–34.
cated MEMS problems. Shown in Fig. 15 is a circular [9] Atluri SN. The meshless local Petrov-Galerkin (MLPG) method. Tech
Science Press; 2002.
conductor located above a wave shaped ground plane. The [10] Mukherjee YX, Mukherjee S. The boundary node method for
surface charge density distribution obtained by the varying potential problems. Int J Numer Meth Engng 1997;40:797–815.
basis BCM is shown in Fig. 16. The computed surface [11] Chati MK, Mukherjee S. The boundary node method for three-
charge densities at points p and e are 0.18 and 0.198, dimensional problems in potential theory. Int J Numer Meth Engng
respectively. 2000;47:1523 –47.
[12] Chati MK, Mukherjee S. The boundary node method for three-
dimensional linear elasticity. Int J Numer Meth Engng 1999;46:
1163– 84.
[13] Zhang J, Yao Z, Li H. A hybrid boundary node method. Int J Numer
6. Conclusion Meth Engng 2002;53(4):751 –63.
[14] Chen W. Symmetric boundary knot method. Engng Anal Bound Elem
A new implementation of the BCM is presented in this 2002;26(6):489–94.
paper. The new implementation uses a varying polynomial [15] Lancaster P, Salkauskas K. Surface generated by moving least squares
basis when constructing the interpolation functions. methods. Math Comput 1981;37:141– 58.
[16] Jackson JD. Classical electrodynamics, 3rd ed. New York: Wiley;
Depending on whether the cloud is singular or non- 1999.
singular, the base interpolating polynomial is defined to [17] Shi F, Ramesh P, Mukherjee S. On the application of 2D potential
satisfy the linear completeness condition. The new theory to electrostatic simulation. Commun Numer Meth Engng 1995;
implementation can handle points at corners and edges 11:691–701.
and is much simpler to implement and faster compared to [18] Onate E, Idelsohn S, Zienkiewicz OC, Taylor RL, Sacco C. A
stabilized finite point method for analysis of fluid mechanics
our original implementation. Numerical results indicate problems. Comput Meth Appl Mech Engng 1996;139:315– 46.
that the Hermite-type BCM is more accurate and typically [19] Aluru NR, Li G. Finite cloud method: a true meshless technique based
produces higher order convergence rate compared to the on a fixed reproducing kernel approximation. Int J Numer Meth
varying basis BCM. Engng 2001;50(10):2373–410.
[20] Jin X-Z, Li G, Aluru NR. On the equivalence between least-squares
and kernel approximation in meshless methods. Comput Model Engng
Sci 2001;2(4):447–62.
Acknowledgements [21] Kapur S, Long DE. High-order Nyström schemes for efficient 3-D
capacitance extraction. 38th International Conference on Computer-
Aided Design, San Jose, CA, USA; 1998. p. 178– 85.
This work was supported by the National Science [22] Golub GH, Van Loan CF. Matrix computations. Baltimore MD: Johns
Foundation under Grants CCR-9875671 and CCR-0121616. Hopkins University Press; 1989.