Fast algorithms applied to the acoustical energy boundary element method
Matthias RAM(1) , Otto VON ESTORFF(1)(1) Hamburg University of Technology, Institute of Modelling and Computation, Germany, mub@tuhh.de

Abstract
The Energy Boundary Element Method (EBEM) aims at the solution of acoustic high frequency problems,
where the classic BEM becomes inefficient due to the large number of degrees of freedom (DOF) required. A
transition to non-phased energetic state variables removes the correlation between the investigated frequency
range and the element size. In the EBEM the number of DOF is thus determined primarily by the complexity of
the geometry. In order to speed up the computation and to handle structures of very high geometric complexity,
a fast multipole algorithm for the EBEM is investigated. It has a great potential to reduce the numerical effort
in the Boundary Element Method (BEM). Helmholtz and Laplace problems have been solved very efficiently
applying this algorithm. However, the kernels used in the EBEM require an adaptation of the algorithm. In thispaper the development of a fast multipole formulation of the EBEM is presented.
Keywords: Fast multipole, Boundary elements, Energy method

1

INTRODUCTION

Nowadays, the BEM can be considered a well-established solution technique for acoustic problems of manykinds. The number of a models DOF N required for satisfactory solution accuracy, however, does not only
depend on the complexity of the geometry, but also on the investigated frequency range. A commonly used
rule of thumb suggests at least 6-10 elements per wave length. Since the matrices of the system of equations in
the BEM are generally fully populated and non-symmetric, the memory required for the storage and solution, as
well as the solution time, grow very fast with the maximum frequency under consideration. Thus, the BEM ismostly limited to the low and medium frequency ranges. The numerical effort for the conventional BEM grows
proportionally with O N 2 . The application of the Fast Multipole Method (FMM) allows the reduction of thecomplexity to O (N) respectively O (N log N) [1]. Although this improvement means a tremendous accelerationand considerably extends the frequency range where BEM is applicable, the FMM accelerated BEM (FMBEM)
still fails to cover the entire audible frequency range in most practical problems.
The EBEM constitutes a different approach to deal with the particular challenges of high frequency problems[2, 3, 4, 5]. By neglecting interferences between propagative waves and using energetic quantities for thedescription of the sound field, similarly to the statistical energy analysis (SEA), the correlation between theinvestigated frequency and the element size is removed. To satisfy the underlying assumptions, the calculations
have to be carried out within a broader frequency range that contains several eigenfrequencies. The method isthus only applicable for higher frequencies, where the modal overlap is high and yields a frequency-smoothed
response. Unlike in the SEA, the results obtained by the EBEM give detailed insight in the local spatial
distribution of the variables describing the sound field.
This paper aims at the formulation of a FMM for the EBEM (FMEBEM) to speed up the computation andallow for finer discretiziation, where it is required by high geometric complexity of the model. The approach,first presented in [6], is discussed in more detail. Section 2 gives an outline of general FMM procedures. In
section 3 the application to the EBEM is shown followed by numerical experiments to verify the accuracy andperformance of the FMEBEM presented in section 4. The paper closes with a summary and an outlook in
section 5.

5457

2

FUNDAMENTALS OF FMM ACCELERATED BOUNDARY ELEMENT METHODS

The following section is based on references [1, 7, 8] which may be consulted for more in depth information.
The boundary integral equations occurring in the BEM are of the formZf (x) =x ∈ Γ,

K(x, y)φ (y) dy(1)Γwhere x and y denote the position of receiver and source point, respectively. The function f is known on theboundary Γ of the computational domain, φ is an unknown function and K (x, y) is a given kernel function. The
underlying idea of the FMBEM is to avoid the costly evaluation of all the direct interactions between separate
elements when evaluating the kernel function. Instead, several neighboring source elements are clustered andtheir influence on sufficiently far situated receiver elements is evaluated by some suitable approximation resulting
in a "post office scenario" as illustrated in Figure 1.

(a)

(b)

Figure 1. Direct evaluation of interactions (a) and evaluation in a single stage FMM (b) according to [7]2.1 Definitions
The mathematical basis for the FMM is given by a factorization of the kernel function of the form
∞

(1)

(2)

K(x, y) = ∑ kn (x − y0 )kn (y − y0 ).

(2)

n=0

(1)

In Equation 2, y0 denotes an arbitrary point in the R3 . The functions kn are usually singular in the origin,(2)
whereas kn are usually entire functions [1]. The separation of receiver and source in Equation 2 principally
allows the desired clustering after truncation of the infinite series. The expansion, however, is usually valid
under a certain condition only, the most common beingkx − y0 k > ky − y0 k .

(3)To ensure this condition, a hierarchical tree structure is introduced. In 3D, a so called oct-tree is used. To
obtain the tree, a cubic bounding box is defined, so that it contains the entire boundary surface. This boundingbox is said to be on level 0 and contains all the elements of the model. Boxes at higher levels l > 0 are createdby splitting the parent box on level l − 1 into eight equally sized children. The division is repeated until either
a desired level or a limiting number of elements per box is reached. A childless box on the deepest level iscalled leaf. Furthermore the following definitions are introduced: Two boxes are said to be near neighbors ifthey are on the same level and share at least one boundary point. A box is within the interaction list of anotherbox of the same level if their parents are near neighbors but the boxes themselves are not.
Boxes that are not near neighbors are well separated, meaning that Equation 3 is satisfied for all points of thereceiving box when y0 is chosen as the center point of the box containing the source point. Thus, it becomes
obvious that a tree structure in the FMM must always consist of at least three levels, since on level 0 and
1 there are no well separated boxes. At level 2 on the other hand, all well separated boxes are within theinteraction list of a box and their influence can be taken into account by evaluation of Equation 2. If the depthof the tree exceeds level 2, the method at hand is called Multi Level Fast Multipole Method (MLFMM).

5458

Using Equation 1 and Equation 2 a multipole expansion for Γ0 , the part of Γ contained within the box, and apoint x outside the box is obtained as
∞

Z
Γ0

(1)

K(x, y)φ (y) dy = ∑ kn (x − y0 )Mn (y0 ),

(4)

n=0

where the multipole moment centered at y0 is defined byZ

Mn (y0 ) =(2)

kn (y − y0 )φ (y) dy.

(5)

Γ0

In a similar manner, a local expansion is defined for a point x inside a box as
∞

Z
Γ0

(3)

K(x, y)φ (y) dy = ∑ Ln (x0 )kn (x − x0 ),

(6)

n=0

(3)

where Ln (x0 ) denotes the local coefficients centered at x0 and the functions kn are chosen to be a complete
set of interior solutions of the governing partial differential equation. For obvious reasons, the summations inEquation 4 and 6 have to be truncated after a finite number of terms in practical implementations.
To complete the tool set for the MLFMM, three operators are required: An operator to shift the center ofa multipole expansion from a box to its parent, called M2M, an operator to translate a multipole expansion
centered at a box to a local expansion centered at another box in its interaction list, which is called M2L andan operator to shift the center of a local expansion from a box to its child, called L2L.
2.2 Algorithm
The FMM can be employed for the approximation of a matrix-vector product and as such in combination withan iterative solver for the solution of a boundary value problem by means of the BEM. The major steps ofthe procedure are as follows. First, the boundary is discretized and a trial vector is determined as usual. The
elements are organized in a tree structured as described in subsection 2.1. Starting at the deepest level, themultipole moments for all leaf boxes are calculated using Equation 5. The multipole moments are then shiftedto the respective parent boxes using M2M, summing up the contributions of all the children of a parent. This
procedure is repeated following the hierarchy of the tree upwards until all boxes on level 2 posses a multipole
moment. Once this step is completed, the multipole moments on each level are translated to local coefficientsof the interaction list using M2L and summed up for each box. Now, starting from level 2, tracing the tree
downwards, the local coefficients of each parent box are shifted to its children using L2L and are added totheir existing coefficients. Once again, this step is repeated until all leaf boxes posses a local expansion that
accounts for the influence of all well separated boxes. This expansion is then evaluated using Equation 6. The
contributions of all source points within the same box and near neighbor boxes are evaluated directly and are
added to the result from the local expansion.

3

APPLICATION OF THE FMM TO THE EBEM

The boundary integral equation used in the EBEM is [3, 5]Z

H(x, y) · nx σ (y) dy + c(x)σ (x) = In (x).

(7)Γ

In Equation 7, nx denotes the outward normal vector at the receiving point, σ the source strength on theboundary, c(x) the boundary factors and In (x) the sound intensity in normal direction. The kernel function is
given byH(x, y) =e−µr r
,
4π r3

(8)

5459

where r = x − y, r = kx − yk and µ is the attenuation coefficient. The strong singularity occurring in Equation 8
may be handled by a suitable numerical integration scheme [9]. In case of constant elements these singularities
are avoided. To apply the FMM to Equation 7 an expansion in the form of Equation 2 for Equation 8 isrequired. In its general form, no obvious factorization of the kernel function seems to be available. Since thekernel is not oscillatory and smooth outside the singularity, the application of a so called kernel-independent or
black box FMM [11] would be feasible, where the factorization does not need to be explicitly known. For the
special case µ = 0, however, which is a valid assumption for many practical problems, an analytical formulation
as described in section 2 can be found. This analytical FMM for the undamped case will be described in thefollowing. As µ = 0, the following identity is observed
H(x, y) · nx =∂ G(x, y)r · nx
≡−
.
4πr3
∂ nx

(9)

1
denotes the Green’s function of Laplace’s equation. Thus, the scalar product H(x, y) · nx iswhere G(x, y) = 4πrtaken as the negative normal derivative of the Laplace kernel function at the receiver point. With this approach,a FMM based on the factorization of G(x, y), known as [10]G(x, y) =1 ∞ n m∑ ∑ S̄n (x − y0 )Rmn (y − y0 ),4π n=0
m=−n

(10)

m
is possible. In Equation 10 the symbol ¯· indicates the complex conjugate. The functions Rm
n and Sn are known
as regular and singular solid harmonics. They are defined in terms of the polar coordinates r, φ , θ of a vector x
as1
Pm (cos θ )eimφ rn(n + m)! n

(11)

1
1
Pm (cos θ )eimφ n+1 ,(n − m)! n
r

(12)

Rm
n (x) =andSnm (x) =where Pnm is the associated Legendre function of degree n and order m. Here, the definition of Pnm from [10] is
used, where the factor (−1)m , known as Condon-Shortley phase, is omitted from the definition. By comparing
m(l)
m(l)
(2)
(1)
Equation 10 to Equation 2 it becomes obvious that kl (x − y0 ) = S̄n(l) (x − y0 ) and kl (y − y0 ) = Rn(l) (y − y0 )have been chosen, where the original summation index n has been replaced with l to avoid confusion withassociated Legendre functions degree n.
The normal derivative appearing in Equation 9 depends on the receiver position only. Due to the separation ofsource and receiver, its computation requires a modification of the last step of the FMM algorithm only, whichis the evaluation of at the receiver position by Equation 6. Using
m(l)

(3)

kl (x − x0 ) =∂ Rn(l) (x − x0 )

,

(13)

1 ∞ n
∂ Rm (x − x0 )Lnm (x0 ) n.
∑∑
4π n=0 m=−n
∂ nx

(14)

∂ nx

Equation 6 can be rewritten asZ

H(x, y) · nx σ (y) dy = −Γ0Taking the discretization into account, the multipole moment of the leaf boxes may be calculated from theelements i inside the box with their shape functions Ni as
Mnm (y0 ) = ∑Zi∈Box Γi5460

Rmn (y − y0 )Ni (y)σ (y) dy.

(15)

The operators M2M, M2L and L2L can be adopted from the FMM for the Laplace equation. These are [10]Mnm (y00 ) =

n

n0

0

0

m−m
0 (y0 ),
∑ 0 Rmn0 (y0 − y00 )Mn−n
∑
0
0

(16)

n =0 m =−n∞

Lnm (x0 ) = ∑

n0

0

0

m
m+m
0 (x0 − y0 )Mn0 (y0 ),∑ 0 (−1)n S̄n+n0(17)n0 =0 m =−n

Lnm (x1 ) =

∞

n0

0

0

−m
(x1 − x0 )Lnm0 (x0 ).
∑ 0 Rmn0 −n
∑
0
0

(18)

n =n m =−nWhere Equation 16 represents M2M, M2L is given by Equation 17 and Equation 18 defines L2L. For theevaluation of Equation 14 the derivatives of Rmn need to be calculated. They are given by [10]∂ m 1 m−1
Rn = Rn−1 − Rm+1
n−1 ,
∂x
2

∂ m i m−1
R = R
+ Rm+1n−1 ,∂ y n 2 n−1
∂ m
R =Rm
n−1 .
∂z n

(19)
(20)
(21)

As has been mentioned before, the expansion is truncated, thus the outer sum is terminated after p terms,
resulting in a total p2 coefficients for the local and multipole expansion of each box. The number of terms p
is refered to as expansion length.

4

RESULTS

The previously described FMEBEM was implemented within the framework of a C++ library for acoustic
boundary element analyses. The method was implemented as a MLFMM, where the number of levels nl ischosen according to [12] bynl ≈ log8 (N) .

(22)

4.1 Sphere
For a first verification and to perform some numerical experiments concerning the performance and accuracy ofthe method, the geometry of a sphere is analyzed. Boundary conditions for the normal intensity are chosen suchthat they correspond to a point source in the center of the sphere. The sphere is modeled using 9,359 elements
and constant shape functions. In Figure 2, the resulting source strengths on the boundary are shown as computedby the direct evaluation of the boundary integral equation (Figure 2a) and the FMEBEM (Figure 2b-2e).

(a) direct EBEM(b) FMEBEM, p = 3(c) FMEBEM, p = 4(d) FMEBEM, p = 5(e) FMEBEM, p = 6

Figure 2. Source strength on the boundary computed by direct evaluation of Equation 7 (a) and FMM withvarying expansions lengths (b-e)

5461

Using the expansion length p = 3, there are clearly visible artifacts in the solution, showing deviations from thedirect solution within the order of magnitude of roughly 10 %. With increasing expansion length, however, these
deviations decline and for p = 6 the results of the FMEBEM match those of direct computation
very well. The
operators M2M, M2L and L2L that have been implemented are of the complexity O p4 , hence the expansion
length can have great impact on the computation times and thus the efficiency of the FMM. To analyze theeffect of the expansion length on accuracy and speed of the method in more detail, the mean relative error was
calculated asq2∑N (σdirect (ri ) − σFMM (ri ))
q
.
(23)
ε=
2
(ri )
∑N σdirect
In Equation 23, σdirect (ri ) denotes the source strength at collocation point i as computed by the direct EBEM
and σFMM the source strength as a result of the FMM. The resulting error as well as the required computation
times are shown in Figure 3. The behavior of the error for large p is unexpected. From theory, a monotonously
103

10−2

102

10−3
10−4

time [s]error ε [-]

10−1

101
10

5

15

20
25
expansion length p

30

35

40

Figure 3. Computing time and relative error of the FMEBEM over the expansion length

time [s]

102
102

101

101
100

direct EBEMFMEBEM
104

105

100
10−1

106

used memory [GiB]

decreasing error, that eventually converges towards zero would be expected. However, the error in Figure 3
converges towards a finite value. The reason for this behavior is suspected to be related to inaccurate handlingof very large numerical values occurring during the computation of the factorials in Equation 11 and 12 forlarge p. Since the overall accuracy is still good and the effect only occurs for impractically large expansion
lengths, it was not yet further investigated.
To analyze the performance of the method for varying problem sizes, the sphere was discretized with approximately 1,000 elements. This base mesh was refined several times by splitting the element edges until amaximum of approximately one million elements was reached. Figure 4 shows the computing time and required
memory for different numbers of DOF for both, direct EBEM and FMEBEM, using a constant expansion lengthof p = 5.

degrees of freedom N

Figure 4. Computing time and memory consumption during solution for EBEM and FMEBEM for different
problem sizes

5462


It can be observed that the EBEM does scale approximately with the complexity O N 2 while the FMEBEM
behaves as O (N). The reduced effort allows for faster solution and larger problem sizes due to reduced memory
consumption.
4.2 Car interior
Finally, a somewhat more complex model has been analyzed to demonstrate the applicability to more realistic
engineering problems. The mesh shown in Figure 5 was used to model a generic car interior excited by aloudspeaker system. It consists of approximately 40,000 constant elements. The blue circular area in the front
has a unit intensity defined as boundary condition, modeling a uniformly vibrating speaker membrane. For thegreen area at the bottom a high absorption coefficient of α f loor = 0.7 is defined. The remaining yellow elements
have a low absorption of α = 0.05 prescribed. While there are severe simplifications in comparison to an actual
car interior, the complexity of the model is high compared to the sphere since different types of boundary
conditions and a more complex geometry are involved. Solutions computed by EBEM (6a) and FMEBEM
using p = 5 (6b) are shown in Figure 6. The results match very well, though the relative error computed on the

Figure 5. Model of a generic car interior(a) EBEM solution(b) FMEBEM solution(c) error on the boundary(d) error inside the domain after post processing

Figure 6. Comparison of solutions for the car model obtained by EBEM and FMEBEM

5463

boundary (6c) becomes quite large in some areas. However, the areas, where the largest relative error occurs,
are those, where the EBEM solution approaches zero so even very small absolute errors result in large relative
errors. The evaluation of the field variables after post processing the results (6d) shows excellent agreement.

5

SUMMARY AND OUTLOOK

A fast multipole method for the undamped EBEM kernel function has been developed and analyzed. The
method has been verified by results of a direct EBEM and shows very high accuracy. The required memory
during the solution is reduced significantly and speed is vastly increased compared to the direct evaluation.
More efficient operators of complexity O p3 are known for the Laplace kernel [1, 10, 12]. They offer the
potential to improve the current implementation. The application of a black box FMM for the damped andundamped EBEM kernel [11] may also be of interest.

REFERENCES[1] Nishimura, N. Fast multipole accelerated boundary integral equation methods, Applied Mechanics Reviews
55 (4), 2002, pp 299-324[2] Bitsie, F. The Structutral-Acoustic Energy Finite Element Method and Energy Boundary Element Method,
Dissertation, Purdue University, 1996
[3] Le Bot, A. A Vibroacoustic Model for High Frequency Analysis, Journal of Sound and Vibration 211, 1998,
pp 537-554[4] Wang, A.; Vlahopoulos N.; Wu, K. Development of an energy boundary element formulation for computing high-frequency sound radiation from incoherent intensity boundary conditions, Journal of Sound andVibration 278, 2004, pp 413-436[5] Ram, M.; Keuchel, S.; Piepiorka, J.; von Estorff, O. Energiebasierte BEM und Fast-Multipole BEM fürHochfrequenzprobleme: Eine Gegenüberstellung, Fortschritte der Akustik - DAGA 2018, 19.-22. März 2018
in Munich[6] Ram, M.; Keuchel, S.; von Estorff, O. Beschleunigungsverfahren für die Energie-Boundary-ElementeMethode, Fortschritte der Akustik - DAGA 2019, 18.-21. März 2019 in Rostock[7] Gumerov, N. A.; Duraiswami, R. Fast Multipole Methods for the Helmholtz Equation in Three Dimensions,Elsevier Series in Electromagnetism, Elsevier Science, 2004
[8] Keuchel, S. Aufwandsreduzierung in der Fast-Multipole-Boundary-Elemente-Methode, Technische Universität Hamburg-Harburg, 2017[9] Keuchel, S.; Hagelstein, N. C.; Zaleski O.; von Estorff, O. Evaluation of hypersingular and nearly singular integrals in the Isogeometric Boundary Element Method for acoustics, Computer Methods in Applied
Mechanics and Engineering 325, 2017, pp 488-504[10] Yoshida, K.-I. Applications of Fast Multipole Method to Boundary Integral Equation Method, Kyoto University, 2001[11] Fong, W.; Darve, E. The black-box fast multipole method, Journal of Computational Physics 228 (23),2009, pp 8712-8725[12] Greengard, L.; Rokhlin, V. A New Version of the Fast Multipole Method for the Laplace Equation in
Three Dimensions, Research Report YALEU/DCS/RR-1115, 1996

5464