The fast multipole boundary element method performance evaluation for topological optimization procedure

,


INTRODUCTION
Although the Boundary Element Method (BEM) provides some advantages when modelling many problems, its efficiency is not suitable for largescale models. The BEM, in general, produces dense and non-symmetric matrices that, in spite of being smaller in size, require O(N 2 ) operations to compute the coefficients (N is the number of equations of the linear system and O regards to operation). In order to solve the resulting system using direct solvers, other O(N 3 ) operations is also required. In order to overcome this inefficiency a coupling between the Fast Multipole Method and the BEM is proposed. This will allow solving problems with several million degrees of freedom (DOF's). Generally, the Finite Element Method (FEM) is indicated to solve models with several million DOF's. On the other hand, the BEM has been limited for solving problems with a few thousand DOF's for many years. In the last years, great efforts have been made to improve the computational cost of the BEM, maintaining its features, such as, easy of modeling, small matrices and no mesh dependency. The BEM can now be applied to large scale problems, as the topology optimization problem. Because this sort of problem is iterative, the number of elements increases as the material is removed and therefore a significant number of DOF's are introduced. The computational cost is a serious problem especially when the case under investigation is a 3D problem. During the last decades, many efforts have been done in order to accelerate the BEM for large-scale problems. The FMM was firstly presented by [1,2] with the aim of accelerating BEM solutions. The main goal was to reduce the CPU time in FMM accelerated BEM to O(N). Thereafter, this technique was applied for solving elasticity [3] and fluid [4] problems in large-scale. According [5] the FMM was considered one of the top algorithms in scientific computing developed in the 21 th century. In this publication, the authors developed a complete tutorial which presents the basic concept and the main procedures in the FMM for solving boundary integral equations for 2D potential problems. The FMM formulation was extended in [6] for large-scale analysis of two-dimensional (2D) Stokes flow problems. For solving the dual Boundary Integral Equation (BIE) formulation, [6] employed a linear combination for the velocity and the hipersingular BIE for traction to attain a better conditioning for the system of equations. Some examples were presented and showed the good accuracy and efficiency with the proposed approach. The book Fast Multipole Boundary Element Method [7] give many instructions in order to provide fundamentals to implement the method. The FMM was implemented by [8] for solving the effective thermal conductivity (ETC) of random micro-heterogeneous materials using representative elements and the FMBEM. The main goal of this paper is to implement the FMM in a topology optimization code. The idea relies on comparing the performance of both methodologies, i.e., optimization with Direct BEM against FMBEM with respect to CPU time and resulting topologies. This paper is organized as follow: The main idea of Topological Derivative (DT) is discussed and the respective analytical expressions for Poisson problems are presented. The BEM and the FMBEM for 2D potential problem are shown. In the sequence some numerical examples and their results are presented. Finally this work is concluded and some discussions are carrying on.

TOPOLOGICAL DERIVATIVE
A topological derivative (DT) for the Poisson equation is applied in this work for determining the domain sensitivity. A simple example of applicability consists in a case where a small hole of radius (ε) is open inside the domain. The concept of topological derivative consists in determining the sensitivity of a given cost (ψ) function when this small hole is increased or decreased. The local value of the DT at a point ( ) inside the domain for this case is evaluated by equation (1).
where ψ(Ω) and ψ(Ω ε ) are the cost function evaluated for the original and the perturbed domain, respectively, and f is a problem dependent regularizing function. By equation (1), it is not possible to establish an isomorphism between domains with different topologies. This equation was modified introducing a mathematical idea that the creation of a hole can be accomplished by simply perturbing an existing one whose radius tends to zero. This allows the restatement of the problem in such a way that it is possible to establish a mapping between the two domains [9], as equation (2).
where δε is a small perturbation on the holes radius.
In the case of linear heat transfer, the direct problem is stated as equation (3).
The boundary conditions imposed to the external boundaries are subjected to equation (4): where equation (5), Regarding to equation (5), h is a function which takes into account the type of boundary condition on the holes to be created. It means that if h(0, 0,1) a boundary condition of Robin is imposed on the holes. The temperature and flux on the hole boundary are u ε , ∂u ε ∂n = q ε , while u ∞ ε and h c ε are the hole internal convection parameters, respectively. After an intensive analytical work, an explicit expressions for DT were developed for problems governed by equation (3). Table 1 summarizes the final expressions for the topological derivative, considering the three classical cases of boundary conditions on the holes.

THE BOUNDARY ELEMENT METHOD
A brief review on the boundary element method using constant elements is summarized in this work. An initial domain depicted in Figure 1 is established for prescribed boundary conditions, considering a Laplace equation governing a 2D potential problem presented as equation (6). For a potential problem, three kinds of boundary conditions may be imposed: Dirichlet, Neumann and/or Robin. For this presentation, the first and second boundary conditions are imposed as equation (7).
where u is the potential field in domain Ω, Γ is the boundary of Ω, n is the outward normal. Note that the barred quantities are the values imposed by the boundary conditions. The solution of equation (6) under boundary conditions equation (7) is presented as equation (8).  where u * (x, y) and q * (x, y) are the Green's function for 2D problems as equation (9).
∂r ∂η (9) where r represents the distance between the collocation point x and the field point y, as depicted in Figure 1. Taking x to the boundary, the classic BIE formulation of BEM [10] is obtained as equation (10).
If the boundary is smooth at the collocation point x, C(x) = 1/2. The next step consists in discretizing the boundary Γ using N constant elements. The discretized equation of BIE is now presented as equation (11).
where the u j and q j (j = 1,2,…,N) are the nodal values of the u and q at the element ΔS j , respectively. Applying the boundary conditions (7) at each node and switching the columns for grouping the unknown variables, one finds the equation (12).
where A is the coefficient matrix, λ the unknown vector and B the known right-hand side vector.

THE FAST MULTIPOLE BOUNDARY ELEMENT METHOD
The BEM uses the Green's functions as the weighting function on its formulation which increase the accuracy when compared with others numerical techniques [11]. As a result, the spatial dimension is reduced by one. Additionally, the computational cost of a traditional direct BEM can be reduced by using the FMBEM. The goal of FMM relies on translating node-to-node interactions to cell-to-cell interactions. These cells have a hierarchical structure called tree while the small ones are called leaves. The FMM employs iterative equation solvers (GMRES) where matrix-vector multiplications are calculated using fast multipole expansions. As iterative equations are used, some parameters for the FMM, such as maximum number of elements allowed in a leaf (maxl) and in the tree structure (levmax), number of terms in multipole expansion (nexp) and local expansion (ntylr), and also the GMRES tolerance (tol) must be set. Expansions used for 2D potential problem for the FMM are briefly summarized in Table 2. Further details about the analytical derivations can be found in [5,7]. The main idea of the fast multipole BEM can be briefly described as: o Step 1 -Discretization. The domain Ω is discretized into boundary elements. o Step 2 -Construction of the tree structure of the boundary element mesh. A square covering the discretized domain Γ is considered. This square is classified as a cell of level 0. This parent cell is divided into four child cells, now classified as level 1. This procedure is iteratively done until a stop criteria is achieved. This criteria is achieved when the number of elements imposed by the user in that cell is reached. A cell having no child cells is call leaf, which are in grey in

NUMERICAL RESULTS
The high computational effort involved in an optimization process motivates the implementation of the FMM in order to maintain those attractive characteristics when coupling BEM and DT [10]. This section presents one example that demonstrates the application of the proposed method. The results obtained for the FMBEM for each case are compared with Direct BEM. During the optimization process, the computational cost, number of DOFs and volume were taken into account. For a specific iteration, the respective intermediary topology is illustrated. The iterative process was halted when a given amount of material was removed from the original domain. In all cases, the total potential energy was used as the cost function. A regularly-spaced grid of internal points was generated automatically, taking into account the radius of the holes created during each iteration. The radius was obtained as a fraction of a reference dimension of the domain (r = ω l ref ). In all cases l ref = min (height, width) was adopted. The objective in all cases was to minimize the material volume. The current volume of the domain (V f ) was checked at the end of each iteration until a reference value was achieved (V f = φ V 0 , where V 0 represents the initial volume and φ a defined percentage of material to be removed).

Heat Conductor
This example refers to a square domain subjected to low temperature boundary conditions (BC) on its corners and a decentralized high temperature BC at the left surface. The problem is illustrated as Figure 3, where T H is the high temperature (373 K) and T L is the low temperature (273 K).
The remaining boundaries and the holes opened during the optimization process are insulated. The stop criteria was set when V f = 0.6 V 0 is reached. In order to evaluate the different resulting topologies due to the FMBEM parameters, five cases with distinct set up are considered. All cases (case (a), case (b), case (c), case (d) and case (e)) presented in Figure 4 are always compared with the topology resulting by using the direct BEM, namely case (f). For the case Figure 2. FMM Scheme.
both cases, Figure 5. Also, it is important to note that both topologies attained the same volume at the same iteration. The CPU time x DOF for case (e) and case (f) are presented in Figure 6. During the optimization process, a maximum number of approximately 2500 elements were evaluated. It is possible to verify that the performance of the FMBEM is superior when a large number of elements is used. An interesting evaluation relies on determining the intersection between the curves. This intersection allows determining exactly the number of DOF in which the FMBEM overtakes the direct BEM in efficiency.
As the main goal of this work regards to decrease the computational cost, some additional computational artifices in the numerical routine were also employed. One of these artifices regards to reduce the internal grid of internal points. In this sense the present code was written in order to generate internal points only near the boundaries (offset) or in the domain.
Obviously, when dealing with domains with a significant internal sensitivity, an evaluation on all domain is required. A good recommendation is to use both numerical artifices, i.e. some iteration with an offset of internal points and a predefined intermediary iteration which takes into account a  (f) constant elements were used and all integrations were performed analytically, characterizing it as a benchmark test. The comparison on performance between both methods is only carried out when the FMBEM final topology matches with that one stressed in the benchmark test. Figure 4 shows the topological results using different parameters (tol, maxit, maxl, levmx, nexp, ntylr) in the FMBEM. The first four topologies showed a slight difference when compared with case (f), due to the parameters of FMBEM employed. From now on, only case (e) and case (f) will be considered. Case (e) produced a topology that matches perfectly with that resulting by using Direct BEM, and for this reason it is possible to compare the temperature distribution in  complete grid over the domain, see Figure 7. It is also important to remark that this procedure is not possible in the finite element method due to its features of domain mesh. The feature of absence of domain mesh makes BEM more suitable for topological optimization than the FEM.

CONCLUSIONS
The computational cost is an important Issue in any numerical analysis. Regardless of the many attractions of the boundary element method, the technique is not widely used in commercial codes because of the high computational cost for solving large-scale problems. This problem increases when considering an iterative optimization process, where the problem must be evaluated several times. In order to overcome this difficulty, the FMBEM was implemented in a topological optimization code. The resulting topology of a benchmark test using the FMBEM was compared with the final topology obtained by using the direct BEM. The CPU time for both cases was compared until the final topologies have been reached. The final topologies for case (e) and (f) have shown good agreement once the FMBEM parameters were adjusted. The results suggest that the direct BEM is more efficient for problems with a coarse discretization, i.e., smaller number of DOF's. As the iterative process evolves, the number of elements increases significantly and the FMBEM overtakes the direct BEM in efficiency. Some remarks must be taken into account, such as; while the iterative process does not reach around 2500 DOF's, the direct BEM is recommended. When this number of DOF's is exceeded a flag (previously implemented in the code) must be turned on so that the process is carried on using the FMBEM. Another interesting conclusion relies on the fact that the final shape of the resulting topology depends significantly on the parameters set for the FMBEM. Finally, the use of the FMBEM allows very well refined topologies to be treated without needing parallel computation.