FACTA UNIVERSITATIS Series: Mechanical Engineering Vol. 16, N o 3, 2018, pp. 405 - 417 https://doi.org/10.22190/FUME180912034B © 2018 by University of Niš, Serbia | Creative Commons License: CC BY-NC-ND Original scientific paper  EFFICIENT CALCULATION OF THE BEM INTEGRALS ON ARBITRARY SHAPES WITH THE FFT UDC 519.6:539.3 Justus Benad Berlin University of Technology, Berlin, Germany Abstract. This paper builds upon the results of a recent study which illustrates how the Fast Fourier Transformation (FFT) can be used to accelerate the Boundary Element Method (BEM) for arbitrary shapes. In the present work, we further deepen this understanding and focus especially on implementation details in order to calculate the boundary integrals with the FFT. Different numerical techniques are compared for an exemplary shape. Also, additions to the concept are mentioned such as the introduction of a high-resolution grid close to the boundary and a low-resolution grid farther away. Key Words: Laplace Equation, Navier Equation, Boundary Element Method, FFT 1. INTRODUCTION Over the past years, the Boundary Element Method (BEM) has become a highly efficient tool in contact mechanics for the calculation of the deformations and stresses with the half-space approximation [1, 2]. The boundary integral equations are simple convolutions over the surface of a half-space which can be evaluated rapidly with the Fast Fourier Transformation (FFT) [3, 4]. For many practical applications, however, the complex geometry of the contacting bodies exceeds the limitations of the half-space approximation. One of these fields is the modeling of deformations and stresses in bio mechanics, where the complex geometries of joints such as the hip and the knee joint (see for example [5, 6]), clearly exceed the range of application of the half-space theory. Another field is given by a wide range of aerospace applications, such as the simulation of deformations and stresses in highly stressed parts of aircraft engines. Turbine blades and discs are among these critical components [7, 8]. Highly undulating surfaces, cooling holes, and other complex geometrical elements in close proximity to the contact regions of turbine blade fir-tree connections cause difficulties when using the half-space theory to model deformations and stresses [9]. Therefore, other modeling techniques Received September 12, 2018 / Accepted November 10, 2018 Corresponding author: Justus Benad TU Berlin, Institut für Mechanik, FG Systemdynamik und Reibungsphysik, Str. d.17. Juni 135, 10623 Berlin, DE E-mail: mail@jbenad.com 406 J. BENAD have to be used for these kinds of problems, such as the Finite Element Method (FEM) or the classical BEM for arbitrary shapes [10]. Generally, this comes at the cost of a much higher computational complexity. Many valuable methods have been proposed in the past to reduce the computational complexity of both the FEM and the classical BEM for arbitrary shapes [10-14]. A recent study [15] illustrates how the classical BEM for arbitrary shapes can be accelerated with the FFT in a similar way as it is done with great success in contact mechanics within the framework of the half-space approximation. The study highlights that although the boundary integrals of the BEM are not simple convolutions over the boundary of arbitrary shapes, they can still be regarded as convolutions over the space which is one dimension higher than that of the boundary. This makes the evaluation of the boundary integrals with the FFT on the surface of an arbitrary shape possible if the dimension of the FFT is increased. For a three-dimensional shape with n 2 surface points, the computational complexity is O(n 4 ) to invert the classical BEM matrix. When the boundary integrals are evaluated with the FFT, this complexity is reduced to O(n 3 log n 1.5 ). This reduction along with highly efficient implementations which are available for the FFT make this technique appealing. In the present work, we will build upon the results found in [15] and further deepen the understanding of how the technique can be implemented. In an exemplary simulation, different numerical techniques are compared for a curved shape. Also, we mention additions to the method which further decrease the computational complexity. The overall goal of this work and future works will be to strive towards a complexity comparable to that of the FFT on the half-space which is O(n 2 log n). 2. INTERPRETATION OF THE BOUNDARY INTEGRAL FORMULATIONS AS CONVOLUTIONS The boundary integral formulation of the LAPLACE equation ( ) 0u x  (1) is * * 0 0 0 0 ( ) ( ) ( ) ( ) ( ) ( ) S S c x u x q x u x x dS u x q x x dS     . (2) It relates the values for ( )u x and ( ) ( ) ( )q x n x u x  on boundary S of a region with outward normal vector n to the value of u at a certain point 0 x in the region. It is c = 1 when 0 x lies inside the region, and c = 1/2 when point 0 x lies on the smooth boundary. Term * 0 ( )u x x is the fundamental solution of Eq. (1). For instance, it is  * 0 0 1 ( ) ln 2 u x x x x      (3) for the two-dimensional case of x yx xe ye  . For this case boundary S in Eq. (2) is a line, whereas for the three-dimensional case with x y zx xe ye ze   boundary S is a surface. The boundary integral formulation of NAVIER’s equation 1 1 ( ( )) ( ) ( ) 1 2 u x u x b x          (4) is Efficient Calculation of the BEM Integrals on Arbitrary Shapes with the FFT 407 * * 0 0 0 0 ( ) ( ) ( ) ( ) ( ) ( ) , S S c x u x t x u x x dS u x t x x dS       (5) where c is defined as above, u and t n  are the values on boundary S and * t is known from the fundamental solution 0 0 3 0* 0 0 3 4 ( ) ( ) ( ) 16 (1 ) x x x x I x x x x u x x              (6) with the material law. Both of the boundary integral formulations (2) and (5) from above can be interpreted as convolutions over the space which is one dimension higher than the boundary. Recall the two-dimensional convolution one obtains for the half-space and which can be solved rapidly with the FFT. It is 0 0 0 0 ( , ) ( , ) a ab b S u x y K x x y y dx dy   , (7) where u is the deformation,  is the load, a is the direction of the deformation and b is the direction of the load, (a,b)  {x, y, z}. Here, two coordinates x and y lie in the plane of surface S over which the integration is performed. This characteristic makes it possible to perfectly align a uniform grid on which the FFT is performed to calculate the convolution in (7) with the half space surface (see Fig. 1a). For arbitrary shapes, one has to enclose the entire surface with a uniform three-dimensional grid in order to apply the FFT (see Fig 1b). Integral formulations (2) and (5) clearly represent convolutions over this space. One only has to pay attention to appropriately set zeros at the grid points where there is no surface so as not to distort the results. In [15] an exemplary implementation of this technique is presented which illustrates the feasibility of the concept. Using the FFT on the three-dimensional space to obtain the integrals on the surface lowers the computational complexity to O(n 3 log n 1.5 ) from the complexity of O(n 4 ) which is needed to invert the classical BEM matrix. a) b) Fig. 1 A uniform grid aligned with the even surface of a half-space (a), and an arbitrary shape fully enclosed with a uniform three-dimensional grid (b), [15] 408 J. BENAD 3. IMPLEMENTATION DETAILS We will now highlight a few important details one has to consider when the technique described above is implemented. Surface normal vectors The components of the normal surface vectors in integral formulations (2) and (5) depend on x . This characteristic does not occur for the half-space because the surface is even and stretches to infinity. For the arbitrary case, however, the components have to be considered in detail. This can be illustrated with an example for the two-dimensional LAPLACE equation. With 2 2 0 0 ( ) ( )r x x y y    , (8) the fundamental solution (3) can be written as * 0 0 1 ( , ) ln( ) 2 u x x y y r      . (9) With Eq. (9) and q n u  , one obtains 0 0* 2 (( ) ( ) ) 2 x y n x x e y y e q r       . (10) Equations (9) and (10) can now be inserted into Eq. (2) which yields the boundary integral formulation 0 0 0 0 0 0 2 ( ) ( )1 1 ( , ) ( , ) ln( ) + . 2 2 x y S S x x n y y n c x y u x y r qdS u dS r         (11) Now we can turn our attention to the components of outward normal vector n . As stated above, they depend on the vector’s position on the border: nx = nx(x,y), ny = ny(x,y). In order to reveal the convolution terms, one has to split the integral formulation according to the normal vector components: 0 0 0 0 0 0 0 2 0 0 0 2 0 0 1 ( , ) ( , ) ln( ( , )) ( , ) 2 ( )1 ( , ) ( , ) 2 ( , ) ( )1 ( , ) ( , ) . 2 ( , ) S x S y S c x y u x y r x x y y q x y dS x x u x y n x y dS r x x y y y y u x y n x y dS r x x y y                   (12) (Note that in Eq. (12) we have also switched the variables in the kernel functions to highlight the convolution.) The procedure above has to be applied in the same style for the boundary element formulation of the Navier equation in order to reveal the convolution terms. Efficient Calculation of the BEM Integrals on Arbitrary Shapes with the FFT 409 Boundary shape and interpolation When the boundary of the arbitrary shape is not aligned with the uniform grid on which the FFT is performed one simply obtains the values on the boundary where it cuts through the grid (see Fig. 2). Subsequently, the data points on the grid can be interpolated to a certain desired spacing on the surface. Fig. 2 An arbitrary boundary which cuts through a uniform grid. The intersection points are marked with black dots Conjugate gradient method Note that although the FFT gives the results of the boundary integrals, the final results for the potential and the flow, or the deformations and the stresses which we seek still reside within these integral formulations (see Eqs. (2) and (5)). Therefore, an iterative method such as the conjugate gradient method has to be used to obtain them. Note that with the FFT technique used for the half space, only the surface stresses occur within the integral and require a conjugate gradient method, and not the surface deflections (see Eq. (7)). Further reduction of the computational complexity The kernel functions change very little far away from evaluation point 0 x . In close proximity to the point they change rapidly (see Fig. 3). Fig. 3 Kernel functions of the integral formulation (12) Therefore, one may obtain a further reduction in computational complexity when the uniform grid has a finer resolution close to the evaluation point than it has farther away from it. This technique is of particular interest if the values of only very few points of the 410 J. BENAD domain are of interest, which is the case for the method described here. In contrast to the FFT for the half-space where the results of the FFT are of interest on all discretization points in the domain, we seek here only the values which lie directly on or very near a certain boundary shape. With a finer uniform grid in close vicinity of the surface, one may save the effort of computing non-required values at many discretization points which lie within the arbitrary surface or which surround it. In order to realize this concept, one can perform a rough calculation on the entire uniform grid which fully encloses the surface but which has a low resolution. Subsequently, the influence of the rough discretization points close to the evaluation points is deleted and replaced with the influence of many discretization points from a finer local grid. Several of these finer local grids are placed along the surface but not in the areas where we have no interest in the results (see Fig. 4). Note that the convolutions over these finer grids can be obtained rapidly parallel to each other since their results do not depend on each other. Fig. 4 Fine local uniform grids placed along the surface of an arbitrary shape which is fully enclosed by a rough uniform grid 4. SIMULATION EXAMPLE In this section we will present a small example to illustrate how the FFT can be used to obtain the boundary integrals on arbitrary shapes. Contrary to the example given in [15] where a simple rectangular boundary was used, we will here extend the calculation method to more sophisticated curved shapes. We will show different implementation techniques which all have the same order of computational complexity. However, their accuracy varies significantly. Technique I: FFT, nearest-panel boundary Consider the boundary of an arbitrary two-dimensional shape as displayed in Fig. 5 and let it be fully enclosed by a uniform grid. Each panel through which the boundary cuts is highlighted gray in Fig. 5. The two intersection points of the boundary with the edges of a panel shall define the beginning and end of a constant boundary element. This leads to an irregular size distribution of the elements along the boundary. For this first preliminary study, we will accept this. Efficient Calculation of the BEM Integrals on Arbitrary Shapes with the FFT 411 Fig. 5 An arbitrary two-dimensional shape fully enclosed with a uniform grid. The intersection points of its boundary (gray thin line) with the grid (black dots) define start and end points of the boundary elements (bold black straight lines). In the left graph the grid has n = 10 panels on each side. The right graph displays a finer resolution with n = 56 We now seek to calculate the boundary integrals (12) of the LAPLACE equation (1) on the boundary of the chosen shape. In order to use the FFT for the operation, each active panel (marked gray in Fig. 5) is loaded with the potential and flow of the constant boundary element which cuts through it. We can then perform the convolutions with the kernel functions from Eq. (12) using the FFT in the same way as in [15]. In order to test the method, we use a simple analytical solution for the LAPLACE equation. Here, we use the exemplary solution a 1 / 2u x  (13) of Eq. (1). With the known geometry of the chosen shape and Eq. (13), we then have the exemplary boundary values for the potential and the flow. The boundary integrals (12) are now obtained with the FFT and the resulting numerical values ui of the potential are compared with the analytical solution ua on the boundary. The results are displayed in the first graph in Fig. 9. It becomes apparent that this first technique of loading only the nearest boundary panels and taking the results directly from the convolution delivers only poor results. The numerical values oscillate with a great error around the analytical solution. Also, it shall be mentioned that for a finer resolution the numerical values on the boundary obtained with this first method do not converge to the analytical solution. Technique II: FFT, weighted nearest-panel boundary In order to obtain more accurate results than with Technique I, we distribute the potential and flow of one boundary element not only to one panel but over three panels which are weighted so that their center corresponds precisely to the midpoint of the boundary element. The resulting active panels are shown in Fig. 6. To each main panel which is directly cut by the boundary (marked with a square), there is one horizontal 412 J. BENAD neighbor (marked with a cross) and one vertical neighbor (marked with a circle). Depending on the position of the panel midpoint we must choose either the left or the right neighbor for the additional horizontal panel, and either the upper or the lower neighbor for the vertical panel, so that by weighing the elements their center can lie directly at the midpoint of the boundary element. Fig. 6 The boundary elements (solid black lines) of the chosen shape shown together with their active panels (gray) from Technique II: To each main panel which is directly cut by the boundary (marked with a square), there is one horizontal neighbor (marked with a cross) and one vertical neighbor (marked with a circle). The left or the right neighbor for the additional horizontal panel and the upper or the lower neighbor for the additional vertical panel is determined depending on which one is closest to the midpoint of the boundary element. In the left graph the grid has n = 10 panels on each side. The right graph displays a finer resolution with n = 56 After the potential and the flow have been distributed to the active panels, we can proceed just as in Technique I. We perform the convolutions with the kernel functions from (12) with the FFT and resulting numerical values ui of the potential are compared with analytical solution ua on the boundary. The results can be seen in the second diagram in Fig. 8. There is a definite improvement when compared to the first technique; however, there is still a high oscillating error. Technique III: FFT, direct center integral, weighted nearest-panel boundary In order to further improve the results, we now calculate the boundary integrals (12) directly in a region close to a desired boundary value (see Fig. 7). Outside this small region, we continue using Technique II. Efficient Calculation of the BEM Integrals on Arbitrary Shapes with the FFT 413 Fig. 7 Exemplary midpoint of a boundary element (cross) around which a square is placed in which the integral is calculated directly. The grid has a resolution of n = 56 panels on each side and the small square has r = 7 panels in the horizontal and vertical direction Note that the direct integration in a small region around the boundary element does not increase the order of the computational complexity of the method. However, the accuracy is further improved by the technique as can be seen in the third graph in Fig. 8. Here, a small square with r = 7 panels on each side is used. It shall be mentioned that for a finer resolution of the grid the numerical values converge to the analytical solution. Therein, the number of panels r for the direct integration does not increase when the grid is refined. Technique IV: FFT, direct center integral, weighted nearest-panel boundary, interpolation The raw data from Technique III can still be improved. A large portion of the error which remains is due to the alternating boundary element sizes caused by the uniform grid which cuts through the arbitrary boundary. Thus, the data is linearly interpolated to a constant spacing along the boundary and each resulting value is averaged with its two neighbors to eliminate the influence of the alternations. This technique further improves the accuracy of the results as can be seen in the last graph in Fig. 8. The order of the computational complexity remains the same. The last point to make in the paper is that, in addition to the boundary values, also inner values of the problem are calculated with each method above. These inner values are not needed in order to obtain the boundary values but they are simply available after the boundary integral is obtained. This is due to the application of the FFT on the entire domain. At a certain distance to the boundary these values are highly accurate. Fig. 9a displays the inner values as they are obtained with the convolution in Technique I. Close to the boundary the error starts to increase rapidly. Fig. 9b displays the inner values as obtained in Technique III and Technique IV. Here, the graphs show only the portion of the results from the convolutions. The additional portion which leads to the numerical values for u on the boundary (compare Fig. 8) comes from the direct integration which is not shown in Fig. 9b. The last two graphs (see Fig. 9c) show again the inner values as obtained in the Technique III and Technique IV, only now the resolution n of the grid is increased ten 414 J. BENAD times. As the number of panels for the direct integration remains the same (r = 7) while the overall resolution is increased, the values which represent the potential within the shape now stretch all the way to its boundary. Fig. 8 Numerical results for u (black dots) displayed along boundary S of the shape and compared with the analytical solution for u on the boundary (red line). The graphs are displayed for a grid with n = 56 panels on each side. For the chosen shape this corresponds to N = 214 surface points. (See the red line in Fig. 9 for a view of the analytical solution for u displayed over the x-y plane.) Efficient Calculation of the BEM Integrals on Arbitrary Shapes with the FFT 415 Fig. 9 Exemplary results for the inner values available after the convolutions in a) Technique I, n = 56, r = 7, b) Technique III & IV, n = 56, r = 7, c) Technique III & IV, n = 560, r = 7. The red line shows the analytical solution for u on the boundary. (Note that the grid lines in c) are not displayed only for a better visual illustration.) 416 J. BENAD 5. CONCLUSIONS In this study we have further deepened the understanding of how the boundary integrals of the BEM can be obtained on arbitrary shapes with the FFT. We have built upon the results in [15] and presented an exemplary numerical calculation of the boundary integrals of the LAPLACE equation for a chosen curved shape. The accuracy of different implementation techniques is compared. We have also drawn attention to several general aspects of the concept: the importance of the normal surface vectors is highlighted, that is, the one with which one can split the boundary integral formulations on arbitrary shapes to obtain simple convolution integrals. We have further illustrated how the arbitrary boundary cuts through the uniform grid. Then it is highlighted that with the boundary integral formulations of the LAPLACE and Navier equation there is always a conjugant gradient method necessary to obtain the desired values on the boundary. Furthermore, a concept is mentioned to further reduce the computational complexity of the method by introducing a high-resolution grid close to the boundary and a low-resolution grid farther away. It can be concluded that the use of the FFT to obtain the boundary integrals of the BEM is appealing and worthy of further investigation. This is due to the lower computational complexity of the method O(n 3 log n 1.5 ) than the inversion of the full BEM matrix O(n 4 ), and to the potential for further decreasing this complexity. Furthermore, very efficient implementations are available for the FFT. Also, it should be mentioned that the method discussed in this work is very similar to the technique used with great success for the half space in contact mechanics. Thus, many proven techniques from this field could be adopted in the future. Among them is the rapid calculation of the contact area, the acquisition of stick and slip zones, and the inclusion of adhesion in the model [1]. Acknowledgements: The author would like to thank V. L. Popov for many valuable discussions on the topic and critical comments. REFERENCES 1. Popov, V.L., Pohrt, R., Li, Q., 2017, Strength of adhesive contacts: Influence of contact geometry and material gradients, Friction, 5(3), pp. 308-325. 2. Li, Q., Popov, V.L., 2018, Boundary element method for normal non-adhesive and adhesive contacts of power-law graded elastic materials, Computational Mechanics, 61(3), pp. 319-329. 3. Popov, V.L., 2017, Contact Mechanics and Friction, 2 ed, Springer, Berlin. 4. Pohrt, R., Li, Q., 2014, Complete Boundary Element Formulation for Normal and Tangential Contact Problems, Physical Mesomechanics, 17(4), pp. 334-340. 5. Ferguson, S., Bryant, J., Ganz, R., Ito, K., 2000, The influence of the acetabular labrum on hip joint cartilage consolidation: a poroelastic finite element model, Journal of biomechanics, 33(8), pp. 953-960. 6. Pakhaliuk, V., Poliakov, A., Kalinin, M., Pashkov, Y., Gadkov, P., 2016, Modifying and expanding the simulation of wear in the spherical joint with a polymeric component of the total hip prosthesis, Facta Universitatis-Series Mechanical Engineering, 14(3), pp. 301-312. 7. Bräunling, W., 2015, Flugzeugtriebwerke, 3. ed, Springer, Berlin. 8. Torenbeek, E., 1982, Synthesis of Subsonic Airplane Design, Kluwer Academic Publishers, Dordrecht. 9. Benad, J., 2019, Numerical methods for the simulation of deformations and stresses in turbine blade fir-tree connections, Facta Universitatis-Series Mechanical Engineering, accepted for publishing. 10. Gaul, L., Fiedler, C., 2013, Methode der Randelemente in Statik und Dynamik, 2 ed, Springer, Berlin. 11. Phillips, J., White, K., 1997, A precorrected FFT-method for electrostatic analysis of complicated 3-D structures, IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems, 16(10), pp. 1059-1072. 12. Masters, N., Ye, W., 2004, Fast BEM solution for coupled electrostatic and linear elastic problems, NSTI-Nanotech, 2, pp. 426-429. Efficient Calculation of the BEM Integrals on Arbitrary Shapes with the FFT 417 13. Lim, K., He, X., Lim, S., 2008, Fast Fourier transform on multipoles (FFTM) algorithm for Laplace equation with direct and indirect boundary element method, Computational Mechanics, 41, pp. 313-323. 14. Benedetti, I., Aliabadi, M., Davi, G., 2008, A fast 3D dual boundary element method based on hierarchical matrices, International journal of solids and structures, 45(7-8), pp. 2355-2376. 15. Benad, J., 2018, Acceleration of the Boundary Element Method for arbitrary shapes with the Fast Fourier Transformation, arXiv preprint arXiv:1809.00845.