# The FMM-BEM Method for the 3D Particulate Stokes Flow

Authors: Hassib Selmi; Mohamed Abdelwahed; Hatem Hamda; Lassaad El Asmi
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

This work introduces new functions based on the spherical harmonics and the solid harmonics which have been used to construct a multipole development for the 3D Stokes problem in order to reduce the operations costs in the BEM method. We show that the major properties of those functions are inherited from the solid harmonics. The contribution of this paper is the introduction of new formulas that serve to calculate the multipole moments and the transfer functions that are necessary for the schemes of orderO
(
N
log
N
). Moreover, new translation formulas are introduced to obtain an O
(
N
) scheme. The error truncation of the resulting scheme is discussed. In comparison to the BEM that attains a limit storage at O
(
1
0
4
), we present here a method based on FMM-BEM that attains a storage at a limit of O
(
1
0
6
). The implementation of the method achieves a high accuracy level at a reasonable cost.


---

## Body

## 1. Introduction

The studies of complex fluids, such as particulate suspensions, emulsions, and sedimentation problem, remain a great challenge in spite of their omnipresence in many physical, chemical, or biomedical processes and industrial applications. This problem represents physically the creeping flow around several solid particles moving in a viscous fluid. Due to the size of the particles, this kind of problems is modelled by the linear Stokes equation. This linearity allows the application of the surface integral equation method or the boundary element method (BEM) which may be more efficient than the boundary value methods such as the finite element method (FEM) or the finite difference method (FDM). Indeed, the advantage of the BEM is that only the surface is discretized, generating fewer elements than a volume discretization and making the system matrix much smaller. However, the integral equation method is a global approach; every point is affected by the points of the entire system, giving a dense matrix. Since there are efficient algorithms, that is, GMRES, for solving the sparse matrix generated by FEM and FDM, it was believed that BEM was more expensive to use than FEM or FDM. The development of the fast multipole method (FMM) with the aim of accelerating the product matrix vector and reducing the complexity of such operations fromO
(
N
2
) to O
(
N
log
⁡
N
) or even to O
(
N
), and the memory requirements to O
(
N
log
⁡
N
) or O
(
N
) [1–3], brought new life to the classical BEM. Indeed such method approximates the effect of the far field points of the entire system. Like in FEM, only the terms in a certain band of the matrix are stored in an exact way and all the others are approximated and compressed.We introduce in this paper new multipole moments and transfer functions forO
(
N
log
⁡
N
) scheme and new formulas for the translations from multipole to multipole, multipole to local, and local to local that are necessary to construct the O
(
N
) scheme for solving the 3D particulate Stokes problem. The multipole and local expansions can be useful to develop an efficient method based on the FMM to determine a real time solution of the dynamics of rigid bodies in a 3D infinite or semi-infinite viscous fluid. This solution is obtained by determining the stresses exerted by the fluid on each body and then the velocities of the bodies. These stresses are represented by Stokeslets distributed over the surfaces of the bodies. We begin in the next section by presenting a brief overview of the literature. In Section 3, we present the principle of the fast multipole method. Section 4 presents the two representations used in FMM formulas for the Laplace problem: the first is based on spherical harmonic [2] and the second on the solid harmonic [4]. We present also in this section the new functions that will be used to represent differently the Green function of the Laplace equation. For the sake of conciseness, clarity, and complexity of the algorithm, we choose to work with this new representation with the aim of developing the FMM formulas that permit solving the Stokes problem in an O
(
N
log
⁡
N
) scheme in the first step and in O
(
N
) scheme in the final step. In Section 5, we present a detailed description of the FMM formulas founded to be applied to the Green function of the Stokes equation. We also present a discrete algorithm as well as the study of the error bound. Finally, in Section 6, we present, discuss, and interpret the numerical results obtained from a cluster.


## 2. Overview

FMM reduces the cost of the product matrix vector fromO
(
N
2
) to O
(
N
log
⁡
N
) or even O
(
N
) by approximating the calculation by handling numerical series. It uses a hierarchical subdivision of space into panels or clusters of sources calculating the multipolar moment that are used in the evaluation of the far field expansion, and it reduces the space of storage from O
(
N
2
) for the dense matrix of size N to O
(
N
log
⁡
N
) or even to O
(
N
). From the algorithmic point of view, the FMM works down and up the TREE constructed by the recursive subdivision of the initial domain. FMM was initially introduced by Rokhlin [5] as a fast solution method for integral equations for two-dimensional Laplace equation. Barnes and Hut [6] developed an O
(
N
log
⁡
N
) algorithm. The FMM was then developed by Greengard and Rokhlin [1–3], as a fast evaluation method for the pairwise force calculation in multibody problems with Coulombic potential. FMM has been applied to problems in various fields such as boundary integral equation method (BIEM) and molecular dynamics (MD). Rokhlin [7] and Coifman et al. [8] used the FMM to solve the Helmholtz equation in 2D. Greengard et al. [9] and Mammoli and Ingber [10, 11] applied the FMM to the Stokes equations in fluid mechanic and Greengard and Helsing [12] for the problem of elasticity. Nishimura et al. [13, 14] studied this topic in order to apply fast multipole boundary integral equation method (FM-BIEM) to practical problems in fracture mechanics and earthquake engineering. Gumerov and Duraiswami used the FMM to solve the biharmonic equation [15]. In this paper, we discuss the applications of the FMM-BIEM to the fundamental boundary value problems in three dimensions given by Stokes equations and show its efficiencies.


## 3. The FMM Principles for the 3D Stokes Problem

The FMM method replaces the classical product matrix vector by an operation called the multipole product that allows obtaining, by a very fast algorithm, a good approximation of the solution. The complexity is ofO
(
N
log
⁡
N
) or O
(
N
) in place of O
(
N
2
) for the direct calculus. The basic idea of the multipole method is the separation of the variables in the Green function. This Green function is rewritten differently in a new expression that determines the form of the moments and the transfer functions. To develop the FMM method to the kind of problems of our interest, we will consider the integral equation of the problem to be solved for the unknown values of the intensities of Stokeslets (see Pozrikidis [16]) located on the surfaces of the particles and will consider that this equation has been already discretized by the boundary element method. The problem is then supposed to be reduced to a linear system, of order N, of the unknown values of stress at the collocation points of the particles surfaces. N is related to the number N
p of the particles by N
=
3
N
p
N
T, where N
T is the number of the triangles per particle. The linear system is solved by an iterative method. At each iteration, we start from a great imaginary cube containing the whole collocations points. This cube constitutes level zero of the subdivision. The n
th subdivision is obtained from the (
n
-
1
)
th one by subdividing each cube into 8 small cubes (cf. Figure 1).Figure 1
Recursive subdivision of the initial cube.The level of the constructed tree is chosen after fixing a desired accuracy of computationɛ. The set of collocation points is then distributed into several subdomains. For each one a moment centered at the middle of the subdomain is calculated as well as the transfer function for the far field interaction. To evaluate the integral equation at a given point x using the O
(
N
log
⁡
N
) scheme, we have to work in a tree given by the recursive subdivision of the cubic initial domain and evaluate the far field interaction when necessary. The near interactions are evaluated at the finest level. If we have to use the O
(
N
), we must translate the moment of each box at each level to the far field box and convert the information of the far field subdomain at a local moment centered at the middle of the box. Beginning at the coarsest level, these local moments are shifted to the children level (cf. Figure 2).Figure 2
Transmission of the far field expansion to the boxes in its interaction list and translation to the children.To evaluate the integral at a given pointx, we go directly to the box containing this point at the finest level and calculate the far field expansion due to all points outside the neighbors boxes. The near interactions are evaluated directly in a similar manner as in O
(
N
log
⁡
N
) scheme.


## 4. The Harmonic Functions and Laplace Equation

It is well known that the fundamental solution of the Laplace equation referred to as the Green function of the Laplacian operator can be written using the spherical harmonicY
n
m [2]. The potential field due to a source point q located at y
(
ρ
,
α
,
β
) and calculated at a point x
(
r
,
θ
,
ϕ
) is given by
(1)
Φ
(
x
,
y
)
=
1
|
x
-
y
|
=
∑
n
=
0
∞
‍
∑
m
=
-
n
n
‍
M
n
m
r
n
+
1
Y
n
m
(
θ
,
ϕ
)
,
where
(2)
M
n
m
=
q
ρ
n
Y
n
-
m
(
α
,
β
)
,
(3)
Y
n
m
(
θ
,
ϕ
)
=
(
n
-
|
m
|
)
!
(
n
+
|
m
|
)
!
P
n
m
(
cos
⁡
(
θ
)
)
e
i
m
ϕ
.
This potential can be written differently using the solid harmonics R
n
m and S
n
m [4]:
(4)
Φ
(
x
,
y
)
=
∑
n
=
0
∞
‍
∑
m
=
-
n
n
‍
R
n
m
(
o
y
→
)
¯
S
n
m
(
o
x
→
)
=
∑
n
=
0
∞
‍
∑
m
=
-
n
n
‍
R
n
m
(
o
y
→
)
S
n
m
(
o
x
→
)
¯
.
In this paper, we will use new functions that we define from spherical harmonic as follows:
(5)
S
y
n
m
=
α
n
m
1
r
n
+
1
Y
n
m
(
θ
,
ϕ
)
,
R
y
n
m
=
1
α
n
m
ρ
n
Y
n
m
(
α
,
β
)
,
where
(6)
α
n
m
=
(
n
+
|
m
|
)
!
(
n
-
|
m
|
)
!
.
From (3) and (5), we obtain the following expression for the functions introduced above:
(7)
S
y
n
m
(
o
x
→
)
=
(
n
-
|
m
|
)
!
P
n
|
m
|
(
cos
⁡
θ
)
e
i
m
θ
1
r
n
+
1
,
R
y
n
m
(
o
y
→
)
=
1
(
n
+
|
m
|
)
!
P
n
|
m
|
(
cos
⁡
α
)
e
i
m
β
ρ
n
.
Here, the special functions P
n
m are the associated Legendre functions, defined in [17] as
(8)
P
n
m
(
x
)
=
(
-
1
)
m
(
1
-
x
2
)
m
/
2
d
m
d
x
m
P
n
(
x
)
,
where P
n
(
x
) denotes the Legendre polynomial of degree n, defined by Rodrigues formula [17].The first function coincides with solid harmonicS
n
m for m
≥
0, and the second one coincides with solid harmonic S
n
m for m
≥
0.For the negative values ofm, we have the following propriety:
(9)
S
n
-
m
(
o
x
→
)
=
(
-
1
)
m
S
n
m
(
o
x
→
)
¯
;
R
n
-
m
(
o
x
→
)
=
(
-
1
)
m
R
n
m
(
o
x
→
)
¯
,
S
y
n
-
m
(
o
x
→
)
=
S
y
n
m
(
o
x
→
)
¯
;
R
y
n
-
m
(
o
x
→
)
=
R
y
n
m
(
o
x
→
)
¯
S
y
n
-
m
(
o
x
→
)
=
(
-
1
)
m
S
n
-
m
(
o
x
→
)
;
R
y
n
-
m
(
o
y
→
)
=
(
-
1
)
m
R
n
-
m
(
o
y
→
)
;
m
>
0
.


### 4.1. Representation of the Laplacian Kernel

Like in (1) we have the following new representation of the Green function of the Laplace equation:
(10)
Φ
(
x
,
y
)
=
∑
n
=
0
∞
‍
∑
m
=
-
n
n
‍
R
y
n
m
(
o
y
→
)
¯
S
y
n
m
(
o
x
→
)
=
∑
n
=
0
∞
∑
m
=
-
n
n
R
y
n
m
(
o
y
→
)
S
y
n
m
(
o
x
→
)
¯
.
Let us suppose that the point y is located inside a sphere D of radius a with center o and x is outside the sphere D. If we approach Φ by a finite sum obtained from (10) at the order p, we obtain an error bound which is the same error obtained by truncating the sums in (1) and (4):
(11)
E
p
=
1
r
-
a
(
a
r
)
p
+
1
.
The major advantage of this new representation is that we only have to work with positive index m, which decreases the complexity of the algorithms in a remarkable way.


## 4.1. Representation of the Laplacian Kernel

Like in (1) we have the following new representation of the Green function of the Laplace equation:
(10)
Φ
(
x
,
y
)
=
∑
n
=
0
∞
‍
∑
m
=
-
n
n
‍
R
y
n
m
(
o
y
→
)
¯
S
y
n
m
(
o
x
→
)
=
∑
n
=
0
∞
∑
m
=
-
n
n
R
y
n
m
(
o
y
→
)
S
y
n
m
(
o
x
→
)
¯
.
Let us suppose that the point y is located inside a sphere D of radius a with center o and x is outside the sphere D. If we approach Φ by a finite sum obtained from (10) at the order p, we obtain an error bound which is the same error obtained by truncating the sums in (1) and (4):
(11)
E
p
=
1
r
-
a
(
a
r
)
p
+
1
.
The major advantage of this new representation is that we only have to work with positive index m, which decreases the complexity of the algorithms in a remarkable way.


## 5. The Schemes for 3D Stokes Problem

We will develop here the two numerical schemes of theN
log
⁡
(
N
) and the N order multipole expansion associated with the problem of several solid arbitrarily shaped particles moving in a creeping flow, which are based on a suitable translation formulas of S
y
n
m and R
y
n
m.The fundamental solution or the Green function of the free-space Stokes problem [18] is
(12)
G
k
j
(
x
,
y
)
=
-
1
8
π
ν
[
δ
k
j
|
x
-
y
|
+
(
x
j
-
y
j
)
(
x
k
-
y
k
)
|
x
-
y
|
3
]
.
Except for the multiplicative constant, the Green function may be written as follows:
(13)
G
k
j
(
x
,
y
)
=
[
δ
k
j
1
|
x
-
y
|
-
(
o
x
→
)
k
∂
∂
x
j
(
1
|
x
-
y
|
)
vvvvivivvvv
+
(
o
y
→
)
k
∂
∂
x
j
(
1
|
x
-
y
|
)
]
,
and the integral equation for the velocity on a surface S
l of a particle l is
(14)
v
k
(
x
)
=
∫
S
‍
G
k
j
(
x
,
y
)
q
j
(
y
)
d
s
;
x
∈
S
l
,
where S
=
∪
i
=
1
N
p
S
i, N
p is the number of particles, S
i is the surface of the ith particle, and q
j
(
y
) is jth component of the surface stress.The discrete problem associated with (14) is presented as follows:
(15)
v
k
(
x
G
l
′
,
s
′
)
=
∑
l
=
1
N
p
‍
∑
s
=
1
N
t
‍
(
∫
T
s
l
‍
G
k
j
(
x
G
l
′
,
s
′
,
y
)
d
s
(
y
)
)
q
j
s
,
l
,
where T
s
l is the surface element, x
G
l
′
,
s
′ is the gravity center of a given surface, N
p and N
t denote, respectively, the number of particle and surface element per particle, and q
j
s
,
l is the stress applied on the surface element. Since we are interested here in the far field calculus, the treatment of the singularity problem will not be considered.


### 5.1.O
(
N
log
⁡
N
) Scheme

Let us suppose that we have uniform distribution points on a certain surfaceS located inside a sphere D of radius a with center o. From (10) and (13), the velocity at a point x outside the sphere D distant of r from o is given by the following equation:
(16)
v
k
(
x
)
=
∫
S
‍
G
k
j
(
x
,
y
)
q
j
(
y
)
d
s
=
∑
n
=
0
∞
∑
m
=
-
n
n
[
F
y
n
m
j
k
(
o
x
→
)
¯
M
n
,
j
m
(
1
)
(
o
)
v
v
z
v
v
v
v
v
v
i
v
+
H
y
n
m
j
(
o
x
→
)
¯
M
n
,
j
m
,
k
(
2
)
(
o
)
]
,
where
(17)
M
n
,
j
m
,
(
1
)
(
o
)
=
∫
s
‍
R
y
n
m
(
o
y
→
)
q
j
(
y
)
d
s
,
M
n
,
j
m
,
k
(
2
)
(
o
)
=
∫
s
‍
(
o
y
→
k
)
R
y
n
m
(
o
y
→
)
q
j
(
y
)
d
s
,
F
y
n
m
j
k
(
o
x
→
)
=
δ
k
j
S
y
n
m
(
o
x
→
)
-
(
o
x
→
)
k
∂
∂
x
j
(
S
y
n
m
(
o
x
→
)
)
,
H
y
n
,
m
j
o
x
→
=
∂
∂
x
j
(
S
y
n
m
(
o
x
→
)
)
.
One can notice that, for a given pair of n and m
,
M
n
,
j
m
,
k
(
2
)
(
o
) has 9 components and M
n
,
j
m
(
1
)
(
o
) has 3 components, and therefore the number of multipole moments in this formulation is 12. The derivatives that appear in the relations above are calculated in an exact way using the following formulas.


#### 5.1.1. Derivatives Formulas

According to the relation between Cartesian coordinates and polar coordinates and some chain rules [19], the derivatives of S
y
n
m can be obtained as follows:
(18)
∂
∂
x
1
(
S
y
n
m
(
o
x
→
)
)
=
{
1
2
(
S
y
n
+
1
m
-
1
(
o
x
→
)
-
S
y
n
+
1
m
+
1
(
o
x
→
)
)
,
if
m
>
0
;
-
Real
(
S
y
n
+
1
1
(
o
x
→
)
)
,
if
m
=
0
;
∂
∂
x
1
(
S
y
n
-
m
(
o
x
→
)
)
¯
,
if
m
<
0
.
∂
∂
x
2
(
S
y
n
m
(
o
x
→
)
)
=
{
i
2
(
S
y
n
+
1
m
-
1
(
o
x
→
)
+
S
y
n
+
1
m
+
1
(
o
x
→
)
)
,
if
m
>
0
;
i
Imag
(
S
y
n
+
1
1
(
o
x
→
)
)
,
if
m
=
0
;
∂
∂
x
2
(
S
y
n
-
m
(
o
x
→
)
)
¯
,
if
m
<
0
.
∂
∂
x
3
(
S
y
n
m
(
o
x
→
)
)
=
{
-
S
y
n
+
1
m
(
o
x
→
)
,
if
m
>
0
;
∂
∂
x
3
(
S
y
n
-
m
(
o
x
→
)
)
¯
,
else
.
In a similar manner, we can obtain the derivatives formulas of R
y
n
m for n
≥
1 as follows:
(19)
∂
∂
x
1
(
R
y
n
m
(
o
x
→
)
)
=
{
1
2
(
R
y
n
+
1
m
-
1
(
o
x
→
)
-
R
y
n
+
1
m
+
1
(
o
x
→
)
)
,
if
m
>
0
;
-
Real
(
R
y
n
+
1
1
(
o
x
→
)
)
,
if
m
=
0
;
∂
∂
x
1
(
R
y
n
-
m
(
o
x
→
)
)
¯
,
if
m
<
0
,
∂
∂
x
2
(
R
y
n
m
(
o
x
→
)
)
=
{
i
2
(
R
y
n
+
1
m
-
1
(
o
x
→
)
+
R
y
n
+
1
m
+
1
(
o
x
→
)
)
,
if
m
>
0
;
i
Imag
(
R
y
n
+
1
1
(
o
x
→
)
)
,
if
m
=
0
;
∂
∂
x
2
(
R
y
n
-
m
(
o
x
→
)
)
¯
,
if
m
<
0
,
∂
∂
x
3
(
R
y
n
m
(
o
x
→
)
)
=
{
-
R
y
n
+
1
m
(
o
x
→
)
,
if
m
>
0
;
∂
∂
x
3
(
R
y
n
-
m
(
o
x
→
)
)
¯
,
else
,
where i
2
=
-
1. It should be noted that (
∂
/
∂
x
j
)
(
R
y
0
0
(
o
x
→
)
)
=
0 for all j
∈
{
1,2
,
3
}.


#### 5.1.2. Error Bound

By truncating (16) at the order p, and basing on the error given by (11), the error bound is given by
(20)
E
p
′
=
C
(
δ
j
k
r
-
a
+
1
r
-
a
+
a
r
2
-
r
a
)
(
a
r
)
p
+
1
,
where C
=
∫
S
‍
|
q
(
y
)
|
d
s.


### 5.2.O
(
N
) Scheme

To develop a scheme of orderN, we introduce three basic translations adequate to the functions S
y
n
m and R
y
n
m which will be useful to build the formulas of translations multipole to multipole (M2M), multipole to local (M2L), and local to local (L2L). For the rest of the paper, we will take k
=
inf
⁡
(
m
,
n
′
).


#### 5.2.1. Relations Formulas betweenS
y
n
m and R
y
n
m

Under the condition|
o
y
→
|
<
|
o
o
′
→
|, the function R
y
n
m satisfies a relationship of the following form:
(21)
R
y
n
m
(
o
′
y
→
)
=
∑
n
′
=
0
n
‍
{
∑
m
′
=
0
n
′
‍
(
-
1
)
m
′
R
y
n
′
m
′
(
o
′
o
→
)
¯
R
y
n
-
n
′
m
+
m
′
(
o
y
→
)
∑
m
′
=
0
k
‍
R
y
n
′
m
′
(
o
′
o
→
)
R
y
n
-
n
′
m
-
m
′
(
o
y
→
)
vvvvvivv
+
∑
m
′
=
k
+
1
n
′
‍
(
-
1
)
m
′
-
m
R
y
n
′
m
′
(
o
′
o
→
)
R
y
n
-
n
′
m
′
-
m
(
o
y
→
)
¯
vvvvivvv
+
∑
m
′
=
0
n
′
‍
(
-
1
)
m
′
R
y
n
′
m
′
(
o
′
o
→
)
¯
R
y
n
-
n
′
m
+
m
′
(
o
y
→
)
}
which allows a pole shift that will permit us to construct the multipolar to multipolar translation.Moreover, and under the condition|
x
0
x
→
|
<
|
o
x
0
→
|, the function S
y
n
m can be further expressed in the following form that will be used later to construct the conversion multipolar to local moment:
(22)
S
y
n
m
(
o
x
→
)
=
∑
n
′
=
0
∞
‍
{
∑
m
′
=
0
n
′
‍
(
-
1
)
n
′
R
y
n
′
m
′
(
x
0
x
→
)
S
y
n
+
n
′
m
+
m
′
(
o
x
0
→
)
vvvvvivv
+
∑
m
′
=
0
k
‍
(
-
1
)
n
′
+
m
′
R
y
n
′
m
′
(
x
0
x
→
)
¯
S
y
n
+
n
′
m
-
m
′
(
o
x
0
→
)
vvvvvivv
+
∑
m
′
=
k
+
1
n
′
‍
(
-
1
)
n
+
2
m
′
-
m
R
y
n
′
m
′
(
x
0
x
→
)
¯
S
y
n
+
n
′
m
′
-
m
(
o
x
0
→
)
}
.
Just like the first relation, under the condition |
x
0
x
→
|
<
|
x
0
x
1
→
|, we have a third equation that satisfies a pole relationship that will be used later to construct the local to local translation:
(23)
R
y
n
m
(
x
0
x
→
)
=
∑
n
′
=
n
∞
‍
{
∑
m
′
=
0
n
′
‍
(
-
1
)
m
′
R
y
n
′
m
′
(
x
1
x
→
)
¯
R
y
n
-
n
′
m
+
m
′
(
x
0
x
1
→
)
∑
m
′
=
0
k
‍
R
y
n
′
m
′
(
x
1
x
→
)
R
y
n
-
n
′
m
-
m
′
(
x
0
x
1
→
)
vvvvvivv
+
∑
m
′
=
k
+
1
n
′
‍
(
-
1
)
m
′
-
m
R
y
n
′
m
′
(
x
1
x
→
)
R
y
n
-
n
′
m
′
-
m
(
x
0
x
1
→
)
¯
vvvvvivv
+
∑
m
′
=
0
n
′
‍
(
-
1
)
m
′
R
y
n
′
m
′
(
x
1
x
→
)
¯
R
y
n
-
n
′
m
+
m
′
(
x
0
x
1
→
)
}
.


#### 5.2.2. Translation of a Multipole Expansion

Let us suppose that the surfaceS contains uniform distribution points included inside a sphere D of radius a centered at a point o at a distance ρ from the origin o
′ with ρ
>
(
c
+
1
)
a and c
>
1. The integral equation is given by the multipole expansion (16):
(24)
v
k
(
x
)
=
∑
n
=
0
∞
‍
∑
m
=
-
n
n
‍
[
F
y
n
m
j
k
(
o
x
→
)
¯
M
n
,
j
m
(
1
)
(
o
)
vvvvvivv
+
H
y
n
m
j
(
o
x
→
)
¯
M
j
,
m
m
,
k
(
2
)
(
o
)
]
.
If the point x is outside the sphere D
0 of radius a
+
ρ, then the integral equation takes the following form:
(25)
v
k
(
x
)
=
∑
n
=
0
∞
‍
∑
m
=
-
n
n
[
F
y
n
m
j
k
(
o
′
x
→
)
¯
M
n
,
j
m
(
1
)
(
o
′
)
‍
vvvvvivv
+
H
y
n
m
j
(
o
′
x
→
)
¯
M
j
,
m
m
,
k
(
2
)
(
o
′
)
]
,
where
(26)
M
n
,
j
m
(
1
)
(
o
′
)
=
∑
n
′
=
0
n
‍
{
∑
m
′
=
1
n
′
‍
(
-
1
)
m
′
R
y
n
′
m
′
(
o
′
o
→
)
¯
M
n
-
n
′
,
j
m
+
m
′
(
1
)
(
o
)
∑
m
′
=
0
k
‍
R
y
n
′
m
′
(
o
′
o
→
)
M
n
-
n
′
,
j
m
-
m
′
(
1
)
(
o
)
vvvvvivv
+
∑
m
′
=
k
+
1
n
′
‍
(
-
1
)
m
′
-
m
R
y
n
′
m
′
(
o
′
o
→
)
M
n
-
n
′
,
j
m
′
-
m
(
1
)
(
o
)
¯
vvvvvivv
+
∑
m
′
=
1
n
′
‍
(
-
1
)
m
′
R
y
n
′
m
′
(
o
′
o
→
)
¯
M
n
-
n
′
,
j
m
+
m
′
(
1
)
(
o
)
}
,
M
n
,
j
m
,
k
(
2
)
(
o
′
)
=
(
o
′
o
→
)
k
M
n
,
j
m
(
1
)
(
o
′
)
+
∑
n
′
=
0
n
‍
{
∑
m
′
=
1
n
′
‍
(
-
1
)
m
′
R
y
n
′
m
′
(
o
′
o
→
)
¯
M
n
-
n
′
,
j
m
+
m
′
,
k
(
2
)
(
o
)
∑
m
′
=
0
k
‍
R
y
n
′
m
′
(
o
′
o
→
)
M
n
-
n
′
,
j
m
-
m
′
,
k
(
2
)
(
o
)
vvvvvvivv
+
∑
m
′
=
k
+
1
n
′
‍
(
-
1
)
m
′
-
m
R
y
n
′
m
′
(
o
′
o
→
)
M
n
-
n
′
,
j
m
′
-
m
,
k
(
2
)
(
o
)
¯
vvvvvvivv
+
∑
m
′
=
1
n
′
‍
(
-
1
)
m
′
R
y
n
′
m
′
(
o
′
o
→
)
¯
M
n
-
n
′
,
j
m
+
m
′
,
k
(
2
)
(
o
)
}
.


#### 5.2.3. Conversion of a Multipole Expansion into a Local Expansion

Let a surfaceS containing uniform distribution points included inside a sphere D of radius a centered at a point o at a distance ρ from a given point x
0 with ρ
>
(
c
+
1
)
a and c
>
1. The corresponding multipole expansion (16) converges for any points x inside the sphere D
0 of radius a centered at x
0. The velocity is given by
(27)
v
k
(
x
)
=
∫
S
‍
G
k
j
(
x
,
y
)
q
j
(
y
)
d
s
=
∑
n
=
0
∞
‍
∑
m
=
-
n
n
‍
[
F
y
R
n
m
j
k
(
x
0
x
→
)
¯
L
n
,
j
m
(
1
)
(
x
0
)
vvvvvwiivv
+
H
y
R
n
m
j
(
x
0
x
→
)
¯
L
n
,
j
m
,
k
(
2
)
(
x
0
)
]
,
where
(28)
L
n
,
j
m
(
1
)
(
x
0
)
=
∑
n
′
=
0
∞
{
∑
m
′
=
0
n
′
‍
(
-
1
)
n
M
n
′
,
j
m
′
(
1
)
(
o
)
S
y
n
+
n
′
m
+
m
′
(
o
x
0
→
)
¯
‍
vvvvvivv
+
∑
m
′
=
1
k
‍
(
-
1
)
n
+
m
′
M
n
′
,
j
m
′
(
1
)
(
o
)
¯
S
y
n
+
n
′
m
-
m
′
(
o
x
0
→
)
¯
vvvvvivv
+
∑
m
′
=
k
+
1
n
′
‍
(
-
1
)
n
-
m
M
n
′
,
j
m
′
(
1
)
(
o
)
¯
S
y
n
+
n
′
m
′
-
m
(
o
x
0
→
)
}
,
L
n
,
j
m
,
k
(
2
)
(
x
0
)
=
-
(
o
x
0
→
)
k
L
n
,
j
m
(
1
)
(
x
0
)
+
∑
n
′
=
0
∞
‍
{
∑
m
′
=
0
n
′
‍
(
-
1
)
n
M
n
′
,
j
m
′
,
k
(
2
)
(
o
)
S
y
n
+
n
′
m
+
m
′
(
o
x
0
→
)
¯
vvvvvvivv
+
∑
m
′
=
1
k
‍
(
-
1
)
n
+
m
′
M
n
′
,
j
m
′
,
k
(
2
)
(
o
)
¯
S
y
n
+
n
′
m
-
m
′
(
o
x
0
→
)
¯
vvvvvvivv
+
∑
m
′
=
k
+
1
n
′
‍
(
-
1
)
n
-
m
M
n
′
,
j
m
′
,
k
(
2
)
(
o
)
¯
S
y
n
+
n
′
m
′
-
m
(
o
x
0
→
)
}
,
F
y
R
n
m
j
k
(
x
0
x
→
)
=
δ
k
j
R
y
n
m
(
x
0
x
→
)
-
(
x
0
x
→
)
∂
∂
x
j
(
R
y
n
m
(
x
0
x
→
)
)
,
H
y
R
n
m
j
(
x
0
x
→
)
=
∂
∂
x
j
(
R
y
n
m
(
x
0
x
→
)
)
.


#### 5.2.4. Translation of a Local Expansion

Suppose thatv
(
x
) can be written as in (29), and suppose that x
1
∈
D
0; then for any point x
∈
D
0, v
(
x
) can be written according to the local moment centered at x
1 as follows:
(29)
v
k
(
x
)
=
∫
S
‍
G
k
j
(
x
,
y
)
q
j
(
y
)
d
s
=
∑
n
=
0
∞
‍
∑
m
=
-
n
n
‍
[
F
y
R
n
m
j
k
(
x
1
x
→
)
¯
L
n
,
j
m
(
1
)
(
x
1
)
vvvvvwiivv
+
H
y
R
n
m
j
(
x
1
x
→
)
¯
L
n
,
j
m
,
k
(
2
)
(
x
1
)
]
,
where
(30)
L
n
,
j
m
,
k
(
1
)
(
x
1
)
=
∑
n
′
=
0
∞
‍
{
∑
m
′
=
1
n
′
‍
(
-
1
)
m
L
n
′
,
j
m
′
(
1
)
(
x
0
)
R
y
n
-
n
′
m
+
m
′
(
x
0
x
1
→
)
¯
∑
m
′
=
0
k
‍
(
-
1
)
m
-
m
′
L
n
′
,
j
m
′
(
1
)
(
x
0
)
R
y
n
-
n
′
m
-
m
′
(
x
0
x
1
→
)
¯
vvvvvivv
+
∑
m
′
=
k
+
1
n
′
‍
L
n
′
,
j
m
′
(
1
)
(
x
0
)
R
y
n
-
n
′
m
′
-
m
(
x
0
x
1
→
)
vvvvvivv
+
∑
m
′
=
1
n
′
‍
(
-
1
)
m
L
n
′
,
j
m
′
(
1
)
(
x
0
)
¯
R
y
n
-
n
′
m
+
m
′
(
x
0
x
1
→
)
}
,
L
n
,
j
m
,
k
(
2
)
(
x
1
)
=
-
(
x
0
x
1
→
)
k
L
n
,
j
m
(
1
)
(
x
1
)
+
∑
n
′
=
0
∞
‍
{
∑
m
′
=
1
n
′
‍
(
-
1
)
m
L
n
′
,
j
m
′
,
k
(
2
)
(
x
0
)
¯
R
y
n
-
n
′
m
+
m
′
(
x
0
x
1
→
)
¯
∑
m
′
=
0
k
‍
(
-
1
)
m
-
m
′
L
n
′
,
j
m
′
,
k
(
2
)
(
x
0
)
R
y
n
-
n
′
m
-
m
′
(
x
0
x
1
→
)
¯
vvvvvvivv
+
∑
m
′
=
k
+
1
n
′
‍
L
n
′
,
j
m
′
,
k
(
2
)
(
x
0
)
R
y
n
-
n
′
m
′
-
m
(
x
0
x
1
→
)
vvvvvvivv
+
∑
m
′
=
1
n
′
‍
(
-
1
)
m
L
n
′
,
j
m
′
,
k
(
2
)
(
x
0
)
¯
R
y
n
-
n
′
m
+
m
′
(
x
0
x
1
→
)
¯
}
.


#### 5.2.5. Error Bound

With a similar reasoning as in the error bound of theO
(
N
log
⁡
N
) scheme and under the same hypothesis of the representation (29), the evaluation of the velocity v
(
x
) at an order p generates an error bound expressed as follows:
(31)
E
p
′′
=
C
(
δ
j
k
ρ
-
r
)
+
(
1
ρ
-
r
)
+
(
r
ρ
2
-
ρ
r
)
(
r
ρ
)
p
+
1
,
where C
=
∫
C
‍
|
q
(
y
)
|
d
s.


## 5.1.O
(
N
log
⁡
N
) Scheme

Let us suppose that we have uniform distribution points on a certain surfaceS located inside a sphere D of radius a with center o. From (10) and (13), the velocity at a point x outside the sphere D distant of r from o is given by the following equation:
(16)
v
k
(
x
)
=
∫
S
‍
G
k
j
(
x
,
y
)
q
j
(
y
)
d
s
=
∑
n
=
0
∞
∑
m
=
-
n
n
[
F
y
n
m
j
k
(
o
x
→
)
¯
M
n
,
j
m
(
1
)
(
o
)
v
v
z
v
v
v
v
v
v
i
v
+
H
y
n
m
j
(
o
x
→
)
¯
M
n
,
j
m
,
k
(
2
)
(
o
)
]
,
where
(17)
M
n
,
j
m
,
(
1
)
(
o
)
=
∫
s
‍
R
y
n
m
(
o
y
→
)
q
j
(
y
)
d
s
,
M
n
,
j
m
,
k
(
2
)
(
o
)
=
∫
s
‍
(
o
y
→
k
)
R
y
n
m
(
o
y
→
)
q
j
(
y
)
d
s
,
F
y
n
m
j
k
(
o
x
→
)
=
δ
k
j
S
y
n
m
(
o
x
→
)
-
(
o
x
→
)
k
∂
∂
x
j
(
S
y
n
m
(
o
x
→
)
)
,
H
y
n
,
m
j
o
x
→
=
∂
∂
x
j
(
S
y
n
m
(
o
x
→
)
)
.
One can notice that, for a given pair of n and m
,
M
n
,
j
m
,
k
(
2
)
(
o
) has 9 components and M
n
,
j
m
(
1
)
(
o
) has 3 components, and therefore the number of multipole moments in this formulation is 12. The derivatives that appear in the relations above are calculated in an exact way using the following formulas.


### 5.1.1. Derivatives Formulas

According to the relation between Cartesian coordinates and polar coordinates and some chain rules [19], the derivatives of S
y
n
m can be obtained as follows:
(18)
∂
∂
x
1
(
S
y
n
m
(
o
x
→
)
)
=
{
1
2
(
S
y
n
+
1
m
-
1
(
o
x
→
)
-
S
y
n
+
1
m
+
1
(
o
x
→
)
)
,
if
m
>
0
;
-
Real
(
S
y
n
+
1
1
(
o
x
→
)
)
,
if
m
=
0
;
∂
∂
x
1
(
S
y
n
-
m
(
o
x
→
)
)
¯
,
if
m
<
0
.
∂
∂
x
2
(
S
y
n
m
(
o
x
→
)
)
=
{
i
2
(
S
y
n
+
1
m
-
1
(
o
x
→
)
+
S
y
n
+
1
m
+
1
(
o
x
→
)
)
,
if
m
>
0
;
i
Imag
(
S
y
n
+
1
1
(
o
x
→
)
)
,
if
m
=
0
;
∂
∂
x
2
(
S
y
n
-
m
(
o
x
→
)
)
¯
,
if
m
<
0
.
∂
∂
x
3
(
S
y
n
m
(
o
x
→
)
)
=
{
-
S
y
n
+
1
m
(
o
x
→
)
,
if
m
>
0
;
∂
∂
x
3
(
S
y
n
-
m
(
o
x
→
)
)
¯
,
else
.
In a similar manner, we can obtain the derivatives formulas of R
y
n
m for n
≥
1 as follows:
(19)
∂
∂
x
1
(
R
y
n
m
(
o
x
→
)
)
=
{
1
2
(
R
y
n
+
1
m
-
1
(
o
x
→
)
-
R
y
n
+
1
m
+
1
(
o
x
→
)
)
,
if
m
>
0
;
-
Real
(
R
y
n
+
1
1
(
o
x
→
)
)
,
if
m
=
0
;
∂
∂
x
1
(
R
y
n
-
m
(
o
x
→
)
)
¯
,
if
m
<
0
,
∂
∂
x
2
(
R
y
n
m
(
o
x
→
)
)
=
{
i
2
(
R
y
n
+
1
m
-
1
(
o
x
→
)
+
R
y
n
+
1
m
+
1
(
o
x
→
)
)
,
if
m
>
0
;
i
Imag
(
R
y
n
+
1
1
(
o
x
→
)
)
,
if
m
=
0
;
∂
∂
x
2
(
R
y
n
-
m
(
o
x
→
)
)
¯
,
if
m
<
0
,
∂
∂
x
3
(
R
y
n
m
(
o
x
→
)
)
=
{
-
R
y
n
+
1
m
(
o
x
→
)
,
if
m
>
0
;
∂
∂
x
3
(
R
y
n
-
m
(
o
x
→
)
)
¯
,
else
,
where i
2
=
-
1. It should be noted that (
∂
/
∂
x
j
)
(
R
y
0
0
(
o
x
→
)
)
=
0 for all j
∈
{
1,2
,
3
}.


### 5.1.2. Error Bound

By truncating (16) at the order p, and basing on the error given by (11), the error bound is given by
(20)
E
p
′
=
C
(
δ
j
k
r
-
a
+
1
r
-
a
+
a
r
2
-
r
a
)
(
a
r
)
p
+
1
,
where C
=
∫
S
‍
|
q
(
y
)
|
d
s.


## 5.1.1. Derivatives Formulas

According to the relation between Cartesian coordinates and polar coordinates and some chain rules [19], the derivatives of S
y
n
m can be obtained as follows:
(18)
∂
∂
x
1
(
S
y
n
m
(
o
x
→
)
)
=
{
1
2
(
S
y
n
+
1
m
-
1
(
o
x
→
)
-
S
y
n
+
1
m
+
1
(
o
x
→
)
)
,
if
m
>
0
;
-
Real
(
S
y
n
+
1
1
(
o
x
→
)
)
,
if
m
=
0
;
∂
∂
x
1
(
S
y
n
-
m
(
o
x
→
)
)
¯
,
if
m
<
0
.
∂
∂
x
2
(
S
y
n
m
(
o
x
→
)
)
=
{
i
2
(
S
y
n
+
1
m
-
1
(
o
x
→
)
+
S
y
n
+
1
m
+
1
(
o
x
→
)
)
,
if
m
>
0
;
i
Imag
(
S
y
n
+
1
1
(
o
x
→
)
)
,
if
m
=
0
;
∂
∂
x
2
(
S
y
n
-
m
(
o
x
→
)
)
¯
,
if
m
<
0
.
∂
∂
x
3
(
S
y
n
m
(
o
x
→
)
)
=
{
-
S
y
n
+
1
m
(
o
x
→
)
,
if
m
>
0
;
∂
∂
x
3
(
S
y
n
-
m
(
o
x
→
)
)
¯
,
else
.
In a similar manner, we can obtain the derivatives formulas of R
y
n
m for n
≥
1 as follows:
(19)
∂
∂
x
1
(
R
y
n
m
(
o
x
→
)
)
=
{
1
2
(
R
y
n
+
1
m
-
1
(
o
x
→
)
-
R
y
n
+
1
m
+
1
(
o
x
→
)
)
,
if
m
>
0
;
-
Real
(
R
y
n
+
1
1
(
o
x
→
)
)
,
if
m
=
0
;
∂
∂
x
1
(
R
y
n
-
m
(
o
x
→
)
)
¯
,
if
m
<
0
,
∂
∂
x
2
(
R
y
n
m
(
o
x
→
)
)
=
{
i
2
(
R
y
n
+
1
m
-
1
(
o
x
→
)
+
R
y
n
+
1
m
+
1
(
o
x
→
)
)
,
if
m
>
0
;
i
Imag
(
R
y
n
+
1
1
(
o
x
→
)
)
,
if
m
=
0
;
∂
∂
x
2
(
R
y
n
-
m
(
o
x
→
)
)
¯
,
if
m
<
0
,
∂
∂
x
3
(
R
y
n
m
(
o
x
→
)
)
=
{
-
R
y
n
+
1
m
(
o
x
→
)
,
if
m
>
0
;
∂
∂
x
3
(
R
y
n
-
m
(
o
x
→
)
)
¯
,
else
,
where i
2
=
-
1. It should be noted that (
∂
/
∂
x
j
)
(
R
y
0
0
(
o
x
→
)
)
=
0 for all j
∈
{
1,2
,
3
}.


## 5.1.2. Error Bound

By truncating (16) at the order p, and basing on the error given by (11), the error bound is given by
(20)
E
p
′
=
C
(
δ
j
k
r
-
a
+
1
r
-
a
+
a
r
2
-
r
a
)
(
a
r
)
p
+
1
,
where C
=
∫
S
‍
|
q
(
y
)
|
d
s.


## 5.2.O
(
N
) Scheme

To develop a scheme of orderN, we introduce three basic translations adequate to the functions S
y
n
m and R
y
n
m which will be useful to build the formulas of translations multipole to multipole (M2M), multipole to local (M2L), and local to local (L2L). For the rest of the paper, we will take k
=
inf
⁡
(
m
,
n
′
).


### 5.2.1. Relations Formulas betweenS
y
n
m and R
y
n
m

Under the condition|
o
y
→
|
<
|
o
o
′
→
|, the function R
y
n
m satisfies a relationship of the following form:
(21)
R
y
n
m
(
o
′
y
→
)
=
∑
n
′
=
0
n
‍
{
∑
m
′
=
0
n
′
‍
(
-
1
)
m
′
R
y
n
′
m
′
(
o
′
o
→
)
¯
R
y
n
-
n
′
m
+
m
′
(
o
y
→
)
∑
m
′
=
0
k
‍
R
y
n
′
m
′
(
o
′
o
→
)
R
y
n
-
n
′
m
-
m
′
(
o
y
→
)
vvvvvivv
+
∑
m
′
=
k
+
1
n
′
‍
(
-
1
)
m
′
-
m
R
y
n
′
m
′
(
o
′
o
→
)
R
y
n
-
n
′
m
′
-
m
(
o
y
→
)
¯
vvvvivvv
+
∑
m
′
=
0
n
′
‍
(
-
1
)
m
′
R
y
n
′
m
′
(
o
′
o
→
)
¯
R
y
n
-
n
′
m
+
m
′
(
o
y
→
)
}
which allows a pole shift that will permit us to construct the multipolar to multipolar translation.Moreover, and under the condition|
x
0
x
→
|
<
|
o
x
0
→
|, the function S
y
n
m can be further expressed in the following form that will be used later to construct the conversion multipolar to local moment:
(22)
S
y
n
m
(
o
x
→
)
=
∑
n
′
=
0
∞
‍
{
∑
m
′
=
0
n
′
‍
(
-
1
)
n
′
R
y
n
′
m
′
(
x
0
x
→
)
S
y
n
+
n
′
m
+
m
′
(
o
x
0
→
)
vvvvvivv
+
∑
m
′
=
0
k
‍
(
-
1
)
n
′
+
m
′
R
y
n
′
m
′
(
x
0
x
→
)
¯
S
y
n
+
n
′
m
-
m
′
(
o
x
0
→
)
vvvvvivv
+
∑
m
′
=
k
+
1
n
′
‍
(
-
1
)
n
+
2
m
′
-
m
R
y
n
′
m
′
(
x
0
x
→
)
¯
S
y
n
+
n
′
m
′
-
m
(
o
x
0
→
)
}
.
Just like the first relation, under the condition |
x
0
x
→
|
<
|
x
0
x
1
→
|, we have a third equation that satisfies a pole relationship that will be used later to construct the local to local translation:
(23)
R
y
n
m
(
x
0
x
→
)
=
∑
n
′
=
n
∞
‍
{
∑
m
′
=
0
n
′
‍
(
-
1
)
m
′
R
y
n
′
m
′
(
x
1
x
→
)
¯
R
y
n
-
n
′
m
+
m
′
(
x
0
x
1
→
)
∑
m
′
=
0
k
‍
R
y
n
′
m
′
(
x
1
x
→
)
R
y
n
-
n
′
m
-
m
′
(
x
0
x
1
→
)
vvvvvivv
+
∑
m
′
=
k
+
1
n
′
‍
(
-
1
)
m
′
-
m
R
y
n
′
m
′
(
x
1
x
→
)
R
y
n
-
n
′
m
′
-
m
(
x
0
x
1
→
)
¯
vvvvvivv
+
∑
m
′
=
0
n
′
‍
(
-
1
)
m
′
R
y
n
′
m
′
(
x
1
x
→
)
¯
R
y
n
-
n
′
m
+
m
′
(
x
0
x
1
→
)
}
.


### 5.2.2. Translation of a Multipole Expansion

Let us suppose that the surfaceS contains uniform distribution points included inside a sphere D of radius a centered at a point o at a distance ρ from the origin o
′ with ρ
>
(
c
+
1
)
a and c
>
1. The integral equation is given by the multipole expansion (16):
(24)
v
k
(
x
)
=
∑
n
=
0
∞
‍
∑
m
=
-
n
n
‍
[
F
y
n
m
j
k
(
o
x
→
)
¯
M
n
,
j
m
(
1
)
(
o
)
vvvvvivv
+
H
y
n
m
j
(
o
x
→
)
¯
M
j
,
m
m
,
k
(
2
)
(
o
)
]
.
If the point x is outside the sphere D
0 of radius a
+
ρ, then the integral equation takes the following form:
(25)
v
k
(
x
)
=
∑
n
=
0
∞
‍
∑
m
=
-
n
n
[
F
y
n
m
j
k
(
o
′
x
→
)
¯
M
n
,
j
m
(
1
)
(
o
′
)
‍
vvvvvivv
+
H
y
n
m
j
(
o
′
x
→
)
¯
M
j
,
m
m
,
k
(
2
)
(
o
′
)
]
,
where
(26)
M
n
,
j
m
(
1
)
(
o
′
)
=
∑
n
′
=
0
n
‍
{
∑
m
′
=
1
n
′
‍
(
-
1
)
m
′
R
y
n
′
m
′
(
o
′
o
→
)
¯
M
n
-
n
′
,
j
m
+
m
′
(
1
)
(
o
)
∑
m
′
=
0
k
‍
R
y
n
′
m
′
(
o
′
o
→
)
M
n
-
n
′
,
j
m
-
m
′
(
1
)
(
o
)
vvvvvivv
+
∑
m
′
=
k
+
1
n
′
‍
(
-
1
)
m
′
-
m
R
y
n
′
m
′
(
o
′
o
→
)
M
n
-
n
′
,
j
m
′
-
m
(
1
)
(
o
)
¯
vvvvvivv
+
∑
m
′
=
1
n
′
‍
(
-
1
)
m
′
R
y
n
′
m
′
(
o
′
o
→
)
¯
M
n
-
n
′
,
j
m
+
m
′
(
1
)
(
o
)
}
,
M
n
,
j
m
,
k
(
2
)
(
o
′
)
=
(
o
′
o
→
)
k
M
n
,
j
m
(
1
)
(
o
′
)
+
∑
n
′
=
0
n
‍
{
∑
m
′
=
1
n
′
‍
(
-
1
)
m
′
R
y
n
′
m
′
(
o
′
o
→
)
¯
M
n
-
n
′
,
j
m
+
m
′
,
k
(
2
)
(
o
)
∑
m
′
=
0
k
‍
R
y
n
′
m
′
(
o
′
o
→
)
M
n
-
n
′
,
j
m
-
m
′
,
k
(
2
)
(
o
)
vvvvvvivv
+
∑
m
′
=
k
+
1
n
′
‍
(
-
1
)
m
′
-
m
R
y
n
′
m
′
(
o
′
o
→
)
M
n
-
n
′
,
j
m
′
-
m
,
k
(
2
)
(
o
)
¯
vvvvvvivv
+
∑
m
′
=
1
n
′
‍
(
-
1
)
m
′
R
y
n
′
m
′
(
o
′
o
→
)
¯
M
n
-
n
′
,
j
m
+
m
′
,
k
(
2
)
(
o
)
}
.


### 5.2.3. Conversion of a Multipole Expansion into a Local Expansion

Let a surfaceS containing uniform distribution points included inside a sphere D of radius a centered at a point o at a distance ρ from a given point x
0 with ρ
>
(
c
+
1
)
a and c
>
1. The corresponding multipole expansion (16) converges for any points x inside the sphere D
0 of radius a centered at x
0. The velocity is given by
(27)
v
k
(
x
)
=
∫
S
‍
G
k
j
(
x
,
y
)
q
j
(
y
)
d
s
=
∑
n
=
0
∞
‍
∑
m
=
-
n
n
‍
[
F
y
R
n
m
j
k
(
x
0
x
→
)
¯
L
n
,
j
m
(
1
)
(
x
0
)
vvvvvwiivv
+
H
y
R
n
m
j
(
x
0
x
→
)
¯
L
n
,
j
m
,
k
(
2
)
(
x
0
)
]
,
where
(28)
L
n
,
j
m
(
1
)
(
x
0
)
=
∑
n
′
=
0
∞
{
∑
m
′
=
0
n
′
‍
(
-
1
)
n
M
n
′
,
j
m
′
(
1
)
(
o
)
S
y
n
+
n
′
m
+
m
′
(
o
x
0
→
)
¯
‍
vvvvvivv
+
∑
m
′
=
1
k
‍
(
-
1
)
n
+
m
′
M
n
′
,
j
m
′
(
1
)
(
o
)
¯
S
y
n
+
n
′
m
-
m
′
(
o
x
0
→
)
¯
vvvvvivv
+
∑
m
′
=
k
+
1
n
′
‍
(
-
1
)
n
-
m
M
n
′
,
j
m
′
(
1
)
(
o
)
¯
S
y
n
+
n
′
m
′
-
m
(
o
x
0
→
)
}
,
L
n
,
j
m
,
k
(
2
)
(
x
0
)
=
-
(
o
x
0
→
)
k
L
n
,
j
m
(
1
)
(
x
0
)
+
∑
n
′
=
0
∞
‍
{
∑
m
′
=
0
n
′
‍
(
-
1
)
n
M
n
′
,
j
m
′
,
k
(
2
)
(
o
)
S
y
n
+
n
′
m
+
m
′
(
o
x
0
→
)
¯
vvvvvvivv
+
∑
m
′
=
1
k
‍
(
-
1
)
n
+
m
′
M
n
′
,
j
m
′
,
k
(
2
)
(
o
)
¯
S
y
n
+
n
′
m
-
m
′
(
o
x
0
→
)
¯
vvvvvvivv
+
∑
m
′
=
k
+
1
n
′
‍
(
-
1
)
n
-
m
M
n
′
,
j
m
′
,
k
(
2
)
(
o
)
¯
S
y
n
+
n
′
m
′
-
m
(
o
x
0
→
)
}
,
F
y
R
n
m
j
k
(
x
0
x
→
)
=
δ
k
j
R
y
n
m
(
x
0
x
→
)
-
(
x
0
x
→
)
∂
∂
x
j
(
R
y
n
m
(
x
0
x
→
)
)
,
H
y
R
n
m
j
(
x
0
x
→
)
=
∂
∂
x
j
(
R
y
n
m
(
x
0
x
→
)
)
.


### 5.2.4. Translation of a Local Expansion

Suppose thatv
(
x
) can be written as in (29), and suppose that x
1
∈
D
0; then for any point x
∈
D
0, v
(
x
) can be written according to the local moment centered at x
1 as follows:
(29)
v
k
(
x
)
=
∫
S
‍
G
k
j
(
x
,
y
)
q
j
(
y
)
d
s
=
∑
n
=
0
∞
‍
∑
m
=
-
n
n
‍
[
F
y
R
n
m
j
k
(
x
1
x
→
)
¯
L
n
,
j
m
(
1
)
(
x
1
)
vvvvvwiivv
+
H
y
R
n
m
j
(
x
1
x
→
)
¯
L
n
,
j
m
,
k
(
2
)
(
x
1
)
]
,
where
(30)
L
n
,
j
m
,
k
(
1
)
(
x
1
)
=
∑
n
′
=
0
∞
‍
{
∑
m
′
=
1
n
′
‍
(
-
1
)
m
L
n
′
,
j
m
′
(
1
)
(
x
0
)
R
y
n
-
n
′
m
+
m
′
(
x
0
x
1
→
)
¯
∑
m
′
=
0
k
‍
(
-
1
)
m
-
m
′
L
n
′
,
j
m
′
(
1
)
(
x
0
)
R
y
n
-
n
′
m
-
m
′
(
x
0
x
1
→
)
¯
vvvvvivv
+
∑
m
′
=
k
+
1
n
′
‍
L
n
′
,
j
m
′
(
1
)
(
x
0
)
R
y
n
-
n
′
m
′
-
m
(
x
0
x
1
→
)
vvvvvivv
+
∑
m
′
=
1
n
′
‍
(
-
1
)
m
L
n
′
,
j
m
′
(
1
)
(
x
0
)
¯
R
y
n
-
n
′
m
+
m
′
(
x
0
x
1
→
)
}
,
L
n
,
j
m
,
k
(
2
)
(
x
1
)
=
-
(
x
0
x
1
→
)
k
L
n
,
j
m
(
1
)
(
x
1
)
+
∑
n
′
=
0
∞
‍
{
∑
m
′
=
1
n
′
‍
(
-
1
)
m
L
n
′
,
j
m
′
,
k
(
2
)
(
x
0
)
¯
R
y
n
-
n
′
m
+
m
′
(
x
0
x
1
→
)
¯
∑
m
′
=
0
k
‍
(
-
1
)
m
-
m
′
L
n
′
,
j
m
′
,
k
(
2
)
(
x
0
)
R
y
n
-
n
′
m
-
m
′
(
x
0
x
1
→
)
¯
vvvvvvivv
+
∑
m
′
=
k
+
1
n
′
‍
L
n
′
,
j
m
′
,
k
(
2
)
(
x
0
)
R
y
n
-
n
′
m
′
-
m
(
x
0
x
1
→
)
vvvvvvivv
+
∑
m
′
=
1
n
′
‍
(
-
1
)
m
L
n
′
,
j
m
′
,
k
(
2
)
(
x
0
)
¯
R
y
n
-
n
′
m
+
m
′
(
x
0
x
1
→
)
¯
}
.


### 5.2.5. Error Bound

With a similar reasoning as in the error bound of theO
(
N
log
⁡
N
) scheme and under the same hypothesis of the representation (29), the evaluation of the velocity v
(
x
) at an order p generates an error bound expressed as follows:
(31)
E
p
′′
=
C
(
δ
j
k
ρ
-
r
)
+
(
1
ρ
-
r
)
+
(
r
ρ
2
-
ρ
r
)
(
r
ρ
)
p
+
1
,
where C
=
∫
C
‍
|
q
(
y
)
|
d
s.


## 5.2.1. Relations Formulas betweenS
y
n
m and R
y
n
m

Under the condition|
o
y
→
|
<
|
o
o
′
→
|, the function R
y
n
m satisfies a relationship of the following form:
(21)
R
y
n
m
(
o
′
y
→
)
=
∑
n
′
=
0
n
‍
{
∑
m
′
=
0
n
′
‍
(
-
1
)
m
′
R
y
n
′
m
′
(
o
′
o
→
)
¯
R
y
n
-
n
′
m
+
m
′
(
o
y
→
)
∑
m
′
=
0
k
‍
R
y
n
′
m
′
(
o
′
o
→
)
R
y
n
-
n
′
m
-
m
′
(
o
y
→
)
vvvvvivv
+
∑
m
′
=
k
+
1
n
′
‍
(
-
1
)
m
′
-
m
R
y
n
′
m
′
(
o
′
o
→
)
R
y
n
-
n
′
m
′
-
m
(
o
y
→
)
¯
vvvvivvv
+
∑
m
′
=
0
n
′
‍
(
-
1
)
m
′
R
y
n
′
m
′
(
o
′
o
→
)
¯
R
y
n
-
n
′
m
+
m
′
(
o
y
→
)
}
which allows a pole shift that will permit us to construct the multipolar to multipolar translation.Moreover, and under the condition|
x
0
x
→
|
<
|
o
x
0
→
|, the function S
y
n
m can be further expressed in the following form that will be used later to construct the conversion multipolar to local moment:
(22)
S
y
n
m
(
o
x
→
)
=
∑
n
′
=
0
∞
‍
{
∑
m
′
=
0
n
′
‍
(
-
1
)
n
′
R
y
n
′
m
′
(
x
0
x
→
)
S
y
n
+
n
′
m
+
m
′
(
o
x
0
→
)
vvvvvivv
+
∑
m
′
=
0
k
‍
(
-
1
)
n
′
+
m
′
R
y
n
′
m
′
(
x
0
x
→
)
¯
S
y
n
+
n
′
m
-
m
′
(
o
x
0
→
)
vvvvvivv
+
∑
m
′
=
k
+
1
n
′
‍
(
-
1
)
n
+
2
m
′
-
m
R
y
n
′
m
′
(
x
0
x
→
)
¯
S
y
n
+
n
′
m
′
-
m
(
o
x
0
→
)
}
.
Just like the first relation, under the condition |
x
0
x
→
|
<
|
x
0
x
1
→
|, we have a third equation that satisfies a pole relationship that will be used later to construct the local to local translation:
(23)
R
y
n
m
(
x
0
x
→
)
=
∑
n
′
=
n
∞
‍
{
∑
m
′
=
0
n
′
‍
(
-
1
)
m
′
R
y
n
′
m
′
(
x
1
x
→
)
¯
R
y
n
-
n
′
m
+
m
′
(
x
0
x
1
→
)
∑
m
′
=
0
k
‍
R
y
n
′
m
′
(
x
1
x
→
)
R
y
n
-
n
′
m
-
m
′
(
x
0
x
1
→
)
vvvvvivv
+
∑
m
′
=
k
+
1
n
′
‍
(
-
1
)
m
′
-
m
R
y
n
′
m
′
(
x
1
x
→
)
R
y
n
-
n
′
m
′
-
m
(
x
0
x
1
→
)
¯
vvvvvivv
+
∑
m
′
=
0
n
′
‍
(
-
1
)
m
′
R
y
n
′
m
′
(
x
1
x
→
)
¯
R
y
n
-
n
′
m
+
m
′
(
x
0
x
1
→
)
}
.


## 5.2.2. Translation of a Multipole Expansion

Let us suppose that the surfaceS contains uniform distribution points included inside a sphere D of radius a centered at a point o at a distance ρ from the origin o
′ with ρ
>
(
c
+
1
)
a and c
>
1. The integral equation is given by the multipole expansion (16):
(24)
v
k
(
x
)
=
∑
n
=
0
∞
‍
∑
m
=
-
n
n
‍
[
F
y
n
m
j
k
(
o
x
→
)
¯
M
n
,
j
m
(
1
)
(
o
)
vvvvvivv
+
H
y
n
m
j
(
o
x
→
)
¯
M
j
,
m
m
,
k
(
2
)
(
o
)
]
.
If the point x is outside the sphere D
0 of radius a
+
ρ, then the integral equation takes the following form:
(25)
v
k
(
x
)
=
∑
n
=
0
∞
‍
∑
m
=
-
n
n
[
F
y
n
m
j
k
(
o
′
x
→
)
¯
M
n
,
j
m
(
1
)
(
o
′
)
‍
vvvvvivv
+
H
y
n
m
j
(
o
′
x
→
)
¯
M
j
,
m
m
,
k
(
2
)
(
o
′
)
]
,
where
(26)
M
n
,
j
m
(
1
)
(
o
′
)
=
∑
n
′
=
0
n
‍
{
∑
m
′
=
1
n
′
‍
(
-
1
)
m
′
R
y
n
′
m
′
(
o
′
o
→
)
¯
M
n
-
n
′
,
j
m
+
m
′
(
1
)
(
o
)
∑
m
′
=
0
k
‍
R
y
n
′
m
′
(
o
′
o
→
)
M
n
-
n
′
,
j
m
-
m
′
(
1
)
(
o
)
vvvvvivv
+
∑
m
′
=
k
+
1
n
′
‍
(
-
1
)
m
′
-
m
R
y
n
′
m
′
(
o
′
o
→
)
M
n
-
n
′
,
j
m
′
-
m
(
1
)
(
o
)
¯
vvvvvivv
+
∑
m
′
=
1
n
′
‍
(
-
1
)
m
′
R
y
n
′
m
′
(
o
′
o
→
)
¯
M
n
-
n
′
,
j
m
+
m
′
(
1
)
(
o
)
}
,
M
n
,
j
m
,
k
(
2
)
(
o
′
)
=
(
o
′
o
→
)
k
M
n
,
j
m
(
1
)
(
o
′
)
+
∑
n
′
=
0
n
‍
{
∑
m
′
=
1
n
′
‍
(
-
1
)
m
′
R
y
n
′
m
′
(
o
′
o
→
)
¯
M
n
-
n
′
,
j
m
+
m
′
,
k
(
2
)
(
o
)
∑
m
′
=
0
k
‍
R
y
n
′
m
′
(
o
′
o
→
)
M
n
-
n
′
,
j
m
-
m
′
,
k
(
2
)
(
o
)
vvvvvvivv
+
∑
m
′
=
k
+
1
n
′
‍
(
-
1
)
m
′
-
m
R
y
n
′
m
′
(
o
′
o
→
)
M
n
-
n
′
,
j
m
′
-
m
,
k
(
2
)
(
o
)
¯
vvvvvvivv
+
∑
m
′
=
1
n
′
‍
(
-
1
)
m
′
R
y
n
′
m
′
(
o
′
o
→
)
¯
M
n
-
n
′
,
j
m
+
m
′
,
k
(
2
)
(
o
)
}
.


## 5.2.3. Conversion of a Multipole Expansion into a Local Expansion

Let a surfaceS containing uniform distribution points included inside a sphere D of radius a centered at a point o at a distance ρ from a given point x
0 with ρ
>
(
c
+
1
)
a and c
>
1. The corresponding multipole expansion (16) converges for any points x inside the sphere D
0 of radius a centered at x
0. The velocity is given by
(27)
v
k
(
x
)
=
∫
S
‍
G
k
j
(
x
,
y
)
q
j
(
y
)
d
s
=
∑
n
=
0
∞
‍
∑
m
=
-
n
n
‍
[
F
y
R
n
m
j
k
(
x
0
x
→
)
¯
L
n
,
j
m
(
1
)
(
x
0
)
vvvvvwiivv
+
H
y
R
n
m
j
(
x
0
x
→
)
¯
L
n
,
j
m
,
k
(
2
)
(
x
0
)
]
,
where
(28)
L
n
,
j
m
(
1
)
(
x
0
)
=
∑
n
′
=
0
∞
{
∑
m
′
=
0
n
′
‍
(
-
1
)
n
M
n
′
,
j
m
′
(
1
)
(
o
)
S
y
n
+
n
′
m
+
m
′
(
o
x
0
→
)
¯
‍
vvvvvivv
+
∑
m
′
=
1
k
‍
(
-
1
)
n
+
m
′
M
n
′
,
j
m
′
(
1
)
(
o
)
¯
S
y
n
+
n
′
m
-
m
′
(
o
x
0
→
)
¯
vvvvvivv
+
∑
m
′
=
k
+
1
n
′
‍
(
-
1
)
n
-
m
M
n
′
,
j
m
′
(
1
)
(
o
)
¯
S
y
n
+
n
′
m
′
-
m
(
o
x
0
→
)
}
,
L
n
,
j
m
,
k
(
2
)
(
x
0
)
=
-
(
o
x
0
→
)
k
L
n
,
j
m
(
1
)
(
x
0
)
+
∑
n
′
=
0
∞
‍
{
∑
m
′
=
0
n
′
‍
(
-
1
)
n
M
n
′
,
j
m
′
,
k
(
2
)
(
o
)
S
y
n
+
n
′
m
+
m
′
(
o
x
0
→
)
¯
vvvvvvivv
+
∑
m
′
=
1
k
‍
(
-
1
)
n
+
m
′
M
n
′
,
j
m
′
,
k
(
2
)
(
o
)
¯
S
y
n
+
n
′
m
-
m
′
(
o
x
0
→
)
¯
vvvvvvivv
+
∑
m
′
=
k
+
1
n
′
‍
(
-
1
)
n
-
m
M
n
′
,
j
m
′
,
k
(
2
)
(
o
)
¯
S
y
n
+
n
′
m
′
-
m
(
o
x
0
→
)
}
,
F
y
R
n
m
j
k
(
x
0
x
→
)
=
δ
k
j
R
y
n
m
(
x
0
x
→
)
-
(
x
0
x
→
)
∂
∂
x
j
(
R
y
n
m
(
x
0
x
→
)
)
,
H
y
R
n
m
j
(
x
0
x
→
)
=
∂
∂
x
j
(
R
y
n
m
(
x
0
x
→
)
)
.


## 5.2.4. Translation of a Local Expansion

Suppose thatv
(
x
) can be written as in (29), and suppose that x
1
∈
D
0; then for any point x
∈
D
0, v
(
x
) can be written according to the local moment centered at x
1 as follows:
(29)
v
k
(
x
)
=
∫
S
‍
G
k
j
(
x
,
y
)
q
j
(
y
)
d
s
=
∑
n
=
0
∞
‍
∑
m
=
-
n
n
‍
[
F
y
R
n
m
j
k
(
x
1
x
→
)
¯
L
n
,
j
m
(
1
)
(
x
1
)
vvvvvwiivv
+
H
y
R
n
m
j
(
x
1
x
→
)
¯
L
n
,
j
m
,
k
(
2
)
(
x
1
)
]
,
where
(30)
L
n
,
j
m
,
k
(
1
)
(
x
1
)
=
∑
n
′
=
0
∞
‍
{
∑
m
′
=
1
n
′
‍
(
-
1
)
m
L
n
′
,
j
m
′
(
1
)
(
x
0
)
R
y
n
-
n
′
m
+
m
′
(
x
0
x
1
→
)
¯
∑
m
′
=
0
k
‍
(
-
1
)
m
-
m
′
L
n
′
,
j
m
′
(
1
)
(
x
0
)
R
y
n
-
n
′
m
-
m
′
(
x
0
x
1
→
)
¯
vvvvvivv
+
∑
m
′
=
k
+
1
n
′
‍
L
n
′
,
j
m
′
(
1
)
(
x
0
)
R
y
n
-
n
′
m
′
-
m
(
x
0
x
1
→
)
vvvvvivv
+
∑
m
′
=
1
n
′
‍
(
-
1
)
m
L
n
′
,
j
m
′
(
1
)
(
x
0
)
¯
R
y
n
-
n
′
m
+
m
′
(
x
0
x
1
→
)
}
,
L
n
,
j
m
,
k
(
2
)
(
x
1
)
=
-
(
x
0
x
1
→
)
k
L
n
,
j
m
(
1
)
(
x
1
)
+
∑
n
′
=
0
∞
‍
{
∑
m
′
=
1
n
′
‍
(
-
1
)
m
L
n
′
,
j
m
′
,
k
(
2
)
(
x
0
)
¯
R
y
n
-
n
′
m
+
m
′
(
x
0
x
1
→
)
¯
∑
m
′
=
0
k
‍
(
-
1
)
m
-
m
′
L
n
′
,
j
m
′
,
k
(
2
)
(
x
0
)
R
y
n
-
n
′
m
-
m
′
(
x
0
x
1
→
)
¯
vvvvvvivv
+
∑
m
′
=
k
+
1
n
′
‍
L
n
′
,
j
m
′
,
k
(
2
)
(
x
0
)
R
y
n
-
n
′
m
′
-
m
(
x
0
x
1
→
)
vvvvvvivv
+
∑
m
′
=
1
n
′
‍
(
-
1
)
m
L
n
′
,
j
m
′
,
k
(
2
)
(
x
0
)
¯
R
y
n
-
n
′
m
+
m
′
(
x
0
x
1
→
)
¯
}
.


## 5.2.5. Error Bound

With a similar reasoning as in the error bound of theO
(
N
log
⁡
N
) scheme and under the same hypothesis of the representation (29), the evaluation of the velocity v
(
x
) at an order p generates an error bound expressed as follows:
(31)
E
p
′′
=
C
(
δ
j
k
ρ
-
r
)
+
(
1
ρ
-
r
)
+
(
r
ρ
2
-
ρ
r
)
(
r
ρ
)
p
+
1
,
where C
=
∫
C
‍
|
q
(
y
)
|
d
s.


## 6. Numerical Results

We present here some numerical results given by applying the FMM-BEM to the Stokes problem. The method uses several parameters on which depend the speed and the precision of the calculation, that is, the size of the cells, the level number of the grids, the order of truncation of the multipole series expansion, and so forth. The optimal values to obtain a good equilibrium between effectiveness and precision are given for the mono- and multilevel approaches. Theoretical complexities are checked numerically. Problems of sizeN
=
O
(
1
0
6
) are solved in about two hours by iterations on a PC Pentium 4, 3.2GHZ with 4Go-RAM.


### 6.1. Configuration Test

#### 6.1.1. Sphere Triangulation

The starting point is to generate a grid of many spheres that we put in a cube that represents the initial domain at the initial level of the tree (cf. Figure3). To generate the collocation points on the spheres, we used a simple technique, presented by Pozrikidis [16] for the discretization of a sphere. Beginning from an octahedron, we developed a triangulation by subdividing recursively each triangle into four subtriangles until obtaining the desired level of discretisation. At each step, the triangle is projected on the surface of the sphere (cf. Figure 4).Initial configuration.
(a)
(b)Figure 4
Sphere triangulation.


#### 6.1.2. Numerical Parameters

For the numerical consideration, we considered that two boxes are far away if they are separated by a third box. In addition, we tookr
>
3
a for the O
(
N
log
⁡
N
) scheme and ρ
>
3
r for the O
(
N
) scheme. The error given by (20) is thus controlled by (
3
C
/
2
L
)
2
l
(
1
/
3
)
p
+
1, where l is the level of the tree on which one approximates the far field interactions, and L represents the dispatcher of the initial cube. In the results of the tests presented below, we took p
=
12 and we worked on a tree of depth equal to 5 which give a plug of error equal to O
(
1
0
6
). We obtained the same thing for the error given by (31).


### 6.2. Presentation and Interpretation of the Results

It is worth to be pointed out that the maximum value of the error obtained in the different tests is about1
0
-
6. These tests permit us to validate the method and show its precision. We remark from the calculation that the error introduced by the FMM compared to the traditional BEM does not have a practical incidence on the quality of the result.


#### 6.2.1.O
(
N
log
N
)

We present here the results obtained with theO
(
N
log
⁡
N
) scheme. Figure 5 that represents the CPU time versus the number N
c of collocation point at different levels (
N
c
=
2
(
N
t
+
1
)
) shows that the method FMM in the case of the scheme O
(
N
log
⁡
N
) is useful not only to minimize the memory storage of the matrix due to the compression carrying out but also to enable us to gain in terms of computing time. This method appears to be much better than the direct calculation for the large values of N where N
=
3
N
p
(
N
2
c
1
). Indeed we can see that the effectiveness of this method appears from a certain number N
c evaluated by the results between 1
0
4 and 3
×
1
0
4. For N
c
≤
9202 the method FMM does not have any influence from the computing times point of view and from N
c
=
27702 the efficiency of the FMM method becomes to be visible; see Figures 5 and 8. To be more precise, we consider two boxes, shown, respectively, in Figures 6 and 7, that are obtained from Figure 5 by zooming on the regions that seem to be interesting for the interpretation. For N
c
=
27702 the efficiency of the method starts from level 2 (see Figure 6) and then for N
c
=
351918 the efficiency becomes to be visible from level 3 (see Figure 7). We can also see from Figure 7 that level 4 becomes better than level 2 from a certain number N evaluated to be between 2.75
×
1
0
5 and 3
×
1
0
5. Thus, as a general rule about the potential of the FMM method, we can conclude that the more the size of the problem N increases, the more we have to work down in the hierarchical tree of the code.Figure 5
CPU time versus the numberN
c of collocation points at different levels.Figure 6
Zoom 1 of Figure5.Figure 7
Zoom 2 of Figure5.CPU time versus the level number for different collocation pointsN
c.
(a)
(b)


#### 6.2.2.O
(
N
)

The interpretation of the results found for theO
(
N
) scheme is similar to that obtained for the O
(
N
log
⁡
N
) one. The same remarks given from all figures for the O
(
N
log
⁡
(
N
)
) scheme may be reproduced here from Figures 9, 10, 11, and 12 which correspond to the O
(
N
). The comparison of the results of the two schemes O
(
N
log
⁡
N
) and O
(
N
) is illustrated in Figures 13, 14, 15, and 16. It should be noted that the O
(
N
) scheme becomes more efficient than the O
(
N
log
⁡
N
) scheme after a certain number N; see Figure 16. However, one has to work under a certain depth of the tree; otherwise, the O
(
N
log
⁡
N
) becomes more efficient from the CPU time point of view; see Figure 15.Figure 9
CPU time versus particles number for different levels.Figure 10
Zoom 1 of Figure9.Figure 11
Zoom 2 of Figure9.CPU time versus the level for several values ofN.
(a)
(b)Figure 13
CPU time of the two schemes versus the level forN
c
=
351918.Figure 14
CPU time of the direct calculus and the two schemesO
(
N
log
⁡
N
) and O
(
N
) versus the collocation points N
c for level = 3.Figure 15
CPU time of the direct calculus and the two schemesO
(
N
log
⁡
N
) and O
(
N
) versus the level for several values of N
c.Figure 16
CPU time of the direct calculus and the two schemesO
(
N
log
⁡
N
) and O
(
N
) versus the level for several values of N
c.


## 6.1. Configuration Test

### 6.1.1. Sphere Triangulation

The starting point is to generate a grid of many spheres that we put in a cube that represents the initial domain at the initial level of the tree (cf. Figure3). To generate the collocation points on the spheres, we used a simple technique, presented by Pozrikidis [16] for the discretization of a sphere. Beginning from an octahedron, we developed a triangulation by subdividing recursively each triangle into four subtriangles until obtaining the desired level of discretisation. At each step, the triangle is projected on the surface of the sphere (cf. Figure 4).Initial configuration.
(a)
(b)Figure 4
Sphere triangulation.


### 6.1.2. Numerical Parameters

For the numerical consideration, we considered that two boxes are far away if they are separated by a third box. In addition, we tookr
>
3
a for the O
(
N
log
⁡
N
) scheme and ρ
>
3
r for the O
(
N
) scheme. The error given by (20) is thus controlled by (
3
C
/
2
L
)
2
l
(
1
/
3
)
p
+
1, where l is the level of the tree on which one approximates the far field interactions, and L represents the dispatcher of the initial cube. In the results of the tests presented below, we took p
=
12 and we worked on a tree of depth equal to 5 which give a plug of error equal to O
(
1
0
6
). We obtained the same thing for the error given by (31).


## 6.1.1. Sphere Triangulation

The starting point is to generate a grid of many spheres that we put in a cube that represents the initial domain at the initial level of the tree (cf. Figure3). To generate the collocation points on the spheres, we used a simple technique, presented by Pozrikidis [16] for the discretization of a sphere. Beginning from an octahedron, we developed a triangulation by subdividing recursively each triangle into four subtriangles until obtaining the desired level of discretisation. At each step, the triangle is projected on the surface of the sphere (cf. Figure 4).Initial configuration.
(a)
(b)Figure 4
Sphere triangulation.


## 6.1.2. Numerical Parameters

For the numerical consideration, we considered that two boxes are far away if they are separated by a third box. In addition, we tookr
>
3
a for the O
(
N
log
⁡
N
) scheme and ρ
>
3
r for the O
(
N
) scheme. The error given by (20) is thus controlled by (
3
C
/
2
L
)
2
l
(
1
/
3
)
p
+
1, where l is the level of the tree on which one approximates the far field interactions, and L represents the dispatcher of the initial cube. In the results of the tests presented below, we took p
=
12 and we worked on a tree of depth equal to 5 which give a plug of error equal to O
(
1
0
6
). We obtained the same thing for the error given by (31).


## 6.2. Presentation and Interpretation of the Results

It is worth to be pointed out that the maximum value of the error obtained in the different tests is about1
0
-
6. These tests permit us to validate the method and show its precision. We remark from the calculation that the error introduced by the FMM compared to the traditional BEM does not have a practical incidence on the quality of the result.


### 6.2.1.O
(
N
log
N
)

We present here the results obtained with theO
(
N
log
⁡
N
) scheme. Figure 5 that represents the CPU time versus the number N
c of collocation point at different levels (
N
c
=
2
(
N
t
+
1
)
) shows that the method FMM in the case of the scheme O
(
N
log
⁡
N
) is useful not only to minimize the memory storage of the matrix due to the compression carrying out but also to enable us to gain in terms of computing time. This method appears to be much better than the direct calculation for the large values of N where N
=
3
N
p
(
N
2
c
1
). Indeed we can see that the effectiveness of this method appears from a certain number N
c evaluated by the results between 1
0
4 and 3
×
1
0
4. For N
c
≤
9202 the method FMM does not have any influence from the computing times point of view and from N
c
=
27702 the efficiency of the FMM method becomes to be visible; see Figures 5 and 8. To be more precise, we consider two boxes, shown, respectively, in Figures 6 and 7, that are obtained from Figure 5 by zooming on the regions that seem to be interesting for the interpretation. For N
c
=
27702 the efficiency of the method starts from level 2 (see Figure 6) and then for N
c
=
351918 the efficiency becomes to be visible from level 3 (see Figure 7). We can also see from Figure 7 that level 4 becomes better than level 2 from a certain number N evaluated to be between 2.75
×
1
0
5 and 3
×
1
0
5. Thus, as a general rule about the potential of the FMM method, we can conclude that the more the size of the problem N increases, the more we have to work down in the hierarchical tree of the code.Figure 5
CPU time versus the numberN
c of collocation points at different levels.Figure 6
Zoom 1 of Figure5.Figure 7
Zoom 2 of Figure5.CPU time versus the level number for different collocation pointsN
c.
(a)
(b)


### 6.2.2.O
(
N
)

The interpretation of the results found for theO
(
N
) scheme is similar to that obtained for the O
(
N
log
⁡
N
) one. The same remarks given from all figures for the O
(
N
log
⁡
(
N
)
) scheme may be reproduced here from Figures 9, 10, 11, and 12 which correspond to the O
(
N
). The comparison of the results of the two schemes O
(
N
log
⁡
N
) and O
(
N
) is illustrated in Figures 13, 14, 15, and 16. It should be noted that the O
(
N
) scheme becomes more efficient than the O
(
N
log
⁡
N
) scheme after a certain number N; see Figure 16. However, one has to work under a certain depth of the tree; otherwise, the O
(
N
log
⁡
N
) becomes more efficient from the CPU time point of view; see Figure 15.Figure 9
CPU time versus particles number for different levels.Figure 10
Zoom 1 of Figure9.Figure 11
Zoom 2 of Figure9.CPU time versus the level for several values ofN.
(a)
(b)Figure 13
CPU time of the two schemes versus the level forN
c
=
351918.Figure 14
CPU time of the direct calculus and the two schemesO
(
N
log
⁡
N
) and O
(
N
) versus the collocation points N
c for level = 3.Figure 15
CPU time of the direct calculus and the two schemesO
(
N
log
⁡
N
) and O
(
N
) versus the level for several values of N
c.Figure 16
CPU time of the direct calculus and the two schemesO
(
N
log
⁡
N
) and O
(
N
) versus the level for several values of N
c.


## 6.2.1.O
(
N
log
N
)

We present here the results obtained with theO
(
N
log
⁡
N
) scheme. Figure 5 that represents the CPU time versus the number N
c of collocation point at different levels (
N
c
=
2
(
N
t
+
1
)
) shows that the method FMM in the case of the scheme O
(
N
log
⁡
N
) is useful not only to minimize the memory storage of the matrix due to the compression carrying out but also to enable us to gain in terms of computing time. This method appears to be much better than the direct calculation for the large values of N where N
=
3
N
p
(
N
2
c
1
). Indeed we can see that the effectiveness of this method appears from a certain number N
c evaluated by the results between 1
0
4 and 3
×
1
0
4. For N
c
≤
9202 the method FMM does not have any influence from the computing times point of view and from N
c
=
27702 the efficiency of the FMM method becomes to be visible; see Figures 5 and 8. To be more precise, we consider two boxes, shown, respectively, in Figures 6 and 7, that are obtained from Figure 5 by zooming on the regions that seem to be interesting for the interpretation. For N
c
=
27702 the efficiency of the method starts from level 2 (see Figure 6) and then for N
c
=
351918 the efficiency becomes to be visible from level 3 (see Figure 7). We can also see from Figure 7 that level 4 becomes better than level 2 from a certain number N evaluated to be between 2.75
×
1
0
5 and 3
×
1
0
5. Thus, as a general rule about the potential of the FMM method, we can conclude that the more the size of the problem N increases, the more we have to work down in the hierarchical tree of the code.Figure 5
CPU time versus the numberN
c of collocation points at different levels.Figure 6
Zoom 1 of Figure5.Figure 7
Zoom 2 of Figure5.CPU time versus the level number for different collocation pointsN
c.
(a)
(b)


## 6.2.2.O
(
N
)

The interpretation of the results found for theO
(
N
) scheme is similar to that obtained for the O
(
N
log
⁡
N
) one. The same remarks given from all figures for the O
(
N
log
⁡
(
N
)
) scheme may be reproduced here from Figures 9, 10, 11, and 12 which correspond to the O
(
N
). The comparison of the results of the two schemes O
(
N
log
⁡
N
) and O
(
N
) is illustrated in Figures 13, 14, 15, and 16. It should be noted that the O
(
N
) scheme becomes more efficient than the O
(
N
log
⁡
N
) scheme after a certain number N; see Figure 16. However, one has to work under a certain depth of the tree; otherwise, the O
(
N
log
⁡
N
) becomes more efficient from the CPU time point of view; see Figure 15.Figure 9
CPU time versus particles number for different levels.Figure 10
Zoom 1 of Figure9.Figure 11
Zoom 2 of Figure9.CPU time versus the level for several values ofN.
(a)
(b)Figure 13
CPU time of the two schemes versus the level forN
c
=
351918.Figure 14
CPU time of the direct calculus and the two schemesO
(
N
log
⁡
N
) and O
(
N
) versus the collocation points N
c for level = 3.Figure 15
CPU time of the direct calculus and the two schemesO
(
N
log
⁡
N
) and O
(
N
) versus the level for several values of N
c.Figure 16
CPU time of the direct calculus and the two schemesO
(
N
log
⁡
N
) and O
(
N
) versus the level for several values of N
c.


## 7. Concluding Remarks

We presented in this work the fast multipole method applied to the Stokes problem in a 3D infinite viscous fluid in the presence of several rigid bodies. The formulas obtained for both schemesO
(
N
log
⁡
N
) and O
(
N
), where N is the size of the matrix obtained from the discretization of the problem, are implemented numerically. The numerical results obtained are very satisfactory and incite us to use the method in the future in a traditional IBEM code that will allow for an efficient resolution of the linear system by iterative methods such as GMRES. This work makes it possible to investigate in future works more complex and realistic configurations such as the sedimentation of dilute suspensions composed by a great number of particles or the movements of deformable particles. Such configurations could enable us to compare between the numerical results and the experimental ones. For the case of deformable bodies, the integral equation must take into account the double potential layer. The development of an adequate FMM-BEM based on formulae given previously constitutes an interesting perspective for a future area of research.



---

*Source: 563603-2014-06-03.xml*