Massimiliano Margonari
email: margonar@ing.unitn.it

Marc Bonnet
email: bonnet@lms.polytechnique.fr

Fast multipole method applied to elastostatic BEM-FEM coupling

Keywords: BEM-FEM coupling, Fast multipole method, sparse approximate inverse

Fast multipole method applied to the coupling of elastostatic BEM with FEM

M. Margonari, Marc

Introduction

The finite element method (FEM) and the boundary element method (BEM) are important numerical tools for computing the solutions of many engineering problems. FEM is appropriate for very large classes of situations, including e.g. those with heterogeneous or non-linear constitutive properties, or finite deformations. On the other hand, BEM is useful for modelling special situations such as very large or unbounded domains, geometrical singularities (e.g. cracks) or to obtain very accurate results in regions of complicated shape (see e.g. [START_REF] Aliabadi | Boundary element formulations in fracture mechanics[END_REF][START_REF] Bonnet | Boundary integral equation method for solids and fluids[END_REF][START_REF] Frangi | 3D fracture analysis by symmetric Galerkin BEM[END_REF]). Coupling the BEM and the FEM allows to exploit their complementary advantages when the geometrical configuration warrants it.

The topic of BEM-FEM coupling has been studied since a long time, and many such coupled formulations have been proposed and analysed [START_REF] Hsiao | A hybrid coupled finite-boundary element method in elasticity[END_REF]. In particular, since the traditional collocation BEM (CBEM) formulations lead to unsymmetric systems of coupled BEM-FEM equations, a number of investigations have been directed towards either forcing the symmetry of the CBEM-FEM equations (like in e.g. [START_REF] Beer | Programming the boundary element method[END_REF][START_REF] Li | A new method for the coupling of finite element and boundary element discretized subdomains of elastic body[END_REF]), or use a symmetric Galerkin BEM (SGBEM) formulation in order to obtain naturally a symmetric system of BEM-FEM equations (see e.g. [START_REF] Ganguly | Symmetric coupling of multi-zone curved Galerkin boundary elements with finite elements in elasticity[END_REF][START_REF] Haas | Mixed-dimensional, symmetric coupling of FEM and BEM[END_REF][START_REF] Polizzotto | Variational formulations for coupled BE/FE methods in elastostatics[END_REF][START_REF] Zeng | A symmetric variational finite element-boundary integral equation coupling method[END_REF]). The latter approach is well suited to optimally exploit direct solvers.

As the problem size grows, direct solvers applied to coupled BEM-FEM equations become impractical or infeasible with respect to both computing time and storage, even using specific implementation strategies such as out-of-core procedures, mainly because of the fully-populated nature of the BEM matrices, whose build-up computational cost and storage requirement are both of order O(N 2 BEM ), where N BEM denotes the number of degrees of freedom (DOFs) supported by the BEM mesh, not to mention the O(N 3 BEM ) growth of the solution time. To overcome these pitfalls, one needs to resort to iterative solution algorithms for linear systems, together with an acceleration technique for computing the BEM contribution to the residual of the matrix BEM-FEM equation.

Coupled SGBEM-FEM formulations usually lead to governing matrices that are symmetric but not sign-definite. In such cases, iterative solvers do not take advantage of the symmetry (in contrast with e.g. the conjugate-gradient technique applied to positive definite problems). Hence the final symmetry of the coupled problem is not as important as in connection with direct solvers. Since CBEM is simpler and less costly to set-up, a good case can be made for considering the unsymmetric CBEM-FEM approach.

In this article, a simple CBEM-FEM coupled approach leading to a system of equations solved by means of the generalized minimal residual (GMRES) iterative algorithm [START_REF] Greenbaum | Iterative methods for solving linear systems[END_REF][START_REF] Saad | GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems[END_REF] is presented. The BEM part of the calculation is accelerated by means of the Fast Multipole Method (FMM), a method originally introduced by Rokhlin [START_REF] Rokhlin | Rapid solution of integral equations of classical potential theory[END_REF] and further discussed in e.g. [START_REF] Greengard | A new version of the fast multipole method for the Laplace equation in three dimensions[END_REF] and in the recent review article by Nishimura [START_REF] Nishimura | Fast multipole accelerated boundary integral equation methods[END_REF]. When applied to elastostatic BEM, it provides a reduction of both storage requirements and computational cost to O(N BEM ). These improvements make BEM a viable tool (either on a stand-alone basis or coupled with FEM) for large problems. In addition, a preconditioning technique known as the SParse Approximate Inverse (SPAI) technique is implemented for improving the convergence (i.e. reducing the number of iterations) of GMRES. The article is organized as follows. In section 2, CBEM and FEM formulations are outlined and the coupled problem is presented. Then, the FMM treatment of the BEM equations is presented in section 3. The solution technique, and especially the preconditioning strategy, is discussed in section 4. Finally, numerical examples are examined in section 5.

Coupled CBEM-FEM formulation

Consider a solid occupying a three-dimensional region Ω. A coupled BEM-FEM model of the solid (Figure 1) is defined on the basis of a partition Ω = Ω B ∪Ω F , where ∂Ω B (the boundary of Ω B ) and Ω F respectively support boundary element and finite element discretizations. Let S I = ∂Ω B ∩ ∂Ω F denote the BEM-FEM interface, while S B and S F are the remaining surfaces such that ∂Ω B = S I ∪ S B and ∂Ω F = S I ∪ S F . Both subregions Ω B , Ω F are here endowed with linear elastic properties:

σ ij = C ijkℓ u k,ℓ (in Ω B ) σ ij = C F ijkℓ u k,ℓ (in Ω F ) (1)
where u i and σ ij denote the Cartesian components of displacements and stresses and C ijkℓ , C F ijkℓ are the components of the fourth-order elasticity tensors in each region. Besides, homogeneous and isotropic properties are assumed in Ω B , i.e.

C ijkℓ = 2µ ν 1 -2ν δ ij δ kℓ + δ ik δ jℓ + δ iℓ δ jk (2)
(where µ and ν are the shear modulus and the Poisson ratio) whereas the C F ijkℓ need not be subjected to such restrictions. The isotropy assumption in Ω B is also not mandatory, but is kept here for simplicity and definiteness. The interface need not be simply connected, and S B or S F may be empty (the case of empty S B corresponds to a FEM region embedded in an infinite medium).

CBEM equations.

The governing equation for Ω B is taken here as the usual Somigliana displacement integral equation

c ij (x)u j (x) + ∂Ω B T ij (x, y)u j (y) dS y - ∂Ω B U ij (x, y)t j (y) dS y = ∂Ω B P i (x, y) dS y (3)
where t i = σ ij n j denote the Cartesian components of traction vector defined in terms of the components n j of the unit normal vector directed away from Ω B .

3

The kernels U ij (x, y) and T ij (x, y) are the components of the Kelvin fundamental displacement and traction, whose expressions are given later in equations ( 14) and [START_REF] Mura | Micromechanics of defects in solids[END_REF]. The last integral with kernel P i (x, y) arises when a known gravitational load is considered. The free-term c ij (x) is c ij (x) = δ ij /2 if ∂Ω B is smooth at the collocation point x, and its value is also known if x is located on an edge or corner of ∂Ω B [START_REF] Guiggiani | A general algorithm for multidimensional cauchy principal value integrals in the boundary element method[END_REF][START_REF] Mantic | A new formula for C-matrix in the Somigliana identity[END_REF]. On introducing isoparametric boundary elements to model geometry and boundary fields, considering in turn all mesh nodes as collocation points, and performing the necessary numerical quadratures [START_REF] Beer | Programming the boundary element method[END_REF][START_REF] Bonnet | Boundary integral equation method for solids and fluids[END_REF] and singular integrations using a direct algorithm [START_REF] Guiggiani | A general algorithm for multidimensional cauchy principal value integrals in the boundary element method[END_REF], equation (3) yields the collocation BEM (CBEM) matrix equation

[H]{u} -[G]{t} = {P} (4)
where [H] and [G] are coefficient matrices, {u} and {t} collect nodal values of boundary displacements and tractions, and {P} contains the contribution of the gravitational load.

FEM equations.

In Ω F the standard set of displacement-based FEM discretized equations is considered, i.e.

[K]{u} = {F } + {T }

where u collects all displacement DOFs in Ω F (including those on S I ); [K] is the elastic stiffness matrix; {F } gathers the nodal generalized forces associated with known loads applied to Ω F (e.g. gravitational body forces, prescribed tractions on S F ); and {T } collects the nodal generalized forces associated with tractions along the BEM-FEM interface, i.e. the values of

{T } = S I t i (y)φ a (y) dS y (1 ≤ i ≤ 3 , a ∈ N I ) (6)
where φ a is the trace on S I of the FEM shape function associated with node a of S I and the set N I collects the numbers of all FEM nodes lying on S I . The stiffness matrix [K] is symmetric, positive semidefinite and sparse; it is positive definite if the boundary conditions on S F do not allow any rigid-body motion of Ω F considered in isolation. Here, a skyline storage of [K] has been used to take advantage of symmetry and sparsity. No attempt at optimizing the storage of [K] through node renumbering has been made, although this would certainly be useful.

Coupled formulation. Finally, the coupled formulation must include relationships between the displacement and tractions on S I associated with Ω B and Ω F .

Here, only the usual perfect bonding condition is considered, i.e.:

u B i (y) = u F i (y) and t B i (y) = t F i (y) (y ∈ S I ) (7)
having conventionally defined the traction vector from both sides of S I in terms of the same unit normal vector, namely that pointing away from Ω B . In the discretized formulations, these conditions can be imposed in a weak sense through the relations

S I [u B i (y) -u F i (y)]φ a (y) dS y = 0 S I [t B i (y) -t F i (y)]φ a (y) dS y = 0 (a ∈ N I ) (8)
allowing the use of FEM and BEM interpolations whose traces on S I are not conforming, whereas the strong coupling conditions on the nodal values

u B i (y a ) = u F i (y a ) , t B i (y a ) = t F i (y a ) (a ∈ N I ) (9)
require conforming interpolations on S I (interpolations are said to be conforming on S I if the BEM mesh and the trace on S I of the FEM mesh (i) coincide, and (ii) are associated with the same interpolation functions). In this article, the strong coupling conditions (9) are adopted, but weak coupling of type ( 8) could be easily considered as well.

Gathering equations ( 4), ( 5) and ( 9) leads to the coupled BEM-FEM set of linear equations

       A BB -G BI H BI 0 A IB -G II H II 0 0 D II K II K FI 0 0 K IF K FF                     z B t I u I u F              =              B B B I F I F F              (10)
having introducing partitions of the DOFs so that {t I , u I } gathers the unknown traction and displacement DOFs on S I , {z B } collects all unknown traction and displacement DOFs on S B (i.e. those not prescribed by the boundary conditions), and {u F } gathers the displacements at all nodes of Ω F not lying on S I . The number N BEM of unknown DOFs on the BEM mesh (including the BEM-FEM interface) is the cumulated length of {z B , t I , u I }, while the number N F of internal DOFs on the FEM mesh is the length of {u F }. The first two rows of [START_REF] Greengard | A new version of the fast multipole method for the Laplace equation in three dimensions[END_REF] are obtained by (i) partitioning the CBEM equation (4) according to whether the collocation point x belongs to S B or S I , and (ii) performing the usual column-switching according to the boundary conditions on S B , with the right-hand sides B B , B I gathering all contributions from prescribed data on S B . The last two rows of (10) correspond to the subsets of FEM equations obtained by taking the shape function associated with nodes on S I and inside Ω F , respectively, as trial functions. In addition, the strong coupling assumption implies that tractions on S I are modelled using the shape functions φ a associated with the BEM discretization of S I , so that the load vector {T } in (5) becomes (in terms of the above partition) {T } = [D II ]{t I }, where the entries of the matrix [D II ] are

[D II ] = δ ij S I φ a (y)φ b (y) dS y (1 ≤ i, j ≤ 3 , (a, b) ∈ N I × N I ) (11)
Solution strategy. In the system [START_REF] Greengard | A new version of the fast multipole method for the Laplace equation in three dimensions[END_REF], the blocks generated by the CBEM equations are fully populated and nonsymmetric, whereas the blocks contributed by the FEM equations are sparse and symmetric. For problems of moderate size, direct solvers may be applied to the system [START_REF] Greengard | A new version of the fast multipole method for the Laplace equation in three dimensions[END_REF]. This usually entails some sort of condensation, in order not to build explicitly the whole matrix (which features blocks of zeros). Alternative BEM-FEM formulations based on the symmetric Galerkin BEM instead of the CBEM lead to symmetric variants of the system [START_REF] Greengard | A new version of the fast multipole method for the Laplace equation in three dimensions[END_REF].

As the problem size grows, direct solvers become impractical or infeasible with respect to both computing time and storage, mainly due to the fully-populated nature of the BEM blocks [G] and [H], and iterative solvers are used instead. In the present case of non-symmetric coupled equations, the system (10) is solved using the generalized minimal residual (GMRES) algorithm [START_REF] Greenbaum | Iterative methods for solving linear systems[END_REF][START_REF] Saad | GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems[END_REF], which is applicable to general invertible square matrices. Such algorithms are based on matrix-vector evaluations, and therefore do not require actual storage of the matrix in [START_REF] Greengard | A new version of the fast multipole method for the Laplace equation in three dimensions[END_REF].

In this context, traditional quadrature methods mean that each evaluation of the BEM contribution to the matrix-vector product in [START_REF] Greengard | A new version of the fast multipole method for the Laplace equation in three dimensions[END_REF], i.e. of the residuals

[A BB ] {z B } -[G BI ]{t I } + [H BI ]{u I } [A IB ] {z B } -[G II ]{t I } + [H II ]{u I } (12)
is a O(N 2 BEM ) task, and BEM numerical quadratures are usually time-consuming (i.e. the coefficient of O(N 2 BEM ) is expected to be fairly large). To accelerate the computation of the residuals [START_REF] Guiggiani | A general algorithm for multidimensional cauchy principal value integrals in the boundary element method[END_REF], the Fast Multipole Method (FMM) is adopted. In the present context of elastostatics, choosing this methodology allows to reduce the computational burden of this task to O(N BEM ). Detailed descriptions of FMMbased algorithms can be found in [START_REF] Greengard | A fast algorithm for particle simulations[END_REF][START_REF] Greengard | A new version of the fast multipole method for the Laplace equation in three dimensions[END_REF][START_REF] Nishimura | Fast multipole accelerated boundary integral equation methods[END_REF]. Besides, the FEM stiffness matrix [K] is evaluated only once, stored (in sky-line fashion) and then invoked at every GMRES iteration. The computation of the sparse matrix [D II ] is very inexpensive and is therefore done for every GMRES iteration (i.e. [D II ] is not stored).

The methodology and formulae pertaining to FMM applied to the Somigliana integral equation ( 3) are concisely reviewed in section 3. Then, section 4 will address the preconditioning strategy used in connection with the GMRES algorithm.

Fast Multipole Method (FMM)

Multipole expansions for 3D elasticity

The FMM used in this article follows closely the treatment presented in [START_REF] Yoshida | Application of fast multipole Galerkin boundary integral equation method to elastostatic crack problems in 3D[END_REF] for solving the hypersingular CBEM for elastostatic crack problems. The latter is based on the following series expressing the inverse of the distance r = |y -x| between two points y and x, expanded about two poles y 0 and x 0 :

1 r = +∞ n=0 n m=-n R n,m (ŷ) +∞ n ′ =0 n ′ m ′ =-n ′ (-1) n S n+n ′ ,m+m ′ (r 0 )R n ′ ,m ′ (x) (13)
having put ŷ = yy 0 , x = xx 0 and r 0 = y 0x 0 , and where the overbar indicates complex conjugation. The (complex-valued) solid harmonics R n,m (z), S n,m (z) can be evaluated by means of the recursive formulae (A.1) to (A.5) given in the appendix. For fixed poles y 0 and x 0 , expansion ( 13) is convergent for all y and x such that |ŷ| < |xy 0 | and |x| < |yx 0 | (Figure 2).

To take advantage of expansion [START_REF] Haas | Mixed-dimensional, symmetric coupling of FEM and BEM[END_REF], it is convenient to formulate the elastostatic Kelvin fundamental displacement in the form

U ij (x, y) = 1 2µA (3 -4ν)δ ij 1 r -(x j -r 0 j -ŷj ) ∂ ∂x i 1 r (14)
where A = 8π(1ν). Introducing (13) into the above formula yields the series expansion

U ij (x, y) = 1 2µA +∞ n=0 n m=-n F m,n ij (x-r 0 )R n,m (ŷ)+G m,n i (x-r 0 )ŷ j R n,m (ŷ) (15)
having put

F m,n ij (x -r 0 ) = +∞ n ′ =0 n ′ m ′ =-n ′ (-1) n S n+n ′ ,m+m ′ (r 0 ) × (3 -4ν)δ ij R n ′ ,m ′ (x) -(x j -r 0 j ) ∂ ∂x i R n ′ ,m ′ (x) G m,n i (x -r 0 ) = +∞ n ′ =0 n ′ m ′ =-n ′ (-1) n S n+n ′ ,m+m ′ (r 0 ) ∂ ∂x i R n ′ ,m ′ (x) (16)
Then, the Kelvin traction vector T ij (x, y) is defined by

T ij (x, y) = C jpkℓ n p (y) ∂ ∂y ℓ U ik (x, y) (17)
and therefore admits the series representation

T ij (x, y) = 1 2µA C jpkℓ n p (y) +∞ n=0 n m=-n F m,n ik (x -r 0 ) ∂ ∂y ℓ R n,m (ŷ) + G m,n i (x -r 0 ) ∂ ∂y ℓ ŷk R n,m (ŷ) (18)
The iterative GMRES solver requires repeated evaluations of the residual of the displacement integral equation (3) for known trial solutions (u j , t j ). With that purpose in mind, let the surface ∂Ω B be split into two complementary parts S(x 0 ) and S(x 0 ), where S(x 0 ) contains x 0 as an interior point. The residual of equation (3) for a trial solution (u j , t j ) is recast in the form

c ij (x)u j (x) + S(x 0 ) T ij (x, y)u j (y) dS y - S(x 0 ) U ij (x, y)t j (y) dS y + Īi (x) (19)
having omitted for simplicity the body force term and where Īi (x), defined by

Īi (x) = S(x 0 ) U ij (x, y)t j (y) dS y - S(x 0 ) T ij (x, y)u j (y) dS y ( 20
)
collects the far-field contributions relative to x or x 0 . With respect to decomposition [START_REF] Polizzotto | Variational formulations for coupled BE/FE methods in elastostatics[END_REF], near-field element integrals (i.e. those for elements lying in S(x 0 )) are to be evaluated by means of the usual numerical integration techniques classically used for BEMs [START_REF] Beer | Programming the boundary element method[END_REF][START_REF] Bonnet | Boundary integral equation method for solids and fluids[END_REF], while the far-field contribution Īi (x) will be evaluated with the help of the multipole expansions. Accordingly, decomposition ( 19) is only considered for collocation points x lying inside S(x 0 ) in such a way that S(x 0 ) includes the support of all singular element integrals. The direct algorithm [START_REF] Guiggiani | A general algorithm for multidimensional cauchy principal value integrals in the boundary element method[END_REF] is chosen here for the evaluation of the latter, because decomposition ( 19) is not well suited to using an indirect regularization based on a rigid-body identity.

On introducing the series representations ( 15) and ( 18) of the fundamental kernels into (20), Īi (x) takes the form

Īi (x) = 1 2µA +∞ n=0 n m=-n F m,n ij (x -r 0 ) M t1 j;n,m (y 0 ) -M u1 j;n,m (y 0 ) + G m,n i (x -r 0 ) M t2 n,m (y 0 ) -M u2 n,m (y 0 ) (21)
where the multipole moments are defined by

M t1 j;n,m (y 0 ) = S(x 0 ) R n,m (ŷ)t j (y) dS y M t2 n,m (y 0 ) = S(x 0 ) ŷj R n,m (ŷ)t j (y) dS y M u1 j;n,m (y 0 ) = C kpjℓ S(x 0 ) ∂ ∂y ℓ R n,m (ŷ)u k (y)n p (y) dS y M u2 n,m (y 0 ) = C kpjℓ S(x 0 ) ∂ ∂y ℓ ŷj R n,m (ŷ) u k (y)n p (y) dS y (22)
Further, introducing expressions ( 16) into [START_REF] Saad | Iterative methods for sparse linear systems[END_REF], rearranging summations by switching dummy summation indices (m, n) and (m ′ , n ′ ), Īi (x) takes the form

Īi (x) = 1 2µA +∞ n=0 n m=-n (3-4ν)δ ij R n,m (x)-(x j -r 0 j ) ∂ ∂x i R n,m (x) L 1 j;n,m (x 0 ) + ∂ ∂x i R n,m (x)L 2 j;n,m (x 0 ) (23)
where the local expansions are related to the multipole moments by the following "multipole-to-local" (M2L) relations:

L 1 j;n,m (x 0 ) = +∞ n ′ =0 n ′ m ′ =-n ′ (-1) n ′ S n+n ′ ,m+m ′ (r 0 ) M t1 j;n ′ ,m ′ -M u1 j;n ′ ,m ′ (y 0 ) L 2 n,m (x 0 ) = +∞ n ′ =0 n ′ m ′ =-n ′ S n+n ′ ,m+m ′ (r 0 ) M t2 n ′ ,m ′ -M u2 n ′ ,m ′ (y 0 ) (y 0 ) (24)

Outline of FMM algorithm

Traditional methods for evaluating the residual of integral equation ( 3) require that all element integrals be recomputed each time a new collocation point x is considered, hence the O(N 2 BEM ) overall complexity of such operation. In contrast, equation [START_REF] Yoshida | Application of fast multipole method to boundary integral equation method[END_REF] shows that the same set of local expansions L 1 j;n,m (x 0 ), L 2 j;n,m (x 0 ) allows the treatment of not one but a cluster of collocation points lying sufficiently close to a pole x 0 .

If there is a positive number R such that x < R and |r 0 + ŷ| > 2R, the truncation error estimate

1 r - p n=0 n m=-n R n,m (ŷ) n n ′ =0 n ′ m ′ =-n ′ S n+n ′ ,m+m ′ (r 0 )R n ′ ,m ′ (x) < 1 2 p+1 R (25)
for the multipole expansion (13) [START_REF] Greengard | A new version of the fast multipole method for the Laplace equation in three dimensions[END_REF] implies that the larger the radius of the cluster of collocation points, the smaller the size of the far-field surface S(x 0 ).

To exploit optimally the acceleration afforded by ( 23), a hierarchical oct-tree structure of elements is introduced. For that purpose, a cube containing the whole boundary ∂Ω B , called 'level-0 cell', is divided into eight cubes (level-1 cells), each of which is divided in the same fashion. A level-ℓ cell is divided into level-(ℓ+1) cells unless it contains less than a preset (relatively small) number M of boundary elements (such cells are termed leaves). The FMM algorithm then consists of:

• An upward pass where multipole moments [START_REF] Saad | GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems[END_REF] are first computed for the lowestlevel cells and then recursively aggregated by moving upward in the tree until level 2 (for which there are 4 × 4 × 4 cells overall) is reached. This operation requires the so-called "multipole-to-multipole" (M2M) identities

M (u,t)1 j;n,m (y 0 ) = n n ′ =0 n ′ m ′ =-n ′ R n ′ ,m ′ (y 0′ -y 0 )M (u,t)1 j;n-n ′ ,m-m ′ (y 0′ ) M (u,t)2 n,m (y 0 ) = n n ′ =0 n ′ m ′ =-n ′ R n ′ ,m ′ (y 0′ -y 0 ) M (u,t)2 n-n ′ ,m-m ′ (y 0′ ) -(y 0 j ′ -y 0 j )M (u,t)1 j;n-n ′ ,m-m ′ (y 0′ ) (26)
to shift the origin from the center y 0′ of a level-(ℓ+1) cell to the center y 0 of a level-ℓ cell. • A downward pass where local expansions are first computed at level ℓ = 2 and then evaluated at selected lower-level cells by tracing the tree structure downwards. This operation requires the so-called "local-to-local" (L2L) identities

L 1 j;n,m (x 0′ ) = n n ′ =0 n ′ m ′ =-n ′ R n-n ′ ,m-m ′ (x 0′ -x 0 )L 1 j;n ′ ,m ′ (x 0 ) L 2 n,m (x 0′ ) = n n ′ =0 n ′ m ′ =-n ′ R n-n ′ ,m-m ′ (x 0′ -x 0 ) L 2 n ′ ,m ′ (x 0 ) -(x 0 j ′ -x 0 j )L 1 j;n ′ ,m ′ (x 0 ) (27)
to shift the origin from the center x 0 of a level-(ℓ) cell to the center x 0′ of a level-(ℓ+1) cell. At all levels, only interaction between well-separated cells are so evaluated. • A direct calculation, which includes all singular integrations, where the near-field contributions to [START_REF] Polizzotto | Variational formulations for coupled BE/FE methods in elastostatics[END_REF] are evaluated using conventional integration methods. For BEM formulations of static problems, one computation of the residual of integral equation (3) the FMM is found to be a O(N BEM ) task.

The near-field integrations, as they involve all singular and nearly-singular integrations, are a costly (albeit O(N B )) task. For that reason, they were not recomputed for each iteration (i.e. each evaluation of the residual ( 12)) in the present implementation. Instead, that part of the matrix block is precomputed and assembled into a matrix, prior to invoking GMRES. The near-field direct influence matrix thus obtained is sparse, with O(N BEM ) entries, and is of course stored accordingly.

Preconditioning strategy

As previously mentioned, the GMRES iterative algorithm is used for solving the system (10) of CBEM-FEM equations. The convergence rate of iterative solvers depends strongly on the spectral properties of the coefficient matrix. These spectral properties can be improved by means of suitably chosen linear transformations, i.e. preconditioning. Preconditioning iterative solvers is thus an important practical issue, to which a great deal of attention is devoted in the literature [START_REF] Greenbaum | Iterative methods for solving linear systems[END_REF][START_REF] Saad | Iterative methods for sparse linear systems[END_REF].

The linear system in equation ( 10) has the form

[A]{X} = {B} (28)
where the N ×N matrix [A] is invertible but not symmetric. Preconditioning of (28) can take any of the forms

[P L ][A]{X} = [P L ]{B} (left) [A][P R ]{Y } = {B} , {X} = [P R ]{Y } (right) [P L ][A][P R ]{Y } = [P L ]{B} , {X} = [P R ]{Y } (two-sided) (29)
where the left preconditioner [P L ] or the right preconditioner [P R ] might ideally be chosen so that the governing matrix of the preconditioned system of equations is the identity matrix. Of course, this would amount to solve the original system by means of a direct solver, and therefore be of little practical value. A good preconditioning strategy is such that the matrix of the resulting linear system of equations is as close as possible to the identity matrix. This objective must in addition be fulfilled within computing time and storage requirements significantly lower, and growing slower with N, than those expected for the original system of equations. Preconditioners are either implicit or explicit. The former require the solution of a linear system for each step of the iterative procedure. The most popular are based on an incomplete LU decomposition (ILU or ILUD) of the original matrix or a sparsified version thereof, whereas the latter are based on an explicit evaluation of the preconditioning matrices [P L ] and/or [P R ]. Such matrices are then usually chosen to be sparse so as to reduce the storage requirement and the computational complexity of the preconditioner evaluation.

In this article, a two-sided preconditioning of equation [START_REF] Greengard | A new version of the fast multipole method for the Laplace equation in three dimensions[END_REF] has been adopted, with the left and right preconditioning matrices defined by

P L =         S 0 0 0 0 Id 0 0 0 0 C -T         P R =        D(A BB ) 0 0 0 0 D(-G II ) 0 0 0 0 D(K II ) 0 0 0 0 C -1        (30)
where D(X) denotes the inverse of the main diagonal of X, Id is the identity block and K FF = CC T denotes the Choleski decomposition of K FF . The block S is constructed by computing, in a manner explained later in this section, the sparse approximate inverse (SPAI) of the

N BEM × N BEM matrix      A BB -G BI H BI A IB -G II H II 0 0 Diag(K II )      (31)
and keeping only the first two rows (in the above block notation sense) of the result.

The blocks C and C T clearly define an implicit preconditioning on {u F }, while the preconditioning acting on {z B , t I , u I } is of explicit type. As mentioned above, the block S in (30) results from the computation of a sparse approximate inverse (SPAI). The concept of SPAI is described in e.g. [START_REF] Benzi | A sparse approximate inverse preconditioner for nonsymmetric linear system[END_REF][START_REF] Grote | Parallel preconditioning with sparse approximate inverses[END_REF][START_REF] Saad | Iterative methods for sparse linear systems[END_REF], where in particular theoretical properties and invertibility of such preconditioning matrices are discussed. This technique appears to be effective even when applied to unsymmetric matrices and with an imposed sparsity of 97-99%. An application of the SPAI preconditioning technique to a BEM formulation for electromagnetic wave problems is presented in [START_REF] Alleon | Sparse approximate inverse preconditioning for dense linear system arising in computational electromagnetics[END_REF].

The left sparse approximate inverse [A ♯ ] of a generic square matrix [A] is defined as the solution of the minimization problem

[A ♯ ] = arg min [E] [I] -[E][A] 2 F [E] ∈ R N ×N sparse (32)
where [E] F = (E ab E ab ) 1/2 denotes the Frobenius norm and the sparsity constraint may either take the form of a predetermined sparsity pattern or be found a posteriori so as to satisfy some supplementary optimality criterion. The former, obviously simpler, approach is adopted here. The choice of the Frobenius norm allows to decompose (32) into N independent minimization problems for the rows of [A ♯ ]:

{A ♯ k } = arg min {E} {e k } -{E}[A] {E} ∈ R 1×N sparse ( 33
)
where {e k } and {A ♯ k } are the k th row of the identity matrix and [A ♯ ] respectively, and the trial row N-vector {E} is constrained by the preset sparsity pattern assigned to the k-th row of [A ♯ ].

The effectiveness of SPAI-based preconditioners largely depends on the choice of sparsity pattern for [A ♯ ]. This choice should be such that [A ♯ ] provides a good approximation of the largest entries in [A -1 ]. A difficulty arises in that the location of such entries in [A -1 ] is usually not known a priori. In [START_REF] Grote | Parallel preconditioning with sparse approximate inverses[END_REF] the sparsity pattern is found for every row (or column) with an adaptive strategy, and the positions of non-zero entries in [A ♯ ] are iteratively modified to minimize the residual in (33). In [START_REF] Alleon | Sparse approximate inverse preconditioning for dense linear system arising in computational electromagnetics[END_REF], four heuristics for selecting a priori the sparsity pattern of [A ♯ ] are considered.

Here, the sparsity pattern is predetermined by means of the following strategy. The number m ≪ N of nonzero entries in each row of As previously mentioned, the technique outlined above is used for computing the sparse approximate inverse of the matrix defined by (31), from which the block S of (30) is then extracted.

[A ♯ ] is chosen. A sparse approximation [ Â] of [A]
{ Â♯ i } = arg min { Ê}∈R 1,m {E}[ Ãi ] 2 -2 trace({E}[ Ãi ]) + 1 (1 ≤ i ≤ N) (34)

Numerical examples

Two numerical examples based on the CBEM-FEM coupling method described above are presented in this section. The first one (section 5.1), which concerns the response of an infinite medium to a uniform temperature applied over an ellipsoidal region, has an exact solution. The second one (section 5.2) concerns a real dam structure, in order to demonstrate the present coupling technique on a realistic problem. For both examples, the stopping criterion for the GMRES algorithm was a backward error less than 10 -5 , i.e.: {B} -[A]{X} / {B} ≤ 10 -5 with the notations of equation [START_REF] Greengard | A new version of the fast multipole method for the Laplace equation in three dimensions[END_REF]. All computations were run on a Linux PC computer equipped with one 3 GHz Pentium 4 CPU unit and 2 GBytes of core memory.

Eshelby solution

This first numerical example involves an infinite elastic body containing an ellipsoidal inclusion E, with principal half axes a 1 , a 2 , a 3 (Figure 4). The Cartesian coordinate system is chosen such that the origin coincides with the center of E while the coordinate axes are directed along the principal axes of E. The inclusion and the surrounding medium have the same elastic moduli (µ, ν); in addition, E is endowed with a constant thermal expansion coefficient α. The elastic response of the infinite body to a uniform temperature increment ∆T ⋆ applied over E is such that the stress inside the inclusion is constant. The exact solution for the stress inside E, a special case of the well-known Eshelby solution (see e.g. [START_REF] Mura | Micromechanics of defects in solids[END_REF]), is known. The nonzero components of the stress inside E are given by

σ ii = E 1 -ν (I i -1)α∆T ⋆ (i = 1, 2, 3) (35)
(no summation over i) with

I i = a 1 a 2 a 3 2 ∞ 0 ds (a 2 i + s)∆(s) and ∆(s) = (a 2 1 + s) 1/2 (a 2 2 + s) 1/2 (a 2 3 + s) 1/2
. This example is used here to demonstrate some of the computational features of the present BEM-FEM coupling method. The FEM domain is taken as the inclusion, i.e. Ω F = E, and the BEM-FEM interface is therefore ∂E, the boundary of the inclusion. Although in principle this example could be treated by taking advantage of the symmetries with respect to the three coordinate planes, the entire domain has to about 80 000 (mesh 7). Details on the meshes, oct-tree parameters and numbers of DOFs are provided in Table 1.

In Figure 6 the non-zero stress components (σ 11 , σ 22 and σ 33 evaluated inside E and along the x-axis are plotted against the exact solution (35). The agreement is reasonable, with relative errors always less than 3%. The comparatively high errors recorded near the interface even with fine meshes have to be ascribed to the FEM and not to the coupling technique or to the FMM. The same results have been in fact obtained with the traditional CBEM approach.

Table 2 displays for each of the seven meshes the computer times spent in (i) the computation of the preconditioning blocks associated with the boundary element and finite element parts of the system (10), (ii) the BEM matrix-vector products performed using FMM for each iteration, and (iii) the overall analysis. In particular, the times recorded for setting up the BEM-based part of the preconditioner include the computation of the SPAI as explained in section 4. The term "Cycle" refers to the overall time spent in one BEM matrix-vector products, while "Direct" stands for the computer time spent in reading the sparse direct influence matrix (see end of section 3) and performing the matrix-vector product with the trial solution. The depth of the oct-tree structure was such that a leaf should contain at least M = 25 elements. Setting up the preconditioners appears to require about 25% of the overall analysis computing time for all meshes, largely because computing the SPAI of ( 31) is a relatively time-consuming operation (note in passing that the recorded SPAI computation times are roughly O(N), as expected). The overall analysis time T (N) also appears to increase in a roughly linear fashion with N. In addition, T (N) obtained using either conventional or FMM-based FEM-BEM coupling is plotted against N in Figure 7. The conventional FEM-BEM coupling entails a O(N 2 BEM ) storage requirement which prevented to actually perform the analyses for meshes 5, 6 and 7. Hence, this part of the graph, materialized with a dashed line, is the O(N 2 BEM ) extrapolation of the values of T (N) recorded for meshes 1 to 4. A cross-over point appears around N c/o ≈ 3 10 4 . The value of N c/o is however expected to be strongly problem-, mesh-and implementation-dependent. In most cases, the total computing time reported in Table 2 exceeds Precond + Cycle × Iters because of other tasks which are not reported in the table (with the exception of mesh 4, most likely as a consequence of elementary timings being rounded off).

The same problem has been run for the seven meshes using a simpler preconditioning strategy involving a diagonal preconditioning of the BEM equations, i.e. with [P L ] defined by

P L =        Id 0 0 0 0 Id 0 0 0 0 Id 0 0 0 0 C -T       
while the definition (30) for [P R ] is retained. Table 3 displays, for each of the seven meshes, the overall computing time spent and the iteration count. A comparison of tables 2 and 3 reveals that the SPAI-based preconditioning lears to reduced overall analysis times and stable iteration counts, whereas the iteration counts entailed by the simpler preconditioner are higher and slowly increase with the problem size. The sharp increase of "Direct" CPU time between meshes 6 and 7 is likely due to the fact that in the latter case the sparse direct influence matrix could not be held in the core memory.

Finally, to further demonstrate that the FMM-based CBEM-FEM coupling procedure can be applied to large problems, an array of 8 × 8 × 8 = 256 identical spherical inclusions (diameter 1 m, centers placed at the nodes of a regular cubical grid with node intervals of 3 m along the three principal directions) em- bedded in the infinite elastic medium, subjected to uniform thermal load, has been considered. The other physical parameters are the same as above. The BEM-FEM mesh, whose main characteristics are reported in Table 4, features 480 three-noded triangular boundary elements on each inclusion-matrix interface, about 1200 linear four-noded tetrahedral domain elements in each inclusion, and N = 467073 DOFs overall. The many-inclusion aspect of this problem makes it especially well-suited to FMM, which is here used for the evaluation of geometrically well-separated contributions (interaction terms between distinct inclusions). In this example, the truncation parameter p = 5 was used in the FMM expansions and the depth of the oct-tree structures was such that a leaf should contain at least 100 elements. Only m = 25 non-zero terms for every row have been adopted for the construction of the sparsified preconditioner. Figure 8 displays a vertical section of the array together with the displacement magnitude displacement. The recorded computing times and iteration count are presented in Table 5 in the same format previously used for the single-inclusion problem.

Pian Telessio dam

This second example involves the Pian Telessio dam, an arc-gravity concrete structure located in the north-west of Italy. The main characteristics of the dam are: maximum height 80 m, crest lenght 515 m, crest maximal altitude 1919 m above sea level, maximum storage level 1917 m above sea level, reservoir volume 23.5 10 6 m 3 . The structure and a large portion of the surrounding rock foundation around the dam have been subdivided into three zones, as sketched in Figure 9: the dam, the concrete layer, and the rock foundation. The boundary conditions are as follows: zero displacements are imposed on the bottom and lateral surfaces of the rock-foundation subdomain, a hydrostatic pressure is applied on the rear face of the dam to take into account the stored water, and the other external surfaces are traction-free. Finite elements were used for the first two zones, and boundary elements for the third zone. All materials are treated as linear elastic, with the following characteristics: E = 33000 MPa, ν = 0.12 and γ = 2500 Kg/m 3 (dam), E = 30000 MPa, ν = 0.12 and γ = 2500 Kg/m 3 (concrete layer) and E = 20000 MPa, ν = 0.2 (rock). Three CBEM-FEM meshes of increasing sizes have been set up, together with a FEM model run with ABAQUS for comparison purposes. Mesh 1 is based on linear three-noded boundary triangular elements and four-noded tetrahedral domain elements, while meshes 2 and 3 are based on

The iteration counts and CPU times required by the various phases of the FMM and by the overall analyses, recorded for the three CBEM-FEM meshes of Table 6, are gathered in Table 7 using the same arrangement as in Table 2. Again, a sharp increase of "Direct" CPU time between meshes 2 and 3 indicates that sparse direct influence matrix for mesh 3 could not be held in the core memory.

The crest displacements computed using the present CBEM-FEM technique, displayed in Figure 11, are seen to agree well with values obtained using the FEM code ABAQUS. Discrepancies between the two computations have been registered for the u z component, with a maximum difference of approximately 4•10 -5 m when using mesh 3. Obviously, if a better accuracy is required the solution procedure should be stopped for a smaller relative error.

The deformed shapes of a transversal section of the dam (at a crest abscissa of 151 m), depicted in Figure 12, are also in good agreement with the ABAQUS solution.

In Figures 13 and14, the stress components evaluated on the interface between the concrete layer and the dam, along the s path, are reported. Meshes 2 and 3 (based on quadratic elements) yield better results than mesh 1, based on linear tetrahedral elements which are not well suited for highly accurate stress analyses.

Summary

In this paper a coupling technique between the finite element method (FEM) and the collocation boundary element method (CBEM) has been presented. Its main feature is the recourse to the Fast Multipole Method in order to both accelerate the BEM contribution to the overall computation and allow BEM meshes involving large numbers of degrees of freedom. In addition, a preconditioning strategy based on the sparse approximate inverse concept has been implemented. All these issues have been demonstrated and validated through numerical experiments, run on a Linux PC computer, involving up to about 4 10 5 BEM unknowns. The accuracy of the numerical solutions is not influenced by the FMM if certain parameters (multipole truncation parameter p, minimum number M of elements in a leaf cell) are set appropriately.

Several implementation areas are believed to have potential for improvement of the overall computational performance. The heuristic for selecting the sparsity pattern in the SPAI algorithm was a rather simple one, and can certainly be improved upon. The preconditioning of the FEM part (a full Choleski decomposition of the FEM stiffness matrix) certainly can be improved upon, with potential savings in core memory space. The computation and handling of the sparse direct influence matrix also has scope for optimization.

The CBEM-FEM coupling strategy was presented here for overall linear problems. However, the extension of the present iterative algorithm to problems involving constitutive nonlinearity (e.g. plasticity) in the FEM domain should not give rise to major difficulties.

Fig. 1 .

 1 Fig. 1. FEM-BEM coupling: geometry and notation.

Fig. 2 .

 2 Fig. 2. Geometrical notation for the multipole expansion.

Fig. 3 .

 3 Fig. 3. M2M, M2L and L2L translations

  is created by assigning to [ Â] the m largest entries of each row of [A], at the same locations, and setting all other entries of [ Â] to zero. The sparsity pattern of [ Â] thus defined is then assigned to [A ♯ ] as well, and the entries of [A ♯ ] are found by solving the minimization problems (33) with [A] replaced with [ Â]. Denoting by I(k) the set of the m column indices of nonzero entries in the k-th row of [A ♯ ] or [ Â], the independent minimization prolems are found to reduce to the minimization of quadratic functionals over m-vectors:

Fig. 4 .

 4 Fig. 4. Inclusion problem: geometry and notation. where { Â♯ i } is the m-vector of the nonzero entries of row i of [A ♯ ] and [ Ãi ] is the m×m matrix with entries [A ab ] a∈I(i),b∈I(a) . The N minimization problems (34) are of small size, allowing for a relatively inexpensive set-up of [A ♯ ]. Computing [A ♯ ] is therefore a O(m 3 N) task, while storing [A ♯ ] requires a O(mN) ≪ N 2 memory space. They are well-posed provided the diagonal of [A] is included in the sparsity pattern of [A ♯ ], i.e. that k ∈ I(k) for every k.As previously mentioned, the technique outlined above is used for computing the sparse approximate inverse of the matrix defined by (31), from which the block S of (30) is then extracted.

Fig. 5 .Fig. 6 .

 56 Fig. 5. Inclusion problem: longitudinal section ((x, z)-plane) of mesh 1.

Fig. 7 .

 7 Fig. 7. Inclusion problem: analysis time T (N ) against problem size N for the traditional and FM-M-based CBEM-FEM coupling methods.

Fig. 8 .Fig. 9 .

 89 Fig. 8. Many-inclusion problem: displacement magnitude drawn on a vertical section of the array of spherical inclusions

Fig. 11 .Fig. 12 .Fig. 13 . 21 -Fig. 14 .

 1112132114 Fig. 11. Pian Telessio dam: crest displacement Ux (a), Uy (b) and Uz (c)

Table 1

 1 Inclusion problem: meshes, oct-tree parameters and numbers of DOFs.

	Mesh Nodes	Elements		Oct-tree		DOFs
			BEM	FEM	Max level Leaves N BEM	N F	N = N BEM +N F
	1	267	346	979	3	42	1050	276	1326
	2	822	1038	3153	3	100	3126	903	4029
	3	1362	1540	5563	3	103	4632	1770	6402
	4	2274	2418	9626	4	301	7266	3189	10455
	5	5881	5200	26602	4	422	15612 9837	25449
	6	12868 9402	61770	5	1175	28218 24495	52713
	7	20258 12842 100200	6	1403	38538 41505	80043

Table 2

 2 Inclusion problem: CPU timing using SPAI-based preconditioning.

	Mesh Precond. (s)		Time (s)		Iters Total time
		BEM FEM Upw. Downw. Direct Cycle	n	(s)
	1	10	<1	1	<1	<1	1	37	43
	2	36	<1	1	1	<1	2	37	154
	3	50	<1	1	2	<1	3	37	202
	4	64	3	3	2	1	6	36	277
	5	169	18	7	3	1	11	37	721
	6	349	101	12	6	2	19	38	1425
	7	512	279	18	6	18	42	38	2913

Table 3

 3 Inclusion problem: CPU timing using diagonal preconditioning of BEM equations.

	Mesh	1	2	3	4	5	6	7
	Iters	48 51	53	55	57	58	60
	Total time (s) 63 156 208 315 835 1715 3426

Table 4

 4 Many-inclusion problem: discretization parameters.

	Nodes	Elements		Oct-tree			DOFs
		BEM	FEM	Max. level Leaves N BEM	N F	N = N BEM + N F
	93227 122880 326493	5		7176	374784 92289	467073
	Table 5
	Many-inclusion problem: CPU timing.
		Precond. (s)		Time (s)		Iters Total time
		BEM FEM Upw. Downw. Direct Cycle	n	(s)
		6609	19	47	48	84	180	147	39656

Table 6 Pian

 6 Telessio dam: discretization parameters for the three CBEM-FEM meshes (1,2,3) and the ABAQUS FEM comparison model (A).

	Mesh Nodes		Elements	Oct-tree			DOFs
			BEM	FEM	Max. level Leaves N BEM	N F	Total
	1	25443 21684 (T3) 73569 (T4)	9	8953	38118 43797	81915
	2	23433	7726 (T6) 10307 (T10)	8	3548	50490 46773	97263
	3	51978 15296 (T6) 14462 (T10)	8	6786	96636 64152 160788
	A	406035	-	279742	-	-	-	-	1218105

Table 7 Pian

 7 Telessio dam: CPU timing for the three CBEM-FEM meshes (1,2,3) and the ABAQUS FEM comparison model (A).

	Mesh Precond. (s)		Time (s)		Iters Total time
		BEs FEs Upw. Down. Direct Cycle	n	(s)
	1	186	26	23	27	24	76	83	7916
	2	328	114	11	21	23	57	82	5818
	3	1215 223	23	36	102	165	85	17775
	A	-	-	-	-	-	-	-	3749

 [START_REF] Aliabadi | Boundary element formulations in fracture mechanics[END_REF]This author was supported in the framework of a research project funded by MIUR (Cofin 2002)

The truncation parameter in multipole expansions and the maximum number of elements in a leaf have been set to p = 10 and M = 30, respectively. GMRES was run without restarts. The SPAI preconditioning technique presented in section 4 has been adopted, the sparsified matrix being build with m = 50 non-zero coefficients on each row. The complete CC t factorization has been chosen as left preconditioner for the finite element coefficient matrix, as explained in the previous section.

A Fast Multipole Method: useful formulae

This appendix provides the formulae for practical computation of the solid harmonics R n,m and S n,m , and of the derivatives of the R n,m , after references [START_REF] Yoshida | Application of fast multipole method to boundary integral equation method[END_REF][START_REF] Yoshida | Application of fast multipole Galerkin boundary integral equation method to elastostatic crack problems in 3D[END_REF] to which the reader is referred for more details. First, the evaluation of R n,m (z) and S n,m (z), where z denotes the generic argument of these functions, is performed as follows:

• The R n,m (z) are computed recursively by setting R 0,0 (z) = 1 and using

• The S n,m (z) are computed recursively by setting S 0,0 (z) = 1/ |z| and using

• Finally, R n,m (z) and S n,m (z) for negative values of m are computed via the identities

Identities (A.2) and (A.4) are used with n ≥ m. In addition, the derivatives of the R n,m are computed by means of the following formulae: