An Advanced ACA/BEM for Solving 2D Large-Scale Elastic Problems with Multi-Connected Domains

T.GortsasS.V.TsinopoulosD.Polyzos

An Advanced ACA/BEM for Solving 2D Large-Scale Elastic Problems with Multi-Connected Domains

T.Gortsas1,S.V.Tsinopoulos2,D.Polyzos1,3

An advanced Boundary Element method(BEM)accelerated via Adaptive Cross Approximation(ACA)and Hierarchical Matrices(HM)techniques is presented for the solution of large-scale elastostatic problems with multi-connected domains like in fiber reinforced composite materials.Although the proposed ACA/BEM is demonstrated for two-dimensional(2D)problems,it is quite general and it can be used for 3D problems.Different forms of ACA technique are employed for exploring their efficiency when they combined with a BEM code.More precisely,the fully and partially pivoted ACA with and without recompression are utilized,while the solution of the final linear system of equations is accelerated via an iterative GMRES solver.The proposed methodology is demonstrated with the solution of large scale,plane strain elastic problems dealing with the bending of unidirectional fiber composite plates with large numbers of periodically or randomly distributed cylindrical elastic fibers embedded in a matrix medium.

Boundary Element Method;Hierarchical Matrices;Adaptive Cross Approximation;ACA;Large Scale Elastic Problems;Composite Materials.

1 Introduction

The Boundary Element Method(BEM)is a well-known and robust numerical tool successfully used for the solution of elastostatic and elastodynamic problems[Beskos(1997);Bonnet(1999);Aliabadi(2002);Beer,Smith,and Duenser(2008);Manolis and Polyzos(2009)].Two remarkable advantages it offers as compared to other numerical methods is the reduction of the dimensionality of the problem by one and its high solution accuracy.Despite the advantages the brutal application of BEM to large-scale problems suffers from very time consuming computations andhigh demands for computer memory capacity.Both problems come from the generation of the non-symmetric coefficient matrix and the solution of the final system of algebraic equations.More precisely,the fully populated matrices produced by BEM requireO(N2)operations for its buildup andO(N3)operations for the solution of the final matrix system through Gaussian elimination or typical LU-decomposition solvers.The use of iterative solvers decreases the operation requirements fromO(N3)toO(M×N2),with M being the number of iterations,but still remains inefficient for large scale problems.To the same conclusion we reach when parallel computing methods are exploited for the solution of the problem.In order to solve the aforementioned drawbacks of BEM and extend its application to real industrial problems,many techniques have been proposed so far in the literature.One can mention the Fast Multipole Boundary Element Method(FM/BEM)[Liu(2009)],the Fast Wavelet/BEM[Bucher,Wrobel,Mansur,and Magluta(2003);Ntalaperas,Tsinopoulos,and Polyzos(2010)],the Precorrected Fast Fourier Transform(FFT)AcceleratedBEM[Xiao,Ye,Cai,andZhang(2012)]and the Hierarchical Matrices-Adaptive Cross Approximation(ACA)/BEM[Rjasanow and Steinbach(2007)].

The Fast Wavelet/BEM relies on the fast wavelet transform so that to compress the dense and fully populated final matrices of BEM and change them to semi-banded and sparse ones.Although its use can be considered as a black box in a BEM code,appears the disadvantage of requiring the knowledge of the final system of algebraic equations the construction of which is in general computational expensive.Details one can find in the representative works of[Bucher,Wrobel,Mansur,and Magluta(2003);Ravnik,Škerget,and Zunic(2009);Ebrahimnejad and Attarnejad(2010);Wang,Ma,Meng,and Huang(2013)].

A very efficient method that accelerates drastically the solution process of a BEM code is that of the FM/BEM.The FM/BEM accelerates the computation of matrix[A]of the final system of algebraic equations[A]{x}={b}toO(N)operations and reduces the memory requirements toO(N)as well.This is accomplished, first,by analytically expanding the fundamental solutions of the differential operator of the problem around the centers of a hierarchical cell structure through fast multipole expansions,second,by performing analytical integrations utilizing constant elements and third,by using an iterative solver(e.g.GMRES)instead of a direct one.The main disadvantage of FM/BEM is its dependence on the fundamental solution of the problem and the use of constant elements in order to perform analytical integrations on flight without storage memory requirements.Representative works are those of[Bapat and Liu(2010);Frangi and Bonnet(2010);Wang,Miao,and Zheng(2010);Wang and Yao(2011);Yusa and Yoshimura(2013)],while the application of FM/BEM to the solution of composite material problems is demonstrated in the works of[Hu,Wang,Tan,Yao,and Yuan(2000);Kong,Yao,and Zheng(20002);Yao,Kong,Wang,and Wang(2004);Liu,Nishimura,Otani,Takahashi,Chen,and Munakata(2005);Lei,Yao,Wang,and Wang(2006);Wang and Yao(2008)].

The precorrected FFT/BEM is a method that accelerates less than the FM/BEM but utilizes the same idea of cells as the FM/BEM does.The idea is to project the boundary fields to a uniform grid of points,then to compute the grid fields by using FFT and finally to interpolate the grid fields to boundary nodes.The nearby interactions are evaluated directly as in conventional BEM.The result is a reduction of the total operations toO(NlogN).Representative works are those of[Phillips and White(1997);Xiao,Ye,Cai,and Zhang(2012)].

An alternative approach to accelerate BEM is the recently proposed adaptive cross approximation algorithm(ACA)along with hierarchical matrices[Bebendorf(2000);Bebendorf and Rjasanow(2003);Rjasanow and Steinbach(2007);Borm,Grasedyck,and Hackbusch(2003,2003a);Bebendorf(2008);Brancati,Aliabadi,and Benedetti(2009)].The acceleration here is achieved because only a small number of elements of the collocation matrix[A]are calculated,while the rest ones are approximated via the already calculated values.In a BEM collocation matrix,this is possible and effective due to the nature of the fundamental solutions,which are functions of the distance between the source and field points.According to ACA/BEM,the matrix[A]is organized into a hierarchical structure of blocks based on the problem’s geometry.Applying a geometrical criterion the blocks are characterized either as non-admissible,where the ACA algorithm is not ef ficient and thus,the conventional BEM is used or admissible,where ACA is effective and is used to approximate them by only calculating a small number of their rows and columns.Each admissible block is represented in a low rank matrix format via the product of two matrices formed by the previously calculated rows and columns,respectively.This low rank format,in conjunction with a usage of an iterative solver,results both signi ficant reductions in memory requirements and a CPU time due to the acceleration of the matrix vector multiplication.The ACA/BEM has already been successfully used for solving static and dynamic linear elastic problems,as it is explained in[Bebendorf and Grzhibovskis(2006);Benedetti,Alliabadi,and Davi(2008);Benedetti,Milazzo,and Alliabadi(2009);Zechner(2012)]for elastostatics and[Benedetti and Alliabadi(2009);Messner and Schanz(2010);Millazzo,Benedetti,and Alliabadi(2012)]for elastodynamics.

Although the FMM/BEM seems to be faster than ACA/BEM,it appears some disadvantages as it is compared to the ACA/BEM,namely,it requires the knowledge of the multipole expansions of the fundamental solution of the problem and signi ficant and complex modi fications in a conventional BEM code in order to be implemented.On the other hand,ACA/BEM is a black box algorithm applied upon the collocation matrix[A]and thus its implementation is the same regardless of the differential operator of the problem.The main structure of a conventional BEM code,including the significant integration library,remains the same and the only changes needed in the code are restricted to the assembly procedures.A comparison between FM/BEM and ACA/BEM can be found in the work of[Brunner,Junge,Rapp,Bebendorf,and Gaul(2010)].

In the present work,a two dimensional hierarchical ACA/BEM formulation is proposed for the solution of elastostatic problems dealing with multi-connected domains like those appearing in the perpendicular to the fibers cross-section of a fiber reinforced composite plate.Although many BEM formulations accelerated by hierarchical matrices and ACA methodologies have appeared in the literature,to our best knowledge,an ACA/BEM methodology for multi-connected domains appears for the first time.According to the proposed formulation the matricesandof each region,formed from the integrals of the elastostatic fundamental displacement and traction,respectively,are approximated by hierarchical matrices,constructed by means of the ACA algorithm.Thus,the global collocation matrixis not constructed explicitly,saving signi ficant amount of memory.The matrix approximation is accomplished by the implementation of four ACA formulations,the fully and partially pivoted ACA with and without recompression in order to study their efficiency.For the solution of the final linear system of equations an iterative solver(GMRES)is employed.The boundary conditions of the problem as well as the interface continuity conditions between the inclusions and the matrix are applied during the matrix vector multiplication procedure.From the approximated operatorsanda block diagonal preconditioner is constructed for the faster convergence of the solution.The efficiency of the proposed formulation is demonstrated by solving a large-scale elastostatic problem dealing with the bending of a fiber composite plate.For a given accuracy of approximation the CPU time as well as the memory demands are compared with the corresponding ones of the classical BEM.The largest problem solved corresponds to 1.2 million degrees of freedom utilizing a desktop pc with 64 GB RAM.

2 Conventional boundary element method

As it has been already mentioned,the boundary element methodology proposed in the present work concerns multi-connected domains and it is demonstrated through the problem described below.Consider a 2Drectangular fibrous composite plate of lengthLand widthD,as shown in Fig.1.The plate occupies a region Ω0of boundaryS0,is made of a matrix material with Young’s modulusEMand Poisson’s ratiovMand is reinforced withNfidentical circular fibers of radiusawith material propertiesEFandvF,respectively.The volume fraction of the fibers with respect to the plate isVf.The fibers are either randomly distributed(Fig.1a)or periodically arranged in a square pattern(Fig.1b).Each fiber occupies a region Ωiof boundarySi,wherei=1,Nf.The plate is fixed at its one end and is subjected to a bending loadPapplied at its free end.

Figure 1:2D rectangular fiber reinforced composite plates with(a)randomly distributed and(b)periodically arranged fibers.

In the present section,the just described 2D elastostatic problem is solved numerically using first the conventional boundary element method for various numbers of fibersNf.

The problem solution can be obtained by solving a combined system of boundary integral equations written for the matrix and each of theNffibers.The boundary integral equations for the matrix and for theithfiber are written as:

local geometry at point x(for a smooth boundarybeing the unity tensor)andare the 2D free space elastostatic fundamental displacement and traction,respectively,given as[Polyzos,Tsinopoulos and Beskos(1998)].

where r is the distance between the points x and y,njthe unit vector normal to the boundary at point y,r,jdenotes spatial derivatives of r and∂r/∂n is the directional derivative with respect to the normal vector at y.

According to conventional BEM formulation,the boundary S is discretized into E three-noded quadratic or two-noded linear isoparametric line boundary elements.For a node k,belonging to the boundary S the discretized integral equation(1)takes the form

where A(e)is the number of nodes of the element e(A=3 for a quadratic and A=2 for a linear element),Nastand for the element shape functions,Jethe Jacobian of the transformation from the global(X1,X2)to the local co-ordinate system ξ andandrepresent the nodal values of the displacement and traction vectors,respectively.Adopting a global numbering for the nodes,each pair(e,a)is associated to a number β and the equation(5)is written as

Collocating Eq(6)at all nodesL,the following linear system of algebraic equations is obtained

where the superscript M indicates matrix medium,are matrices formed by the integrals(7)and(8),respectively,with the term 1/2 added at the diagonal elements ofwhile the in dicesmandγtake values 0,1,2,...Nf,corresponding to the total number of nodesthat the boundariesS0,S1,...,SNfhave been discretized into,respectively.The vectorscontain the nodal displacement and traction vectors of the nodesLm,respectively.

Similarly,collocating the discretized version of the integral equation(2)for all the nodes of the fibers,the following system of algebraic equations can be obtained

Taking into account the continuity conditions at the matrix- fiber interfacesthe system of equations(10)and(11)contains all the unknown displacement and traction in ter facial nodal values,while the half nodal values of u0and t0are known through the boundary conditions at the boundaryS0and the other half ones unknown.

Whenβ/=k,integrals(7)and(8)are non-singular and can be easily computed numerically by Gauss quadrature.In case whereβ=k,the integrals(7)and(8)become singular with singularitiesO(lnr)andO(1/r),respectively.Singular integrals are evaluated with high accuracy via direct integration techniques[Frangi and Guiggiani(2000)].

Combing equations(10)and(11)and rearranging according to the boundary conditions at the boundaryS0,one obtains a linear system of algebraic equations of the form

where the vectors x and b contain all the unknown and known nodal components of the boundary fields,respectively,while the solution of system(12)is usually accomplished through a typical LU decomposition algorithm.

3 Hierarchical and ACA accelerated BEM

In conventional BEM,the matrixof the system(12)is full populated requiringO(N2)operations for its solution withNbeing the DOFs of the problem.Taking also into account that the condition number of the system becomes worse asNincreases,it is apparent that conventional BEM is not able to treat realistic problems with hundreds of thousands DOFs.In order to overcome that difficulty and solve the problem described in the previous section for a large number of fibers,a hierarchical ACA accelerated BEM is proposed.Furthermore,a significant reduction of the problem solution time is also accomplished with the aid of a GMRES iterative solver.

The departure point of the proposed method is that both matricesand,appearing in equations(10)and(11),are represented hierarchically using a block tree structure.By means of simple geometric considerations the blocks corresponding to large distancesrare characterized as far- field blocks(or admissible)and compressed through low rank matrices found by ACA algorithm,while the rest blocks of the tree,which are dominated by the singular behavior of the kernels(3)and(4),are characterized as near field blocks(or non-admissible)and are calculated explicitly as in conventional BEM.

The ACA/BEM methodology proposed here can be analyzed in steps as follows:Step 1:The nodes of each region Ωiare partitioned in a tree structure containing clusters that include only closely located nodes.Starting from a cluster containing all the nodes of the region,two sub-clusters are created.This procedure is repeated recursively for the created sub-clusters and is ended when the number of nodes inside a cluster is equal or less than a prede fined numberLt,which is named the cardinality of the tree.The nodes belonging to the same cluster are renumbered so that the corresponding integrals ofandmatrices to be placed at neighbor rows and columns.The splitting of the clusters can be accomplished with the aid of various bisection techniques like the principal component analysis[Bebendorf(2008)]adopted in the present work.According to that technique,the centroid of each already created cluster of nodes is calculated by using the following equation:

where,Lcis the number of nodes of the cluster anddenote the coordinates of the cluster nodes,andi=1,2.Then,the corresponding cluster’s data matrix of dimensionsd×Lcis formulated in the following way:

wheredis the number of spatial dimensions(d=2 for a 2D problem).

Applying the Singular Value Decomposition(SVD)algorithm for the matrix(14),one obtainswithbeing the matrices provided by SVD,˜U andare orthogonal matrices and˜S is a diagonal matrix.It must be noted that the matrixdoes not need to be calculated explicitly.We define as principal direction the vector represented by the first column of matrixand symbolized as w.Considering the centroid(13)as the base of any position vector,the nodes of the cluster are separated by considering the positive or negative value of their projection along the aforementioned principal vector w.IfCdenotes the parent cluster,the two sub-clusterC1andC2containingLc1andLc2nodes respectively are de fined as

The first four levels of the just described procedure for the region Ω0is depicted in Fig.2.

Figure 2:Clustering of initial domain Ω0to 2,4 and 8 clusters.

Figure 3:Boxes ΩRand ΩCthat enclose the elements belonging to the clusters CR andCC,respectively,used in the admissibility criterion of Eq.1.

The admissibility criterion adopted in the present work can be described as follows:consider a block of the entire matrix the rows and columns of which are represented by the clustersCcandCR,respectively.The elements associated with the nodes,contained inCRandCC,are enclosed in two boxes,ΩRand ΩC,respectively.Each box has its sides parallel to the Cartesian axes and is the smallest possible,with respect to the location of its nodes(Fig.3).The corners are numbered anti-clockwise with the 1stand 3ndcorner corresponding to the minimum and maximum coordinates of the box,respectively,as shown in Fig.3.The block is admissible when the following geometric condition is true:

wheredR,dCare the diagonals of the boxes ΩRand ΩC,respectively,dist(ΩR,ΩC)is the distance of the corresponding bounding boxes[Borm,Grasedyck,and Hack-busch(2003)]andais a positive parameter.In the present workahas been chosen to be 1.2.

Step 3:The admissible blocks of the tree are approximated by low rank matrices with respect to a prescribed accuracyε,by using the ACA algorithm.More precisely,considering an admissible block sub-matrixof matricesor,with dimensionsN×Land a full rankR=min{N,L},the blockcan be written as:

The memory requirements and matrix multiplication CPU cost of a low rank block areO(K(N+L)),while for the corresponding full rank representation areO(N·L).It is obvious that the low rank approximation is ef ficient when the conditionK(N+L)≪N·Lholds.

The best low rank approximation with respect to the Frobenius norm,for a given accuracyε,can be found by means of the SVD(Bebendorf,(2000)).More precisely,according to SVD the matrixreads

and ignoring the singular valuesσi,withi=K+1,...,R.

Although the best low rank approximation is provided by the SVD,its cubic CPU cost renders that technique prohibitive for real world applications.Thus,instead of the SVD the Adaptive Cross Approximation(ACA),which is a special low rank approximation technique,is employed.The main idea of ACA is to construct a representation of(Eq.18)by suitably choosing a small subset of the rows and columns of a matrix.Based on this idea two algorithms have been developed(Bebendorf(2000),Bebendorf and Rjasanow(2003));the ACA with full pivoting,which is anO(K2·N·L)algorithm and requires as starting point the calculation of the entire matrix,and the partially pivoted ACA,which is anO?K2(N+L)?algorithm and requires the calculation of only a small part of.Apparently the partially pivoted ACA is faster and consumes less memory than the fully pivoted one,but the approximation accuracyεis not guaranteed because its stopping criterion is heuristic since thein Eq.(19)cannot be calculated exactly.In the present work the above mentioned drawback is cancelled applying extra convergence checks according to the methodology proposed by Bebendorf(2008).

Both the fully and partially pivoted ACA algorithms,as they are implemented in the present work,are illustrated in the Appendix.

Step 4:Although the ACA provides an efficient approximation,it can be improved further by optimizing the compression of the admissible blocks through recompression of matrix sub-blocks.A number of techniques have been proposed to that purpose and significant reduction of storage memory can be accomplished(Grasedyck(2005)).In the present work the recompression is accomplished by means of the SVD.Each block is recompressed immediately after its approximation by ACA,thus reducing the memory requirements in the process of assembling of hierarchical matrices.

More specifically,in the SVD algorithm described in the previous step the matricesare further decomposed according toQRdecompositionsrespectively.Next,the SVD of the matrixis computed,wherematrices and the optimum rankkSVD(≤k)is found so that the condition(22)to be fulfilled for prescribed accuracyε.Finally,the optimum compressed matricesandN×Kare given by the relations

The non-admissible leaf blocks can also be compressed by using the standard SVD procedure explained in Step 3.

However,for non-admissible leaf blocks positioned exactly on the diagonal of the hierarchical matricesorthe factorization is not efficient.These blocks are full rank blocks and the dimensions ofand,provided by SVD,are of dimensionsN×RandL×R,respectively,which implies that storing the full matrixof dimensionsN×L,consumes less memory.

4 Numerical results

In order to demonstrate the efficiency of the proposed hierarchical ACA/BEM formulations,the 2D elastostatic bending problem,described in section 2,is solved using both ACA/BEM and conventional BEM.To compare similar things,a GMRES iterative solver is also employed for the solution of the system(11)in conventional BEM.Taking into account that fiber regions are uncoupled to each other,the matrix˜A of the final system of algebraic equations is never formed explicitly,thus saving significant amount of memory.The GMRES multiplications are performed directly on Equations(10)and(11),while a block left diagonal preconditioner is used to accelerate the convergence.The dimensions of each block in the preconditioner are chosen to be approximately equal to the number of nodes that each fiber is discretized into.The inversion of each block is accomplished by using the LU decomposition algorithm.The four ACA techniques described in the previous section and illustrated in Appendix,namely the fully and partially pivoted ACA with and without recompression,are used in the ACA/BEM that solves the aforementioned problem.The cardinality of the cluster treeLtand the compression matrix accuracyε,appeared in equations(19),(22)and used in the algorithms of partially and fully pivoted ACA explained in the appendix,are chosenLt=32 andε=10-5,respectively.The GMRES accuracy for solving the algebraic system of equations(12)is chosenεGMRES=10-6,noting that the Iterative Solvers Package(ITSOL)free source routines were used to this end.For the computation of SVD and QR decompositions the free source CLAPACK routines were used.

In table 1,the maximum deflections of the composite plate with fibers periodically arranged in a square pattern are listed for various DOFs N and number of fibersNf.The results were calculated using the above mentioned ACA/BEM formulations and the conventional BEM.The problem parameters,as de fined in section 2,are the following:L=9m,D=3m,EM=66GPa,vM=0.31,EF=360GPa andvF=0.25,uf=0.35 andP=30MN/m.

Table 1:Accuracy in maximum deflection calculation obtained by the partially and fully ACA with and without recompression formulations and the conventional BEM.

Figure 4:Total CPU time(a)and the memory demand(b),using conventional BEM and fully and partially pivoted ACA with and without recompression formulations.

Observing the results depicted in Table 1,one can say that the first general remark is that for a prede fined compression matrix accuracyε=10-5we obtain results with an error being less than to 10-4in the maximum deflection of the composite plate.Another remark is that the partially pivoted ACA without recompression is the most accurate,while the partially pivoted ACA with recompression is the most efficient.

In Figs 4(a)and 4(b)the normalized total CPU time and memory requirements as function of DOFsNare depicted,respectively.In the total CPU time,both the time required for the evaluation of the matrices,and the system solution time have been considered.The total CPU time is normalized by the corresponding time required for the solution of the problem for 100000 DOFs by means of the conventional BEM.Figs 4(a)and 4(b)reveal the well known statement that the total CPU time and the memory demand using conventional BEM is of orderO(N2).Among the four ACA formulations the partially pivoted ACA with recompression is the most efficient.The CPU time and memory requirements for solving a 2D elastic problem with 105DOFs using conventional BEM are about the same with the corresponding ones needed for the same problem with 106DOFs when the partially pivoted ACA is utilized,i.e.,an order of magnitude reduction is achieved.In the sequel,the displacement fields and the corresponding stresses are evaluated for a composite plate with different number of elastic fibers.Figs 8 and 9 depict the contours of displacementuxand normal stressσy,respectively,for the plates with periodically arranged and randomly distributed fibers withNf=48 andVf=0.35.In the same figures,the corresponding displacements and stresses for a composite plate homogenized according to Christensen and Lo(1979)self-consistent homogenization model are also provided.Similarly,in Figs 10 and 11 the corresponding contours forNf=1000 andVf=0.35,are depicted.As it is apparent from all contours,the Christensen model becomes accurate only for large number fibers and when they are randomly distributed.

Figure 5:Contours of the displacement magnitude,for Nf=48 and uf=0.35;(a)homogenized plate(b)periodically arranged fibers and(c)randomly distributed fibers.

Figure6:Contours of normal bending stress,forNf=48anduf=0.35;(a)homogenized plate(b)periodically arranged fibers and(c)randomly distributed fibers.

Figure 7:Contours of the displacement magnitude,for Nf=1083 and uf=0.35;(a)homogenized plate(b)periodically arranged fibers and(c)randomly distributed fibers.

Figure 8:Contours of normal bending stress,for Nf=1083 and uf=0.35;(a)homogenized plate(b)periodically arranged fibers and(c)randomly distributed fibers.

5 Conclusions

An advanced Boundary Element method(BEM)accelerated via Adaptive Cross Approximation(ACA)and Hierarchical Matrices(HM)techniques for the solution of large-scale elastostatic problems with multi-connected domains has been proposed.Fully and partially pivoted ACA techniques with and without recompression have been testing with the most efficient being the partially pivoted ACA with recompression.The solution of the final linear system of equations has been accelerated via an iterative GMRES solver and the total gain in storage memory and solution time renders the proposed ACA/BEM a significant candidate for solving realistic problems.The largest problem solved in the framework of the present work corresponds to 1.2 million degrees of freedom utilizing a desktop PC with 64 GB RAM.The proposed methodology has been demonstrated with the solution of large scale,plane strain elastic problems dealing with the bending of unidirectional fiber composite plates with large numbers of periodically or randomly distributed cylindrical elastic fibers embedded in a matrix medium.The obtained results have been compared to corresponding ones taken for a plate homogenized according to Christensen and Lo(1979)self-consistent homogenization model.An interesting conclusion revealed by the solution of the aforementioned problems is that Christensen’s model becomes accurate only for large number fibers and when they are randomly distributed.

Acknowledgement:This research has been co- financed by the European Union(European Social Fund-ESF)and Greek national funds through the Operational Program “Education and Lifelong Learning”of the National Strategic Framework(NSRF)-Research Funding Program:ARCHIMEDES III.

Appendix

In this appendix the two algorithms dealing with the implementation of the fully and partially pivoted ACA are illustrated.

Algorithm 1:Fully pivoted ACA

Let’s consider a matrix blockof dimensionsN×LwithN≤L.

Algorithm 2:Partially Pivoted ACA

Let’s consider a matrix blockof dimensionsN×LwithN≤L.

1)SetK=1 and form a listZcontaining the indices of the block rows,i.e.,Z={i1,i2,...,iN},noting that the block,as already mentioned,is of dimensionsN×L.The indicesiKare sorted in an ascending order.In the first two positions are placed the indices(one for each DOF)corresponding to the node of theCRcluster which is closer to the center of the ΩCbox and so on.By taking into to account the behavior of the kernels(3)and(4)with respect to the distancer,this sorting increases the probability that the maximum absolute value of the block is located at rowsi1ori2.

12)If the stopping criterion is not satisfied and Z/=∅setK=K+1 return and repeat the algorithm form Step 2.

13)If the stopping criterion is satisfied,the matricesandare formed,via Eqs(19),using all the above calculated vectors aKand bK.

It should be mentioned that(i)in the case ofN＞Lthe above described partially pivoted ACA algorithm is modified by considering that the listZcontains the indices of the block columns i.e.,Z={j1,j2,...,jN}and the followed algorithm is modified similarly.(ii)As it is already mentioned,the stopping criterion in the step 11 is heuristic,which implies that there are cases where although the algorithm terminates the calculated vectorsaKandbKdo not provide an approximantwith the given accuracyε.In order to overcome this drawback extra convergence checks according to the methodologies introduced in[Bebendorf(2008)]are applied.

Aliabadi,M.H.(2002):The boundary element method:applications in solids and structures,vol.2.John Wiley&Sons Ltd.,England.

Bapat,M.S.;Liu,Y.J.(2010):A new adaptive algorithm for the fast multipole boundary element method.CMES:Computer Modeling in Engineering&Sciences,vol.58,no.2,pp.161-184.

Bebendorf,M.(2000): Approximation of boundary element matrices.Numerische Mathematik,vol.86,pp.565-589.

Bebendorf,M.(2008):Hierarchical matrices.a means to efficiently solve elliptic boundary value problems.Berlin:Lecture Notes in Computational Science and Engineering,vol.63,Springer Berlin,Heidelberg.

Bebendorf,M.;Grzhibovskis,R.(2006):Accelerating Galerkin Bem for linear elasticity using adaptive cross approximation,Mathematical Methods in the Applied Sciences,vol.29,pp.1721-1747.

Bebendorf,M.;Rjasanow,S.(2003):Adaptive low-rank approximation of collocation matrices.Computing,vol.70,pp.1-24.

Beer,G.;Smith,I.;Duenser,C.(2008): The boundary element method with programming.Springer.

Benedetti,I.;Alliabadi,M.H.(2009):A fast hierarchical dual boundary element method for three dimensional elastodynamic crack problems.International Journal for Numerical Methods in Engineering,vol.84,no.9,pp.1356-1378.

Benedetti,I.;Alliabadi,M.H.;Davi,G.(2008):A fast 3D dual boundary element method based on hierarchical matrices.International Journal of Solids and Structures,vol.45,pp.2355-2376.

Benedetti,I.;Milazzo,A.;Alliabadi,M.H.(2009):A fast dual boundary element method for 3d anisotropic crack problems.International Journal for Numerical Methods in Engineering,vol.80,no.10,pp.1356-1378.

Beskos,D.E.(1997):Boundary element methods in dynamic analysis:Part II(1986-1996).Applied Mechanics Reviews,vol.50,pp.149-197.

Borm,S.;Grasedyck,L.;Hackbusch,W.(2003):Introduction to hierarchical matrices with applications.Eng Anal Boundary Elem,vol.27,pp.405-422.

Borm,S.;Grasedyck,L.;Hackbusch,W.(2003a):Hierarchical matrices,Lecture notes no.21,Max Planck Institute for Mathematics in the Sciences,Leipzig,Germany.

Bucher,H.F.;Wrobel,L.C.;Mansur,W.J.;Magluta,C.(2003):Fast solution of problems with multiple load cases by using wavelet-compressed boundary element matrices.Commun.Numer.Meth.Engng,vol.19,pp.387-399.

Bonnet,M.(1999):Boundary Integral Equation Method for Solids and Fluids.John Wiley and Sons.

Brancati,A.;Aliabadi,M.H.;Benedetti,I.(2009): Hierarchical adaptive cross approximation GMRES technique for solution of acoustic problems using the boundary element method.CMES:Computer Modeling in Engineering&Sciences,vol.43,no.2,pp.149-172.

Brunner,D.;Junge,M.;Rapp,P.;Bebendorf,M.;Gaul,L.(2010):Comparison of the fast multipole method with hierarchical matrices for the helmholtz-BEM.CMES:Computer Modeling in Engineering&Sciences,vol.58,no.2,pp.131-158.

Christensen,R.M.;Lo,K.H.(1979):Solutions for effective shear properties in three phase sphere and cylinder models.J Mech Phys Solids,vol.27,pp.315-330.

Ebrahimnejad,L.;Attarnejad,R.(2010):Fast solution of BEM systems for elasticity problems using wavelet transforms,vol.87,no.1,pp.77-93.

Frangi,A.;Bonnet,M.(2010):On the application of the fast multipole method to helmholtz-like problems with complex wavenumber.CMES:Computer Modeling in Engineering&Sciences,vol.58,no.3,pp.271-296.

Frangi,A.;Guiggiani,M.(2000):A direct approach for boundary integral equations with high-order singularities.International Journal for Numerical Methods in Engineering,vol.49,pp.871-898.

Golub,G.H.;Van Loan,C.F.(1996): Matrix Computations.John Hopkins University Press,London.

Grasedyck,L.(2005):Adaptive recompression of H-matrices for BEM,Computing,vol.74,pp.205-223.

Hu,N.;Wang,B.;Tan,G.W.;Yao,Z.H.;Yuan,W.F.(2000):Effective elastic properties of 2-D solids with circular holes:numerical simulations.Compos.Sci.Technol.,vol.60,no.9,pp.1811-1823.

Kong,F.Z.;Yao,Z.H.;Zheng,X.P.(2002):BEM for simulation of a 2D elastic body with randomly distributed circular inclusions.Acta Mech.Solida Sinica,vol.15,no.1,pp.81-88.

Lei,T.;Yao,Z.H.;Wang,H.T.;Wang,P.B.(2006):A parallel fast multipole BEM and its applications to large-scale analysis of3-D fiber-reinforcedcomposites.Acta Mech,Sinica,vol.22,no.3,pp.225-232.

Liu,Y.(2009):Fast multipole boundary element method.Campridge University Press,UK.

Liu,Y.J.;Nishimura,N.;Otani,Y.;Takahashi,T.;Chen,X.L.;Munakata,H.(2005):A fast boundary element method for the analysis of fiber-reinforced composites based on a rigid-inclusion model.Journal of Applied Mechanics,vol.72,pp.115-128.

Manolis,G.D.;Polyzos,D.(2009):Recent advances in boundary element methods.A Volume to Honor Professor Dimitri Beskos.Springer.

Messner,M.;Schanz,M.(2010):An accelerated symmetric time-domain boundary element formulation for elasticity.Engineering Analysis with Boundary Elements,vol.34,pp.944-955.

Millazzo,A.;Benedetti,I.;Alliabadi,M.H.(2012): Hierarchical fast BEM for anisotropic time-harmonic 3-D elastodynamics,Computers and Structures,vol.96-97,pp.9-24.

Ntalaperas,D.;Tsinopoulos,S.V.;Polyzos,D.(2010):A Fast Wavelet/BEM for Wave Scattering Solutions.In Mathematical Methods in Scattering Theory and Biomedical Engineering,A.Charalambopoulos,D.Fotiadis,D.Polyzos,Eds,World Scientific Publishing.

Phillips,J.R.and White,J.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,vol.16,no.10,pp.1059-1072.

Polyzos,D.;Tsinopoulos,S.V.;Beskos,D.E.(1998):Static and dynamic boundary element analysis in incompressible linear elasticity.European Journal of Mechanics,A/Solids,vol.17,pp.515-536.

Ravnik,J.;Škerget,L.;Zunic,Z.(2009): Comparison between wavelet and fast multipole data sparse approximations for Poisson and kinematics boundary-domain integral equations.Comput.Methods Appl.Mech.Engng,vol.198,pp.1473-1485.

Rjasanow,S.;Steinbach,O.(2007): The Fast Solution of Boundary Integral Equations.Springer New York.

Wang,L.;Ma,Y.;Meng,Z.;Huang,J.(2013):Wavelet operational matrix method for solving fractional integral and differential equations of Bratu-type.CMES:Computer Modeling in Engineering&Sciences,vol.92,no.4,pp.353-368.

Wang,Q.;Miao,Y.;Zheng,J.(2010):The hybrid boundary node method accelerated by fast multipole expansion technique for 3D elasticity.CMES:Computer Modeling in Engineering&Sciences,vol.70,no.2,pp.123-152.

Wang,H.T.;Yao,Z.H.(2008):A rigid- fiber-based boundary element model for strength simulation of carbon nanotube reinforced composites.Comput.Model.Eng.Sci.,vol.29,no.1,pp.1-13.

Wang,H.T.;Yao,Z.H.(2011):A fast multipole dual boundary element method for the three-dimensional crack problems.CMES:Computer Modeling in Engineering&Sciences,vol.72,no.2,pp.115-148.

Xiao,J.;Ye,W.;Cai,Y.;Zhang,J.(2012):Precorrected FFT accelerated BEM for large-scale transient elastodynamic analysis using frequency-domain approach.International Journal for Numerical Methods in Engineering,vol.90,no.1,pp.116-134.

Yao,Z.H.;Kong,F.;Wang,H.;Wang,P.(2004):2D simulation of composite materials using BEM.Engineering Analysis with Boundary Elements,vol.28,pp.927-935.

Yusa,Y.;Yoshimura,S.(2013): A fast regularized boundary integral method for practical acoustic problems.CMES:Computer Modeling in Engineering&Sciences,vol.91,no.6,pp.463-484.

Zechner,J.(2012): A fast boundary element method with hierarchical matrices for elastostatics and plasticity.http://www.ifb.tugraz.at/juergen/phd_zechner_2012.

1Department of Mechanical Engineering&Aeronautics,University of Patras,Patras,Greece

2Department of Mechanical Engineering,Technical Research Institute,Patras,Greece

3Department of Mechanical Engineering and Aeronautics,University of Patras,26504,Patras,Greece,+0030 2610969442,Email:polyzos@mech.upatras.gr