Flexible multi-level fast multipole BEM with direct solver for industrial
acoustic problems
Yue LI1,2∗ , Onur ATAK1 , Wim DESMET2,3
1 Siemens Industry Software, Interleuvenlaan 68, B-3001 Leuven, Belgium
2 Noise and Vibration Research Group, PMA, KU Leuven, Celestijnenlaan 300 B, B-3001, Heverlee, Belgium
3 DMMS core lab, Flanders Make, Belgium

Abstract
Fast Multipole Boundary Element Method (FMBEM) for acoustics is a very powerful method to solve ultra
large problems, which exploits iterative solvers. However, for Helmholtz problem iterative solvers may fail toconverge in some cases. Whereas FMBEM provides a unique offering for ultra large problems, many industrial
problems in fact can be solved with tens of thousands of Dofs and can be considered as medium to large size. For
such problems, direct solvers are still viable and provide a good option in terms of stability and convenience forsolving many right-hand sides, especially in light of recent many core architectures like GPU. In this paper, we
present a flexible fast multipole accelerated indirect BEM for solving medium to large size acoustic problems.
The assembly process of the indirect BEM is accelerated by a flexible multi-level fast multipole method. The
system matrix is then solved by a direct solver. Numerical results show that this combination offers an efficientand robust approach, where medium to large size acoustic problems can be efficiently solved on a modern
desktop computer.
Keywords: Flexible FMBEM, Acoustics, Direct solver

1

INTRODUCTION

Boundary Element Method (BEM), sometimes referred to as boundary integral equation method (BIEM) or
boundary integral solutions, is a well-known numerical method in many engineering applications. The uniqueadvantage of BEM over other volume based numerical methods, e.g. Finite Element Method (FEM), Finite
Volume Method (FVM), is that BEM reduces the dimensionality of the problem by one. It only requires a
two-dimensional discretization for a three-dimensional problem. This advantage brings great benefit to Computational Aided Engineering (CAE) applications considering that the volume meshing process can be a significant
cost in practice. Besides, unbounded domains are common for exterior acoustic applications. BEM is naturally
convenient in computing unbounded problems as the Kirchhoff-Helmholtz integral equation satisfies the Summerfield radiation condition automatically. The well-known limitation of conventional BEM (CBEM) lies in its
fully populated coefficient matrix. The numerical complexity of CBEM is proportional to O(N 2 ) in assemblyand O(N 3 ) in solving with conventional direct solvers, where N denotes the number of degree of freedoms(Dofs) in the model. Over the years, fast techniques such as fast multipole method [1], H -matrix [2], precorrected FFT, have been developed to accelerate the conventional BEM. These techniques most often exploit
iterative solvers with the consideration of saving memory cost for extremely large cases. Most of the existing
FMBEM literature contribute to the direct BEM (dBEM) formulation [3], which limits the applications to closed
boundary conditions. It should be noted that dBEM can only solve the Helmholtz equation in either interior orexterior domain at a time. Whereas various industrial cases require the modelling of open boundaries and/or
combined interior-exterior problems. Indirect Galerkin BEM (iBEM) formulation is able to handle such cases.
The attractive feature of employing a Galerkin approach is that the system of equations of the boundary element
formulation becomes symmetric which reduces the computational cost by half.
∗ li.yue@siemens.com

7518

Moreover, the existing FMBEM solution coupled with iterative solvers provides a unique offering for ultra
large problems. Many industrial problems in fact can be solved with tens of thousands of Dofs and can beconsidered as medium to large size. For such problems, direct solvers are still viable and provide a good optionin terms of stability and convenience of solving many right-hand sides. In addition, recent advances in many
core architectures, such as GPUs, make dense solvers competitive again [7].
In this paper, we investigate the possibility of applying direct solver to FMBEM. We present a fast multipole
accelerated indirect BEM for solving medium to large size industrial problems. Instead of using iterative solver
in FMBEM, we assembly the full system matrix during the computation in order to benefit from the robust
direct solvers. The contributions of this paper lie in two aspects:
1) For the first time, a flexible multi-level structure is applied to fast multipole BEM using indirect formulations;2) We apply flexible fast multipole method to accelerate BEM matrix assembly. The system matrix is eventually
solved by a robust direct solver. Numerical results have shown the advantages of our approach for the target
problems.

2

THEORY

2.1 Indirect boundary integral formulation for Helmholtz problems
For time harmonic acoustic problems, the steady acoustic pressure p at any location x in a three dimensional
fluid domain V is governed by Helmholtz differential equation∇2 p(xx) + k2 p(xx) = 0, ∀xx ∈ V(1)where k = ω/c is the acoustic wave number depending on angular frequency ω and speed of soundc.
In indirect BEM formula, the boundary variables are defined as double layer potential µ and single layer potential σ , which are defined as the difference of the pressure and the difference of the normal derivative of thepressure respectively:µ(yy) = p1 (yy+ ) − p2 (yy− )σ (yy) =∂ p1 (yy+ )
∂ ny+

−

∂ p2 (yy− )∂ ny−= − jρω[q1 (yy+ ) + q2 (yy− )]

(2)
(3)

where the subscripts 1 and + are associated with the direction of the unit normal, 2 and − are associated withthe opposite direction of the unit normal.
Using the single layer and double layer potentials, we can write the indirect Kirchoff-Helmholtz integral equation as follows:Z∂ G(xx, y )p(xx) = [µ(yy) − G(xx, y )σ (yy)] dSy , ∀xx, y ∈ S
(4)
∂ ny
S
where ny denotes unit normal at a source point y . G(xx, y ) is the Green function of the Helmholtz equation. Its
3D form is given as
G(xx, y ) =eik|xx−yy|
4π | x − y |(5)

The associated boundary conditions of Helmholtz equation can be written for both sides as well. Assuming athin boundary surface, it leads to the following boundary condition expressions:µ(yy) = 0,
σ (yy) = 0,p(yy) = p(yy) on S p(6)∂ p(yy)= − jρ0 ωvn (yy) on Sv
∂n

(7)

7519

σ (yy) = − jkβ µ(yy),∂ p(yy)= − jkβ p(yy) on SZ∂n(8)

At each position on the boundary sureface, either the double layer potential is zero or the single layer potential
is zero, or the single layer potential is related to the double layer potential. From these relations the indirectboundary integral formulation can be reformulated as:p(xx) = −Z
Sp

∂ G(xx, y)∂ G(xx, y)dSv ++ jkβ (yy)G(xx, y )] dSZµ(yy)[
∂n
∂n
SZ
Z

Z

σ (yy)G(xx, y )dS p +

µ(yy)
Sv

(9)

Subsequently, Eq.(9) can be used to rewrite the boundary condition in equations Eq.(6)-Eq.(8). They can be
written in a compact form by using a function f of unknowns σ and µ.
f p (σ , µ) = p on S p(10)

fv (σ , µ) = − jρ0 ωvn on Sv(11)fZ (σ , µ) = 0 on SZ(12)

A variational approach is adopted to avoid Hadamard finite part integrals and to achieve symmetric system
matrices. The equivalent variational statement of the equations Eq.(10)-Eq.(12) can be written as [6]:∀(δ σ , δ µ) :

Z
Sp

Z

f p (σ , µ)δ σ dS p +Sv

Zfv (σ , µ)δ µdSv +Z

SZfZ (σ , µ)δ µdSZ =Sppδ σ dS p −Z
Svjρ0 ωvn δ µdSv (13)2.2 Fast multipole method
The diagonal fast multipole form is proposed by [1]. The basic idea of fast multipole method is to convert the
point-to-point operation to a cluster-to-cluster operation. Considering two clusters of points as shown in Fig.1,
x0 and y0 are the centers of the source cluster C1 and the target cluster C2 respectively. When C1 and C2 are
well separated, one could convert a direct evaluation from x j to yl , to a multipole evaluation which consistsof three main steps: 1) Upwardpass: the contribution of x j is firstly computed to x0 , the center of the source
cluster C1 ; 2) Multipole translation: compute the multipole translation from source cluster’s center to target
cluster’s center, from x0 to y0 in this case; 3) Downwardpass: distribute the contribution from center y0 to yl inthe cluster C2 .

x0

y0

D

yl

xj

C2
C1
Figure 1. Clusters of x j (C1 ) and yl (C2 )In mathematical perspective, the basic solution of the Helmholtz equation (Eq.(5)) can be rewritten using Gegenbauer addition theorem and plane wave expansion form:Dkeikkx j −yl keikkdd +D
ik==
4π k x j − yl k 4π k d + D k 16π 2

The diagonal translation operator ML is defined as

7520

Z
S2

eik(yl −y0 )·ss ML (ss, D )e−ik(x j −x0 )·ss ds

(14)

L−1

(1)

D)
ML (ss, D ) = ∑ (2l + 1)il hl (kD)Pl (ss · D̂

(15)

l=0
(1)

where d = (yl − y0 ) − (x j − x0 ), D = y0 − x0 , d =k d k, D =k D k, hl (·) is a spherical Hankel function of theˆ = (·)/(k · k). s
first kind of order 1, Pl (·) is a Legendre polynomial. Normalized vectors are indicated by (·)
is the unit sphere given by s (θ , φ ) = (sinθ cosφ , sinθ sinφ , cosθ ). 0 ≤ θ ≤ π, 0 ≤ φ ≤ 2π are polar angle andazimuthal angle in the spherical coordinate system.
The integral over the unit sphere surface S can be evaluated numerically: in the θ direction by using a Gausslegendre integration with L points and in the φ direction by using a trapezoidal rule with 2L-points.
Z 1 Z 2πZS2f (s)ds =
−1 0

L 2L

f (s)dφ d(cos(θ )) ≈ ∑ ∑ ωi
i=1 j=1πf (si, j )L(16)with ωi the weight of the Gauss-Legendre integration for the ith integration point in the θ -direction.
In numerical implementation, the sum of l is truncated at l = L. A semi-empirical rule for the expansion length
is given by many papers
L = kdl + p log(kdl + π)

(17)

where p specifies the required precision, dl is the maximum diameter of the clusters.
2.3 Flexible multi-level FMiBEM
When applying FMM to BEM, a partition algorithm is required to cluster the elements of the model into groups.
This can either be done on one level structure resulting in a single-level scheme, or on many levels resulting
in a multi-level scheme. In single-level FMBEM, the number of the clusters should be approximately equal to
the root of the number of elements to ensure the efficiency of the algorithm [5]. As the size of the problem
increases, the number of the elements in the near field clusters increases proportionally to the root of the total
number of elements, which makes the single-level FMBEM less efficient. Multi-level scheme tackles this issue
by employing a hierarchical structure. The number of elements in the cluster on the leaf level (lowest level)
becomes independent of the total number of elements. The size of the model is then reflected by the number oflevels in the tree structure. However, the cluster size increases from lower to higher level. Thus the expansion
order increases due to its dependence on cluster diameter. The computational cost on upper levels increases
quadratically with the expansion order with the consideration of employing more spherical Gauss points in θand φ directions in Eq.(16). This becomes significant when the model size is getting larger, where more levels(5-10) are required.
Many of the existing work tackle this issue by employing an interpolation/filtering algorithm among different
levels. Various schemes have been proposed in the literature [4]. The complexity of such algorithm is crucialin Helmholtz (and Maxwell) problems to ensure that it never overtakes the complexity of fast multipole method
itself. Nevertheless, this adds extra costs to the entire computation. A direct multi-level approach is presented
in [3], where the far field contributions are computed and stored on each level explicitly in order to avoid
the interpolation/filtering procedure. By sacrificing some memory cost, a Burton-Miller collocational BEM withdirect multi-level approach achieves better efficiency over the one with the interpolation/filter algorithm [3].
Inspired by the direct multi-level concept, we further develop a flexible and efficient a priori Octree structure
for symmetric Galerkin indirect BEM.
As used by many FMBEM papers, the original Octree structure always starts with one bounding box, which
generally is an over-determined box. Marked as level 0, the parent box is further divided into 2d child boxes
forming the next level. The procedure is repeated recursively till the predefined criteria is fulfilled. The upper
levels of the original tree are determined by the structured partition. In contrast, a priori Octree starts with a7521

Figure 2. Cube mesh (17k Dofs) with field points (up left); Pressure plot of CBEM (up right); Pressure plot of
FMBEM (down left); Pressure plot of f-FMaccBEM (down right).
level which has been computed a priori. Unlike the immutable structure in the original Octree, both the leaf
level and the top level in a priori Octree are flexible and controllable, which provide an efficient multi-levelstructure for fast multipole computation. The expansion order on the top level is maintained as low as possible
to ensure the efficiency of the ML calculation, and to avoid the interpolation/filtering algorithm at the same time.
The leaf level is computed a priori to fulfill the accuracy requirement from multipole assumption. In the highfrequency fast multipole formula, it is well-known that when kD becomes small (e.g. kD < 0.09 for expansion
order 8), the spherical Hankel function can easily blow up to the magnitude of 1016 , which is the machine
limit that double precision can achieve. A priori Octree allows a flexible construction of the tree structure. By
estimating the box size on the leaf level a priori, it restrains the Hankel function within a reasonable range.
2.4 Hyper-singularity
Hyper-singularity appears in both indirect BEM and Burton-Miller direct BEM formulas. There exists several ways to handle hyper-singularity, which includes singularity cancellation techniques, specialized quadrature
methods and regularization techniques. The presented paper focuses on the application of practical regularization
techniques to FMiBEM. Stokes’ theorem is applied to transform hyper-singularity to singular and weakly singular form. Since the Octree structure has split our problem into near fields and far fields, the hyper-singularity
shall be considered separately. Apparently, there is no singularity exist in the far fields. Thus the far field
contributions are computed directly by the original hyper-singular form, while the near field is transformed
by applying regularization technique. However, the regularization technique is only valid for closed boundary
condition, or when the principle term of hyper-singularity is zero. In practical implementation, this introduces
additional calculations on the edges in the near field clusters. Special care must be paid to the clusters whichcontains shared edges.

7522

2.5 Direct solver
Most of the existing FMBEM papers exploit various iterative solvers to solve ultra large problems. However,for Helmholtz problems, iterative solvers often have slow convergence, even may fail to converge in certain
cases. Finding an efficient preconditioner remains challenging. Whereas FMBEM provides a unique offering forultra large problems, many industrial problems are in fact can be solved with tens of thousands of Dofs. For
such problems, direct solvers are still viable and provides a good option in terms of stability and convenienceof solving many right-hand sides. In addition, dense solvers can now be efficiently accelerated by recent many
core architectures such as GPUs [7]. Therefore, we investigate the possibility of combining FMBEM with direct
solvers. We apply fast multipole method to accelerate the conventional BEM assembly procedure. Instead ofusing low rank approximation and iterative solvers, we compute the far field contributions for each cluster pairs
and assembly back to the global system matrix. In this sense, the global system matrix can be solved by direct
solvers. A dense Cholesky solver is used in our study. Benefiting from the symmetric feature of Galerkinapproach, we save half of the far field computations, as well as the near field computations. Based on our tests,
modern laptops/workstations are able to provide sufficient memory to handle the target size problems.

3

NUMERICAL EXAMPLESAll the following calculations are conducted using single thread CPU on a laptop with i7-7920HQ CPU and
32GB memory in Windows 10. In the following discussions, CBEM denotes conventioanl BEM, FMBEM
denotes conventional fast multipole indirect BEM with iterative solver, f-FMaccBEM denotes flexible fast multipole accelerated indirect BEM with direct solver. We choose Cholesky solver for f-FMaccBEM, and Generalized
Conjugate Residual (GCR) solver with Sparse Approximate Inverse (SPAI) preconditioner for FMBEM.
3.1 3D cube caseFirstly we consider a 3D cube case with prescribed velocity boundary condition. Excitation frequency is 800Hz,
giving kA = 14.8. The cube is discretized by linear triangular elements. Figure 2 shows the mesh of 17k Dofs
together with field points. The contour plots of the field points’ absolute pressure from CBEM, FMBEM andf-FMaccBEM are shown using the same scale in the same figure.

Figure 3. Assembly time of CBEM andf-FMaccBEM, and direct solver time

Figure 4. Total computational time for CBEM,FMBEM and f-FMaccBEM

Figure 3 plots the assembly time for CBEM and FMaccBEM together with the time of Cholesky direct solver.
It can be found that the assembly complexity of f-FMaccBEM follows NlogN compared to N 2 of the CBEM.

7523

For the model size of 38k Dofs, we see a speed up factor of 17 from CBEM (1.1 hours) to f-FMaccBEM(3.7 minutes). It is of interest to observe that the cost of assembling far field matrix and inserting it backinto the global system matrix doesn’t influence the total complexity. For model size up to 20k Dofs, assembly
cost is dominant in the total f-FMaccBEM computation. Significant speedup is achieved by f-FMaccBEM over
CBEM. After 20k Dofs, the cost of Cholesky solver is taking over and becomes the dominant cost. For larger
models, parallel computing techniques [7] can be applied to direct solver to obtain a good overall efficiency.
Figure 4 shows the total computational time of CBEM, FMBEM and f-FMaccBEM. For small to medium size
problem, we find the benefit of combining f-FMaccBEM with direct solver. It outperforms CBEM and theconventional FMBEM with iterative solver. Meanwhile, it should be noted that for a simple geometry in thistest, the iterative solver can easily reach target tolerance within 15-20 iterations. For more complex geometries
at high frequencies, iterative solvers may struggle to converge, which will be demonstrated in the next case.
3.2 Realistic car case
A Chrysler Neon car is studied below. The car mesh consists of 40,234 linear triangular elements and 20,010
nodes. As shown in figure 5, the engine of the car has been modelled in details. In this study we want toevaluate the frequency response function (FRF) of the engine noise. The engine is described by a unit velocity
boundary condition, while the rest of the car mesh is sound hard boundary. Typically, the peak part of engine
noise concentrates in the range 100 − 400Hz. Figure 6 plots the FRF computed by CBEM and f-FMaccBEM.
The receiver point is placed 50cm away from the driver’s door. By taking CBEM as reference, the relative L2error of f-FMaccBEM is well controlled under 0.5% for the entire FRF.

Figure 5. Car mesh
Figure 6. Frequency response function
Table 1. Computational time for the entire FRF
SolverCBEM
FMBEM
f-FMaccBEM

Assembly time
485
91

Solve time
190
190

Total time [mins]
675
400
281

Table 1 gives the total computational time (exclude post processing) for the entire FRF. f-FMaccBEM saves
60% of time compared to CBEM, and 30% of time compared to FMBEM. Table 2 lists some detailed costs of
FMBEM and f-FMaccBEM at different frequencies. f-FMaccBEM outperforms FMBEM at all frequencies. It
can be seen that for higher frequencies, the iterative solver in FMBEM is getting difficult to converge. By using
six element per wave length, the excitation frequency is valid up to 820Hz. Although equipped a preconditioner,7524

the iterative solver requires 549 iterations to converge to the predefined tolerance 0.001 at 800Hz. In contrast,
the direct solver used by f-FMaccBEM remains efficient for all valid frequencies.
Table 2. Computational time of FMBEM and f-FMaccBEM
Solver
FMBEM
FMBEM
FMBEM
FMBEM
f-FMaccBEM
f-FMaccBEM
f-FMaccBEM
f-FMaccBEM

4Frequency
100Hz
200Hz
400Hz
800Hz
100Hz
200Hz
400Hz
800Hz

Assembly time [s]
219
155
164
323

Solve time [s]
358
368
365
342

Total time[s]
863
707
824
1520
577
523
529
665

Nr of iterations
40
75
238
549
-

CONCLUSION

In this paper, we present a fast multipole accelerated indirect BEM for solving medium to large size acoustic
problems. We investigate and successfully combine a flexible multi-level FMBEM with direct solver. The
assembly process of the indirect BEM is accelerated by a flexible multi-level fast multipole method. The system
matrix is then solved by a direct solver. Numerical results have shown that this combination offers a robust andfast solution, where medium to large size acoustic problems can be efficiently solved on a modern desktop
computer. Future work will report on taking advantages of dense solver on modern hardware architectures suchas GPUs.

ACKNOWLEDGEMENTS
This project has received funding from the European Union’s Horizon 2020 research and innovation programme
under grant agreement No 721536.

REFERENCES[1] Rokhlin, Vladimir. Diagonal forms of translation operators for the Helmholtz equation in three dimensions,Applied and Computational Harmonic Analysis, Vol 1 (1), 1993, pp 82-93.
[2] Hackbusch, Wolfgang. Hierarchical matrices: algorithms and analysis. Vol. 49. Springer, Heidelberg, 2015.
[3] Z.-S. Chen, H. Waubke, W. Kreuzer, A formulation of the fast multipole boundary element method (fmbem)
for acoustic radiation and scattering from three-dimensional structures, Journal of Computational Acoustics,
Vol 16, 2008, pp 303–320.
[4] E. Darve, The fast multipole method: numerical implementation, Journal of Computational Physics, Vol 160,
2000, pp 195–240.
[5] K. Giebermann, Fast summations methods of numerical solution of integral equations for scattering problems
in R3, PhD thesis, 1997.
[6] T. Wu, Boundary Element Acoustics, Fundamentals and Computer Codes, WIT Press, 2000.
[7] Dongarra J, Gates M, Haidar A, Kurzak J, Luszczek P, Tomov S, Yamazaki I. Accelerating numerical dense
linear algebra calculations with GPUs. In Numerical computations with GPUs, 2014, pp. 3-28.

7525