Evaluación del desempeño del método de elementos de contorno multipolar rápido para procedimiento de optimización topológica

**Carla T. M. Anflor**^{1} Luciana M. Braga^{1} Eder L. de Albuquerque^{1}

^{1}Universidade de Brasilia. Faculdade do Gama. CEP. 71405-610. PO Box 8114. Brasil. E-mail: anflor@unb.br; lucianam.braga@yahoo.com.br; eder@unb.br

**A****BSTRACT**

The objective of this work is to evaluate the performance of a Boundary Element Method (BEM) accelerated by the Fast Multipole Method (FMM) compared with a direct (BEM) in an optimization topology problem. The formulation of the Fast Multipole Boundary Element Method (FMBEM) is introduced in order to make the optimization process more attractive from the point of view of the computational cost. The fast multipole formulation is briefly summarized. A topological-shape sensitivity approach is used to select the points showing the lowest sensitivities, where material is removed by opening a cavity. As the iterative process evolves, the original domain has holes progressively removed, until a given stop criteria is achieved. For comparison, the topology resulting from a direct BEM optimization process is used. The CPU time x DOF's is also investigated. The accelerated BEM shows good performance in an optimization routine.

**Keywords:** Topology optimization, topological derivative, fast multipole method, boundary element method.

*RESUMEN*

*El objetivo de este trabajo es evaluar el rendimiento de un Método de elementos de contorno (BEM) acelerado por el Método Multipolar Rápido (FMM), en comparación con un BEM directo en un problema de optimización de topología. La formulación del Método de elementos de contorno multipolar rápido (FMBEM) se introduce con el fin de hacer que el proceso de optimización sea más atractivo desde el punto de vista del coste computacional. La formulación del método multipolar rápido se resume brevemente. Un enfoque al respecto de la sensibilidad topológica-forma es empleado para seleccionar los puntos que muestran las sensibilidades más bajas, a la cual se le retirará material mediante la apertura de una cavidad en el mismo. A medida que el proceso iterativo evoluciona, el dominio original tiene agujeros eliminadas progresivamente, hasta que se alcanza un criterio determinado de parada. Para la comparación, se utiliza la topología resultante de un proceso de optimización directa BEM. El tiempo de CPU x GDL también es investigado. El BEM acelerado muestra un buen rendimiento en una rutina de optimización.*

**Palabras clave:** Optimización de la topología, derivada topológica, método multipolar rápido, método de elementos de contorno.

**INTRODUCTION**

Although the Boundary Element Method (BEM) provides some advantages when modelling many problems, its efficiency is not suitable for large-scale 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).

(1)

where *?**/(?)* and *?**(?*_{?}) are the cost function evaluated for the original and the perturbed domain, respectively, andf 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).

(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).

Solve

(3)

The boundary conditions imposed to the external boundaries are subjected to equation (4):

(4)

where equation (5),

(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*,* while * * and 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.

Table 1. Topological derivative for the various boundary conditions prescribed 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).

(6)

Figure 1. Domain *O* and its boundary *r.*

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).

(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).

(8)

where *u*(x, y)* and *q*(x, y)* are the Green's function for 2D problems as equation (9).

(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).

(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).

(11)

where the *u,* and (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).

(12)

where A is the coefficient matrix, *X* 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].

Table 2. M2M, M2L and L2L expansions.

The main idea of the fast multipole BEM can be briefly described as:

- Step 1 - Discretization. The domain fl is discretized into boundary elements.

- Step 2 - Construction of the tree structure of the boundary element mesh. A square covering the discretized domain **r**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 Figure 2.

- Step 3 - Computation of moments on all cells. This step is also known as *upward pass.* The moments are computed on all cells. If a leaf is under consideration, moments are calculated directlyby using * *where *S*_{c} is the set of elements contained in the *leaf* and *z*_{c} is the centroid of the *leaf.* For the parent cell, the M2M translation is applied and the moments are summed on its four child cells. The M2M equation is presented in Table 2 where *z*_{c} is the centroid of the parent cell while the *z**c* represents the centroid of a child cell.

- Step 4 - Determination of cells with interactions. This step is named *Downward pass.* In this step a classification is done in order to define the distribution of all cells around a defined cell C. Adjacent cells are those level *l that* have at least one common vertex. Cells well separated at a level *l* are those who are not adjacent at level *l* but their parent cells are adjacent at level *l-1*.The interaction list of C is a list of all well-separated cells from a level *l*. Far cells of C are those whose parent cells are not adjacent to the parent cell of C. The local expansion associated with a cell C is calculated by the use of the M2L translation. The L2L translation is calculated for the parent cell of C with the expansion point being shifted from the centroid of C parent cell to that of C. Considering the cell C at level 2, the M2L translation is used to compute the coefficients of the local expansion.

- Step 5 - Evaluation of the integrals.

- Step 6 -Iterations of the solutions. The unknown solution vector *?* in the system *A ? = B* is updated by the iterative solver and continues for all levels to evaluate the subsequent matrix and vector multiplication *A* Xuntil the solution *?* converges to a defined tolerance.

Figure 2 depicts the basic idea of the steps involved in the FMM.

Figure 2. FMM Scheme.

**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).

Figure 3. Heat conductor boundary conditions.

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 (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 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.

Figure 4. Comparison between topologies with approximately 60% of volume.

Figure 5. Color map for case (e) and case (f).

Figure 6. Direct BEM and FMBEM CPU times x DOF.

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 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.

Figure 7. Internal Points: Evaluation by offset and over all domain.

**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.

**REFERENCES**

[1] L. Greengard. "The rapid evaluation of potentials fields in particle systems". Ed. The MIT Press, Cambridgue, UK. 1st edition, pp. 1-99. 1988. ISBN: 9780262571920.

[2] V. Roklin. "Rapid Solution of integral equations of classical potential theory". J. Comput. Phys. Vol. 60, Issue 2, pp. 187-207. September, 1985. ISSN: 0021-9991. DOI: 10.1016/0021-9991(85)90002-6.

[3] A. Peirce and J. Napier. "A spectral multipole method for efficient solutions of large scale boundary element models in elastostatics". Int J. Numer. Meth. Eng. Vol. 38, Issue 23, pp. 4009-4034. December, 1995. ISSN: 1097-0207. DOI: 10.1002/nme.1620382307.

[4] A. Mammoli and M. Ingber. "Stokes flow around cylinders in a bounded two-dimensional domain using multipole accelerated boundary element methods". Int J. Numer. Meth. Eng. Vol. 44, Issue 7, pp. 897-917. March, 1999. ISSN:1097-0207. DOI: 10.1002/(SICI)1097-0207 (19990310)44:7<897::AID-NME530>3.0.CO;2-S.

[5] Y. Liu and N. Nishimura. "The fast multipole boundary element method for potential problems: A tutorial". Engineering Analysis with Boundary Elements. Vol. 30, Issue 5, pp. 371-381. May, 2006. ISSN: 0955-7997. DOI:10.1016/j.enganabound.2005.11.006.

[6] Y. Liu. "A new fast multipole boundary element method for solving 2-D Stokes flow problems based on a dual BIE formulation". Engineering Analysis with Boundary Elements. Vol. 32, Issue 2. pp.139151. February, 2008. ISSN:0955-7997. DOI:10.1016/j.enganabound.2007.07.005.

[7] Y. Liu. "Fast Multipole Boundary Element Method: Theory and applications in Engineering". Cambridge University Press. New York, USA. 1^{st} edition, pp. 1-235. 2009. ISBN: 9780521116596.

[8] M. Dondero, A. Cisilino, J. Carella and J. Tomba. "Effective thermal conductivity of functionally graded random micro-heteregeneous materials using representative volume element and BEM." International Journal of Heat and Mass transfer. Vol. 54, Issue 17-18, pp. 3874-3881. August, 2011. ISSN: 0017-9310. DOI: 10.1016/j.ijheatmasstransfer.2011.04.04.

[9] R. Feijóo, A. Novotny, E. Taroco and C. Padra. "The topological derivative for the Poisson's problem". Mathematical Model and Methods in Applied Sciences. Vol. 13, Issue 12, pp. 1825-1844. December, 2003. ISSN: 1099-1476. DOI: 10.1142/S0218202503003136.

[10] L. Wrobel and M. Aliabadi. "The Boundary Element Method. Vol. 1 and Vol. 2: Vol. 1. Applications in Thermo-Fluids and Acoustics". Ed. Wiley, Chichester, UK. 1^{st} edition, pp. 1-468. 2002. ISBN: 978-0-471-72039-3.

[11] C. Anflor and R. Marczak. "Topological optimization of anisotropic heat conducting devices using Bézier-smoothed Boundary Representation". Computer Modeling in Engineering & Sciences (Print). Vol. 78, Issue 3-4, pp. 151-168. December, 2011. ISSN: 1526-1492. DOI: 10.3970/cmes.2011.078.151.

*R**eceived: August 25, 2014 Accepted: January 7, 2015*