Two accelerated isogeometric boundary element method formulations: fast multipole method and hierarchical matrices method

This work presents two fast isogeometric formulations of the Boundary Element Method (BEM) applied to heat conduction problems, one accelerated by Fast Multipole Method (FMM) and other by Hierarchical Matrices. The Fast Multipole Method uses complex variables and expansion of fundamental solutions into Laurant series, while the Hierarchical Matrices are created by low rank CUR approximations from the k−Means clustering technique for geometric sampling. Both use Non-Uniform Rational B-Splines (NURBS) as shape functions. To reduce computational cost and facilitate implementation, NURBS are decomposed into Bézier curves, making the isogeometric formulation very similar to the conventional BEM. A description of the hierarchical structure of the data and the implemented algorithms are presented. Validation is performed by comparing the results of the proposed formulations with those of the conventional BEM formulation. The computational cost of both formulations is analyzed showing the advantages of the proposed formulations for large scale problems.


INTRODUCTION
In general, numerical methods, such as Finite Element Method (FEM) and BEM, are based on the transformation of partial differential equations into integral equations and on the discretization of these integral equations through the creation of elements.There are some characteristics that are common to the numerical methods, among which we can mention the loss of precision of the results with the use of less refined meshes and the increase in processing time with the use of more refined meshes.These mentioned characteristics are partly due to the approximation of the geometry and field variables by the Lagrange polynomials.Lagrange polynomials are not able to accurately represent the geometry of most of the boundary of continuum mechanics problem domains such as circles, ellipses and hyperboles.Furthermore, there is no continuity of the derivatives of functions approximated by Lagrange polynomials between an element and its neighbors.With a view to solving problems like these, the idea was to use the same basis functions used in the Computer Aided Design (CAD) packages, called NURBS, to describe the geometry and to approximate the field variables.Thus arises the isogeometric analysis, which is widely discussed in the literature (Beer et al., 2019;Peigl and Tiller, 1996;Rogers, 2000;Hughes et al., 2005;Cottrell et al., 2009;Kagan and Fisher, 2000).With this new form of analysis, a segment in the parametric space corresponds to an isogeometric element located in the geometric space.Therefore, there is an elimination of the mesh generation step.In turn, the refinement is achieved without much additional effort (Li and Qian, 2011;Shene, 2011;Abramowitz and Stegun, 1972).
It is noteworthy that isogeometric analysis is more adaptable to BEM than to FEM specially for linear problems, as both the CAD system and BEM share surface definitions and do not need the domain discretization for linear problems.Even for some nonlinear problems, as contact mechanics, for example, the BEM doesn't require domain discretization.In the work of Loyola et al. (2022), an isogeometric boundary element formulation was successfully applied to contact mechanics problems.Another difficulty with the joining of CAD and FEM is the decrease in the sparsity of their matrices due to the high continuity of the NURBS basis functions (Collier et al., 2012).Note that the matrices of BEM are already full and this is not an extra problem.Thus, the isogeometric formulation of BEM brought encouraging results in terms of accuracy and efficiency.Some benefits can be observed in the literature (Beer et al., 2019), such as smoother geometries easily obtained, non-linear problems are solved without additional effort, among others.
Despite the known dimension reduction due to the fact that only the boundary is discretized, BEM has disadvantages due to the fact that its matrices are full and not symmetrical, imposing a high cost in the assembly and resolution of the linear system of equations.Several techniques were developed in order to face these problems, such as the Adaptive Cross Approximation (ACA) technique, which approximates the dense matrices of the BEM by a representation through hierarchical matrices.This representation is based on a binary tree structure that partitions the full matrix into smaller blocks, where each block will either be approximated by a low-rank matrix or the original block will be exactly used in the resolution of the linear system (Hackbusch, 2016).Other techniques were also developed, such as wavelets (Bucher et al., 2002), block methods (Rigby and Aliabadi, 1995;Crotty, 1982;Kane et al., 1990), agglutination processes and iterative techniques (Mansur et al., 1992;Barra et al., 1992).
Another well-known technique is the Fast Multipole Method (FMM) which has its roots in gravitational calculus of particle simulation models (Rokhlin, 1985;Greengard and Rokhlin, 1987).FMM improves the performance of the BEM due to the fact that the kernel of the fundamental solution can be expanded in series, which allows the separation of the relationship between the source point and the field point by inserting an intermediate point.There is also a decrease in the number of interactions due to the clustering of boundary elements in cells of a tree structure (Barnes and Hut, 1986).Due to its performance wherever it is applied, FMM is considered one of the top ten algorithms of the 20th century (Dongarra and Sullivan, 2000;Liu, 2009).There are other FMM methods developed for purely numerical kernels or for kernel-independent methods where there is no expansion of analytical kernels.Instead, they use an interpolation method like, for example, Chebyshev.These methods also make use of a matrix compression methods like Singular Value Decomposition (SVD) or Fast Fourier Transform (FFT) (Ying et al., 2004;Fong and Darve, 2009).
Techniques such as hierarchical matrices and FMM make use of iterative methods to solve the linear system of equations because in these techniques matrices of BEM cannot be calculated explicitly in order to save memory and allow large-scale simulations (Greenbaum, 1997).Some well-known iterative methods are the Gauss-Jacobi, Gauss-Seidel, minimum residue (MIN-RES), conjugate gradient, bi-stabilized conjugate gradient (BICGSTAB) and Generalized Minimal Residual Method (GMRES).These methods adapt very well to low rank approximation techniques such as those mentioned at the beginning of this paragraph.
Regarding the accelerated isogeometric formulation by the FMM, there are some related works.For example, in Matsumoto and Takahashi (2012), the FMM is coupled to an isogeometric formulation of the BEM, where B-splines of degree 2 and 3 are used explicitly as basis functions, applied to an infinite domain problem with Neumann boundary conditions and smooth boundary.In Wang et al. (2019), an element-free method (meshless method) coupled with isogeometric analysis and accelerated by the FMM is applied to potential problems.In Simpson and Liu (2016), a FMM coupled to Isogeometric Boundary Element Method (IGABEM) is presented for 3D problems, which uses Chebyshev Latin American Journal of Solids and Structures, 2022, 19(7), e463 3/26 interpolation and the M2L operators are approximated by SVD.In the literature, no other article was found where the FMM, with an analytical kernel, was used together with the Bézier extraction operator for the isogeometric analysis in BEM.
For the accelerated formulation by the hierarchical matrices method, some related works are: Ozdemir and Lee (2004) presents an algorithm IE-QR that builds an approximation QR of low rank using the modified Gram-Schmidt algorithm with cost with m and n being the admissible matrix block dimensions.In Kapur and Long (1998), an algorithm IES3 is presented, which consists of a kernel-independent method for electromagnetic simulations and costs ( log ) O N N ⋅ . An interpolative decomposition method based on QR rankrevealing factorization and costs ( ) where k is the rank of the matrix block (see (Boutsidis et al., 2009;Gu and Eisenstat, 1996;Voronin and Martinsson, 2017)).In Campos et al. (2017), a fast isogeometric formulation of the BEM is adapted to the hierarchical matrices method using ACA for low rank approximation.In Ayala et al. (2020), an approximation method CUR is presented for matrices corresponding to the admissible blocks, through a geometric sampling technique.To the best of the author's knowledge, there is no other article in the literature where the hierarchical matrix method, coupled with IGABEM with Bézier extraction, was subjected to a low-rank matrix compression method using the k -Means clustering technique for geometric sampling.
This paper presents two fast isogeometric formulations of the BEM, one using the FMM (Fast Multipole Accelerated Isogeometric BEM (IGAFMBEM)) and the other using hierarchical matrices with specific matrix compression.Both formulations use NURBS as shape functions and the Bézier decomposition method which, in turn, brings speed in the generation of NURBS and simplicity in the computational implementation.The boundary conditions are generic and are treated as Cabral et al. (1990 and1991) and Campos et al. (2017).

ISOGEOMETRIC FORMULATION
The parametric curves of the isogeometric analysis are related to each other as can be seen in Figure 1.All of them have a common core, which are the Bézier curves, and from then on they specialize according to the properties of each one.

Bézier curves
The i −th Bernstein polynomial base function, , i p B , of degree p is defined by the expression: with . A Bézier curve of degree p can be written as the linear combination of 1 p + Bernstein polynomial basis functions, dependent on the real parameter ξ , with 0 1 ξ ≤ ≤ , like: where P i are control points that form its control polygon.

NURBS
B-spline curves are a generalized form of Bézier curves.It is composed of one or more Bézier curves or segments with a commitment to continuity between the curves.Each control point influences only some segments of the B-spline curve and thus has local control of the curve.A base of univariate B-spline functions is defined from a set of knots or parametric domain U , which is a set of non-descending parametric coordinates written as { } 0 1 2 U= , , , , n p ξ ξ ξ ξ +  where i ξ ∈  is the i −th parametric knot, p is the polynomial degree of the B-splines basis functions and n is the number of basis functions.The B-splines base functions of degree p , , i p N , are defined recursively on U through the Cox-de Boor recursive formula of significant computational cost (Peigl and Tiller, 1996).
A B-spline curve of degree p in 2  is defined by where P i are the control points.
A NURBS curve is also defined from the set of knot U and a set of control points P i , of the form where , i p R are the NURBS basis functions and are defined as and is the weight function with i w being the weight corresponding to the i -th B-spline base function, , i p N , or associated with the i -th control point.

Bézier decomposition
New knots can be inserted in the set of knots, U , without changing the parametric or geometric properties of the curve, provided that for each knot inserted a new control point is added.Thus, inserting a new knot k p > , in the knot set, requires that 2 n + new B-spline basis functions are defined with the new set of knot . So that there is no change in the continuity of the curve, these , are formed from the original control points and are given by: Latin American Journal of Solids and Structures, 2022, 19(7), e463 5/26 where (see (Peigl and Tiller, 1996;Borden et al., 2011)).
The Bézier decomposition is obtained by inserting repeated knot together with all knots inside the knot set, U , until they have a multiplicity equal to the degree of the curve.After the insertions, the NURBS curve is decomposed into a set of Bézier curves, where each curve corresponds to a knot span of U .Considering a curve with n control points and calling j ξ the r we need to perform a decomposition, we can define j α according to Eq. ( 8).One can then write a matrix that relates the new control points to the old ones: Writing 1 P P = , we can rewrite Eq. ( 7) in matrix form in order to represent the sequence of control points created by the h −refinement, ( ) By repeating this operation r times, the final form of the decomposition is obtained: where , P is the original set of control points and P b is the final set, which can be called the Bézier control points.Remember that P has dimension n d be the set of Bernstein basis functions defined by the final set of knots.As the insertion of extra knots does not cause any geometric or parametric changes, the Bézier curve described by the new control points must coincide with the initial B-spline curve.So, from Eq. we have that: Since P is arbitrary: Latin American Journal of Solids and Structures, 2022, 19(7), e463  6/26 where C , called the Bézier extraction operator, relates the B-spline basis functions with the Bézier basis functions and, in this way, the NURBS curve can be generated directly from the Bézier base.

BOUNDARY INTEGRAL EQUATIONS
The proposed formulations are in line with the concept of the collocational boundary element method, but uses NURBS as shape functions for both geometry and numerical analysis.Let 2 Ω ∈  be a bounded domain with the boundary : Γ = ∂Ω , which may contain non-smoothness points.We propose to solve the following mixed boundary value problem for the Laplace equation: where n represents the outer normal to Γ , and u is a function that satisfies Laplace's equation in the domain Ω .After an algebraic reduction process, this problem is solved using the following Boundary Integral Equation (BIE): where ( ) c x is a jump term that arises from the integration process of the integral equation and depends on the geometry at the source point x ∈ Γ .The second integral on the right side of Eq. ( 15) has its value given according to Cauchy's definition of principal value.Furthermore, ( ) For discretization purposes, the continuous fields u and u n

∂ ∂
are written for each Bézier element, whose boundary conditions are interpolated by the corresponding collocation points.In the isogeometric method, the variables are approximated using the same geometry functions, in this case, the NURBS and from Eq. ( 4) we can rewrite Eq. ( 15) in the form: where c i u and


. It is worth noting that in Eq. ( 16) there is no approximation in geometry, but only in the variables u and u n ∂ ∂ . In turn, with the boundary being parameterized by t , we can rewrite Eq. ( 16) as: Latin American Journal of Solids and Structures, 2022, 19(7), e463  8/26 ( ) ( ) ( ) ( ) where the transformation to a standard domain was performed, hence

Γ
is the Jacobian of the transformation from physical coordinate space to parametric space.
It is still necessary one more change of variable in Eq. ( 17), in order to standardize the integration intervals for using the Gauss-Legendre quadrature, obtaining: where dt dξ is the Jacobian of the transformation from parametric space to local coordinate.The Bézier decomposition operation together with the Bézier extraction operator promotes a structure of boundary elements very similar to the elements in the conventional BEM.Each knot span in the parametric space corresponds to an independent isogeometric element of the adjacent isogeometric elements in the physical coordinate space, which correspond to the adjacent parametric knot spans to the initial one.

Numerical calculation of integrals
A very important part of the proposed formulations here, which is responsible for the convergence of the iterative process, refers to the accurate calculation of the integrals in Eq. ( 15) when the source point is located in the element over the which integration takes place.The first integral on the right-hand side has a weak singularity kernel, because of the fundamental solution of the potential, while the second integral has a strong singularity kernel, because of the fundamental solution of the potential derivative, making both integrals singular.It is also important to consider the occurrence of quasi-singular integrals, where the source point is very close to the element where the integration.
To treat singular, quasi-singular, and regular integrals, a transformation based on the hyperbolic sine function was used, which takes into account the position of the orthogonal projection of the source point on the element in which the integration occurs, as well as the distance from the source point to the element (Johnston and Elliott, 2004).The hyperbolic sine transformation has the ability to distinguish and satisfactorily handle each case as will be seen in the numerical results.For regular integrals, the transformation coincides with the Gauss-Legendre quadrature.

Definition of collocation points
The most appropriate way to distribute the collocation points is through the definition of Greville abscissas, because it is on which the control points have maximum influence (Simpson et al., 2012;Farin, 1996;Scott et al., 2013).They are defined as the average of p parametric coordinates: where p is the degree of NURBS and they are in quantity equal to the control points.
Latin American Journal of Solids and Structures, 2022, 19(7), e463 9/26 In the case of problems with boundary presenting corners, the Greville abscissas as defined in Eq. ( 19) have the inconvenience of locating collocation points at these corners.One way to avoid this inappropriate occurrence is to change the positions of the first and last abscissa as follows: ) where s is a coefficient that defines the displacement of the two extreme abscissas and, consequently, of the two corresponding collocation points.Here the value used was 0.5 s = , as recommended by the results presented by Wang and Benson (2015).
Isogeometric formulations have an important particularity, which is the fact that the boundary integral equation given in Eq. ( 15) is written in terms of the control points.In Figure 2a it can be seen, for an example, that the control points may not be on the boundary of the problem and, therefore, they do not constitute a suitable place for the installation of the collocation points.

Treatment of boundary conditions
The variation of the potential and its derivative can also be represented by NURBS, thus gaining smoothness between adjacent isogeometric elements.For this, a concept analogous to that of control points for the representation of geometry is used for the representation of the field variable (Cabral et al., 1990;Campos et al., 1997).
Since the values of the boundary conditions are given at the collocation points that are on the boundary of the problem, then this conditions cannot be directly enforced on Eq. ( 15) still.To solve this impasse, we apply a transformation matrix, , constituted by the NURBS shape functions to relate the values of the boundary conditions at the collocation points with the values at the control points: where the subindex indicates the values at the control points.By inverting the systems of equations given in Eq. ( 21), we can now solve Eq. ( 15).The solution is then obtained at control points and later through the system given in Eq. ( 21), the solution is obtained at the boundary.This procedure is carried out in each patch NURBS that composes the geometric boundary of the problem, each one with its type of boundary condition.

APPLICATION OF FMM TO ISOGEOMETRIC BEM
The Bézier decomposition operation and the Bézier extraction operator provided a simplification such that the FMM coupled to the isogeometric BEM practically does not differ from the FMM coupled to the conventional BEM using traditional elements (Liu, 2009).
The method starts with the construction of the quaternary tree from a square covering the entire boundary of the problem and from then on, each square is subdivided into four subsquares until each square has a number of Bézier elements, C Γ , limited by a prefixed value and each subdivision corresponds to a level change, starting from zero where the square covers the entire original boundary.Each square is known as a cell and a leaf cell is a square that will not be further subdivided.For an illustrative example, Figure 4 shows the transition of operations from level 3 to level 2 and vice versa.A C Γ element will belong to a cell if its center is inside that cell.In Figure 4a we have the upward process, step where the expansions of the multipole moments of the element C Γ around the center of the cell of the lowest level are calculated, which is transferred by translation M2M to the center of the cell level 2 .The downward process starts with the M2L translation transferring the moments accumulated in the center of the cell from the level 2 , C z ′ , to the local expansion point L z , from there to L z ′ by the L2L translation and finally the contributions of the element C Γ , from Latin American Journal of Solids and Structures, 2022, 19(7), e463 10/26 a distant cell, will constitute the actions of the source points 0 z , 1 z , and 2 z , see Figure 4b.The acronyms u FC and u LC ( u meaning upward) stand for parent cell and child cell, respectively.In order to perform the matrix-vector product from Eq. ( 15), basically two expansions are performed, each one referring to the integrals on its right side.To define the algebraic scenario, the following notations and associations must be established: source point z y iy = + ∈ , where i is the imaginary unit, see Figure 5.With the properties of complex variable differential and integral calculus in mind one can write: Latin American Journal of Solids and Structures, 2022, 19(7), e463 11/26 where ( ) ( ) is the fundamental solution of the potential in complex notation and

{ }
Re + represents the real part of the argument + .As well as: where is the fundamental solution of the potential derivative also in complex notation and is the outer normal vector at point z of Γ .
Therefore, the following equivalences are obtained: and where ( ) and ( ) u z come from the interpolation of the values of the boundary conditions corresponding to the element C Γ with the NURBS basis functions from Eq. ( 4), as well as 0 z and z are also given by the same equation.
One of the fundamental ideas that speed up the matrix-vector product is the separation of the relationship between the source point 0 z and the field point z .For this purpose, the kernel function ( ) 0 , u z z * given in Eq. ( 24) is expanded.

Introducing an expansion point C
z close to the point z such that the validity condition, see Figure 5, the following equation is obtained after an algebraic effort: Applying Taylor series to the second logarithm on the right side of Eq. ( 28), we arrive at: where the auxiliary functions ( ) , for 1 and ln

Multipole moments
The integral in Eq. ( 26) can now be calculated by the following expansion in multipole moments: are called moments of the C Γ element around the C z pole and independent of 0 z .
Analogously in Eq. ( 27), the calculation of the expansion at multipole moments is given as follows: where are the moments of C Γ around the C z pole, now referring to the kernel integral u n * ∂ ∂ .

Moment-to-Moment (M2M) translation
This operation is responsible for transporting the moment of C Γ around the center of the lowest level cell until it reaches the center of the cell at the level 2 of the tree.Figure 4a illustrates the translation of the calculated moment in C z , center of the u LC cell at the level 3, to C z ′ , center of the u FC parent cell at level 2, and the calculation is done as follows: ( ) ( ) ( ) Latin American Journal of Solids and Structures, 2022, 19(7), e463 13/26

Moment-to-Local (M2L) translation
The multipole moment of C Γ is carried from the center, C z ′ , from u FC cell, at level 2, to the center L z of d FC cell, also at level 2, by moment-to-local (M2L) translation, where u FC is found in the interaction list of d FC , through the local expansion coefficient given below:

Local-to-Local (L2L) translation
The calculated local expansion coefficient at the center L z of the d FC cell, at the level 2, needs to be translated to the center L z ′ of the d LC child cell, at the level 3, considering the example in Figure 4b.In the general case, the translation must happen to the lowest level, so that the contribution can be attributed to the actions of the source points of the leaf cell.The translation responsible for this step is the L2L given below:

Contribution of distant elements
Finally, the integrals referring to Eq. ( 26) and Eq. ( 27) when C Γ is far from the source points are calculated by FMM as follows: All M2M, M2L and L2L translations can be written using k M  in the same way as done with k M .
FMM is used to calculate the contributions of elements from cells far from the cell holding the source point.The calculation of the contributions of the elements of neighboring cells is performed with the conventional integration of BEM.The number of levels in the quaternary tree favors the validity conditions of the upward and downward processes, expressed in the legend of Figure 4, and contributes to the precision and speed of FMM as can be inferred from the results.

HIERARQUICAL MATRICES METHOD WITH CUR APPROXIMATION
BEM, as is known, produces full and non-symmetric matrices, but due to the characteristic of asymptotic smoothness of the fundamental solutions of the potential and the potential derivative, these matrices have the property that part of their submatrices can be represented by low-rank matrices.These submatrices are the final product of the hierarchical matrices method applied to BEM.The hierarchical form divides the matrices into blocks of sub-matrices and it all starts with structuring the nodes and boundary elements of the problem in the form of two binary trees, one tree for nodes and another for the boundary elements represented by their central points (Hackbusch, 2016;Bebendorf, 2000 and2008;Börm, 2010).From then on, each submatrix block will correspond to a group, which in turn is made up of two subgroups, one of source points and another of boundary elements being central points of the elements.The construction of these groups is done by the hierarchical clustering process from the two binary trees and is governed by the criterion called geometric admissibility condition, which takes into account the size and distance between the subgroups mentioned above and results in the largest possible blocks of low-rank submatrices, as follows: Latin American Journal of Solids and Structures, 2022, 19(7), e463 14/26 If a group is considered admissible, then its submatrix or block will be represented by a low-rank matrix providing the memory saving, otherwise, the submatrix will not be approximated and in this case there will be no savings because the representation will be full.The hierarchical matrices method applied to the matrices of BEM and coupled with a lowrank representation method seeks to obtain the largest possible number of admissible blocks or the largest possible blocks and this makes large-scale problems manageable both in terms of storage costs and in time spent on floating point arithmetic operations, as it can present linear or logarithmic-linear complexity depending on the implementation (Bebendorf and Rjasanow, 2003).
∈  be a matrix block of a BEM matrix corresponding to an admissible group, with m and n being the amounts of source points and boundary elements, respectively.A k −rank CUR approximation (or skeleton) is defined as where, in Julia language notation, [ ] J= , , , k j j j  sets of row and column indexes, respectively, with cardinality k , adaptively selected to ensure that

[ ]
A I,J is invertible and has the largest absolute value of possible determinant of all submatrices k k × of A (Mahoney   and Drineas, 2009; Kumar and Schneider, 2017).There are different methods of selecting J such as the Nearest- Neighbors criterion (NN) and which has been evaluated in high-dimensional problems showing good accuracy (March and Biros, 2017).Another method is Gravity Centers Sampling (GCS) with accuracy comparable to the ACA method (Ayala et al., 2020).The approach of this work is to select J through the k −Means clustering technique (MacQueen,  1967), which presents good performance as will be seen in the results.A brief outline of the k -Means clustering technique is presented in appendix A.
To select I specifically, the k −Means clustering partitions the subgroup of boundary elements of the admissible group into k clusters and then selects the closest k 's boundary elements to the centroids of the k clusters, where k n << .The indices of the selected k 's boundary elements will make up the set J .Then the pivoted partial QR decomposition is used to find I and, in short, proceeds as follows: calculate the factors 1 p , 1 Q and 1 R through the QR factorization of the matrix A[:,J] and through the permutation vector 2 p of the factorization QR of the factor T 1 Q , one finds the indices of the set I referring to the source points of the subgroup of the source points of the admissible group, which in turn correspond to the largest singular values of the A matrix, that is, the main source points (Golub and Loan,  2013).From now on, BEM builds a submatrix with the most significant rows and columns of the A matrix, that is, the

[ ]
A I,J submatrix.There are also other ways to obtain I in linear time, such as the Strong RRQR routines (Gu and Eisenstat, 1996) or maxvol (Goreinov et al., 2008).So getting a good set of column indexes J is critical.The computational cost of CUR approximation with k −Means is given as follows: ( ) O n k t ⋅ ⋅ floating point operations to obtain J , with t being the number of iterations of the k −Means algorithm (Han et al., 2012);   et al., 2020).A summary of the definition of the QR decomposition can be seen in appendix B.
The fact that the kernel of the integrals of the matrix entries H is the fundamental solution of the potential derivative, which in turn is asymptotically smoother than the fundamental solution of the potential, makes its submatrices have the tendency to have rank lower than the submatrices of the G matrix.So it makes sense to base the partial factorization QR pivoted from the submatrices of the H matrix.The qualitative information about the factorizations made on the submatrices of the H matrix, are used to obtain the submatrices of low rank of the G matrix.
Consequently, there is an important reduction in the number of rows and columns in relation to the original submatrix.
The implementation of the isogeometric formulation of BEM coupled with the hierarchical matrices method, proposed in this work, also experienced a simplification due to the use of the Bézier extraction operator derived from the Bézier decomposition operation.The integration did not undergo any changes and was carried out as is done in conventional BEM.Singular, quasi-singular and regular integrals were satisfactorily treated with the hyperbolic sine Latin American Journal of Solids and Structures, 2022, 19(7), e463 15/26 transformation (Johnston and Elliott, 2004).The linear system is also solved with the GMRES iterative method restarted and unpreconditioned (Saad and Schultz, 1986).
5.1 Discussion of the effectiveness of k −Means clustering in the CUR approximation The accuracy of the CUR approximation, of rank k , depends on the good conditioning of the matrix [ ] A I,J which, in turn, depends on adequate choices for the indices J of columns and I of rows, both with cardinality k .Such conditioning is directly related to the maximal volume of the simplex formed by the j −th column vectors, m j c ∈  , of the A matrix (Goreinov and Tyrtyshnikov, 2001).The volume of this simplex is given by the Cayley-Menger determinant as follows: ( ) ( ) where (Sommerville, 1958).
For a simple argument, suppose : u × →    , that is, the fundamental solution of the potential relating real domains and that it is used to calculate the matrix A entries.Without much algebraic effort, it can be deduced from the mean value theorem that ( ) where i x 's and j y ′s are source points and field points corresponding to subgroups of an admissible group, respectively, and lj ψ is a real number between l y and j y .It can be noted that the values of 2 jl d are related to the distance between the field points y 's selected by the set J .Therefore, if the field points are very close to each other, then a small value of k ν will result, which harms the conditioning of the [ ] A I,J matrix and decreases the performance of the CUR approximation.The k −Means clustering, due to its concept of intracluster similarities and intercluster differences, in terms of the distance metric, acts to keep these field points y 's as far apart as possible and, consequently, guarantying 2 jl d that are different from zero preserving the linear independence of the lines as much as possible, which favors the maximal volume of the simplex and the conditioning of the [ ] A I,J matrix.

NUMERICAL RESULTS
In the results that follow, the two proposed formulations are evaluated and compared with the conventional BEM using constant elements.All are applied to problems with known analytical solution, with the exception of the plate with many hole problem where only the processing time was analyzed.The simulation was performed on an Acer notebook with an Intel i5 processor at 2.3 GHz with Turbo Boost up to 2.8 GHz, 8 Gigabytes of RAM.In this work, the Julia programming language in version 1.5.3 was used.

Numerical setup
The number of degrees of freedom is given as a function of the h −refinement, which in turn inserts parametric knots into the U knot set.At each insertion of knots, an additional control point is created respecting Eq. ( 7).So, there is no change in geometry and smoothness between the Bézier curves.The degree p , or order 1 p + , of the NURBS curves will be explained in each case.The npg is the number of parametric Gauss-Legendre coordinates and these will be readapted by the hyperbolic sine transformation for the calculation of integrals.In this work, for the two proposed formulations, it was Latin American Journal of Solids and Structures, 2022, 19(7), e463 16/26 sufficient to set 8 npg = .The tolerance 6 10 ε − = was used as a weight for the stopping criterion for the preconditioned and restarted GMRES method (Saad and Schultz, 1986), which criterion is the defined relative residual where x represents a candidate solution of the linear system in the i −th iteration and 2 * denotes the 2 L norm for a finite dimensional space.
The instrument that will be used to evaluate the performance of the methods is the calculation of the approximation error by the 2 L norm, for the continuous case, between the numerical solution h u and the exact solution u on the boundary Γ = ∂Ω of the problem: 1. nexp is the number of terms of the expansion of moments in multipoles and ntylr is the number of terms in the local expansion.nexp and ntylr , for simplicity, they were considered equal to , which value can be estimated according to the formula ( ) 2. maxl is the maximum number of elements in a leaf cell: 30 maxl = . 3. levmx is the maximum number of levels in the quaternary tree: 15 levmx = .4. To calculate the multipole moments, 3 npg = .
6.1.2Parameters of the hierarchical matrices method: 1. nnucleos is the number of clusters produced by k -Means technique: 5 nnucleos = .

2.
_ max elem is the maximum number of elements in a leaf cell: _ 30 max elem = .

Heat transfer in a hollow cylinder
This is a heat transfer problem, by conduction, in a hollow cylinder modeled as a 2D problem as shown in Figure 6.This problem has an analytical solution that will be used to verify the results obtained by the proposed formulations in this work using NURBS of order 3.The results of conventional BEM using quadratic elements will also be compared.The known boundary conditions are temperature at the inner boundary, S i , and flux at the outer boundary, S e .The analytical solution for the temperature is given by: ( ) and for the flux: ( ) where i T and e q are the temperature and flux in the inner and outer boundary, respectively, and r , i r and e r are given in Figure 6.The data for the simulation are from Table 1.
Latin American Journal of Solids and Structures, 2022, 19(7), e463 17/26 In Figure 7, it is possible to observe the approximation error, by the norm 2 L , for the temperature, calculated along the boundaries S i and S e .It is observed that the three methods converge to the exact solution, however, the isogeometric formulations with fewer degrees of freedom reach greater accuracy in relation to the conventional BEM. Figure 8 presents the result for calculation of the approximation error for the heat flux and it is also possible to observe the convergence of the three methods for the exact solution.

Moulton problem
Consider the problem of heat transfer, by conduction, now in a rectangular plate AOBCD as shown in Figure 9, with constant thermal conductivity of value 1 . where u represents the temperature and q the heat flux in terms of the polar coordinates r and θ : Figure 11.For the two proposed formulations, NURBS of order , degree , were used and for the conventional BEM constant elements were used.It is worth noting that the boundary conditions imposed are not constant, but vary from point to point, including an important singularity at point O for the boundary conditions of heat flux.It should be noted that the boundary geometry is quite simple, and can be represented exactly by low order elements such as NURBS degree for the isogeometric and constant elements for the conventional BEM.
For the temperature result, it can be seen in Figure 10 that the three methods converge asymptotically to the analytical solution.Convergence is a little slower for the case of the heat flux given in Figure 11, which can be explained by the singularity at point O .In terms of accuracy and precision, here we also highlight the isogeometric formulations in relation to the conventional BEM.

Problem of heat transfer in a square plate with many holes
In order to complement the verification of the efficiency of the presented formulations, it is proposed to simulate, in order to access processing time, the problem of heat transfer in a square plate of 2 1 m area and constant thermal conductivity of , with an increasing number of evenly distributed circular holes.The size of the holes is modified so that the sum of the areas of the holes makes up the proportion of 12.47% of the plate area, which proportion is kept constant throughout the simulation.This problem has been analyzed in the reference (Liu, 2009), using Fast Multipole Accelerated BEM (FMBEM) with constant elements.Boundary conditions are applied with thermal insulation on the lower and upper outer boundary and temperatures T=0 and T=1 for the left and right outer boundary, respectively.Zero flux is also defined on the inner boundary.For an example of a plate with 16 holes see Figure 12.The simulation was performed by varying the number of holes between 1 and 10000.It was observed the need to maintain a uniformity in the length of the elements between the NURBS patches of the external and internal boundary.Otherwise, there would be two problems to deal with: first, there would be a need to use different amounts of Gauss-Legendre points on the outer and inner boundary for numerical integration, in order to prevent loss of precision; second, non-uniformity could weaken the validity conditions of FMM.The h −refinement was applied only on the external boundary to guarantee the uniformity mentioned above, while in each hole the amount of source points and boundary elements were kept constant.Thus, the number of degrees of freedom is given as a function of the number of holes and the refinement of the outer boundary.
The processing time was considered from the assembly of the matrices until the end of GMRES.The evaluation involved the conventional BEM, the FMBEM with constant elements (Liu, 2009) and the two proposed formulations in this article.Results can be seen in Figure 13

CONCLUSIONS
In this work two fast isogeometric formulations of the accelerated boundary element method were presented, one by the Fast Multiple Method and the other by the Hierarchical Matrices Method with low rank CUR approximation coupled to a geometric sampling method.The NURBS were used as shape functions and generated by the Bézier extraction operator obtained by the Bézier decomposition which, in turn, reduced the computational cost and brought simplicity in the implementation, making the isogeometric formulations very similar to the formulation of the conventional BEM regarding the structure of elements.
The results show that the proposed formulations have high accuracy and efficiency while dealing with generic geometries and boundary conditions.The effectiveness of the boundary conditions imposition method is also shown.
Regarding the processing time, the method coupled to FMM presented a complexity order close to ( ) O N as it is known in the literature.While the complexity order of operations of the isogeometric formulation of BEM with the Hierarchical Matrices Method was estimated close to ( log ) O N N ⋅ .It is worth mentioning that using 2  matrices, a better result is expected for the hierarchical matrices in terms of memory savings and computational complexity (Hackbusch et al., 2000;Hackbusch and Börm, 2002;Löhndorf, 2003;Börm and Hackbusch, 2004;Börm, 2006 and2013;Börm and Garcke, 2007;Hackbusch, 2016).The use of 2  matrices will be investigated in future work.
Therefore, the proposed formulations in this work show a promising horizon for scientific computing regarding mesh generation, precision, accuracy, storage cost and processing speed of large-scale problems.where D is the Euclidean distance.The S F function represents the sum of all distances between each element and the centroid of its cluster, but it can also be seen as a measure of clustering dispersion.As an ultimate goal, the algorithm seeks to minimize S F by finding a local solution.Given a set of points, the algorithm steps are as follows: 1. Initially, all points in the set P are distributed into random k clusters.2. Calculates through Eq. (A.1), the centroid of each cluster C i .
3. Assigns each point P j x ∈ to a cluster C i * , with centroid i y * closest to the point, that is, 4.After the previous step, many points will have changed groups, so it is necessary to update the centroids of each cluster, so the 2º step is repeated, where we will find a new centroid i y for the cluster C i .
5. The last two steps will be repeated iteratively, until the respective centroids no longer change or satisfy the established precision, then this iteration will be the local minimum.
6.The stopping test is based on the analysis of the differences between the centroids of the current iteration and the previous one, that is, , admits a decomposition A=QR , where Q is a matrix m k × with orthonormal columns and R is a k n × upper quasi-triangular matrix (Golub and Loan, 2013).For R to be upper triangular, the decomposition will require pivoting columns of the matrix A and this introduces a permutation matrix P , hence T AP=QR A=QRP ⇔ . The permutation matrix P is such that the diagonal elements of R are non- increasing, thereby improving precision and providing a basis for a more accurate knowledge of the numerical rank of A with lower computational cost than SVD .
To calculate a pivoted partial QR decomposition, successive orthogonalizations with pivoting over the columns of the A matrix one by one are performed.This task is terminated when the Frobenius norm of the remaining columns is less than a given computational tolerance ε .Let  be the smallest number of steps for which this tolerance is reached, hence the process results in partial factoring for the SVD decomposition (Halko et al., 2011).
fundamental solutions of the potential and the potential derivative, respectively, calculated on the boundary Γ .
of potential and its corresponding derivative at the i −th control point,

Figure 4 :
Figure 4: Validity condition: L L C i z z z z − << − and
to perform the factorization QR truncated over C ; 2 n k ⋅ complexity of operations to perform the factorization QR truncated over T Q .So the total cost is

2
the calculation of the error estimate of the expansion of moments(Liu, 2009).

Figure 9 :
Figure 9: Domain for the Moulton problem.
. It is worth noting the suitability of fast formulations for large problems.The changes in the slopes of the curves show the improvement in speed for the proposed formulations: ( ) O N for Fast Multipole Accelerated Isogeometric BEM (IGAFMBEM) and FMBEM; ( log ) the conventional BEM.The time difference between IGAFMBEM and FMBEM can be explained by the complexity of numerical operations performed.

Figure 13 :
Figure 13: Execution time between assembly and GMRES.
A=QR+E ,whereQ is an orthonormal m ×  matrix, R is an n × upper triangular matrix, and E is a residue matrix satisfying